diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c3e698e7f..cfc9c4606 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -73,6 +73,9 @@ if(ENABLE_OPENMP) find_package(OpenMP REQUIRED COMPONENTS Fortran) endif() +option(RECOM_COUPLED "Use RECOM" OFF) +message(STATUS "RECOM_COUPLED: ${RECOM_COUPLED}") + option(USE_ICEPACK "Use ICEPACK" OFF) message(STATUS "USE_ICEPACK: ${USE_ICEPACK}") @@ -234,7 +237,9 @@ if(OPENMP_REPRODUCIBLE) endif() if(${RECOM_COUPLED}) - target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2 __3Zoo2Det __coccos)# __usetp) +# target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2 __usetp) + target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2) +# target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2 __3Zoo2Det __coccos __usetp) endif() if(${CISO_COUPLED}) diff --git a/src/MOD_PARTIT.F90 b/src/MOD_PARTIT.F90 index 6603b8922..29dc44327 100644 --- a/src/MOD_PARTIT.F90 +++ b/src/MOD_PARTIT.F90 @@ -71,6 +71,14 @@ module MOD_PARTIT integer :: MPI_COMM_FESOM ! FESOM communicator (for ocean only runs if often a copy of MPI_COMM_WORLD) integer :: MPI_COMM_WORLD ! FESOM communicator (for ocean only runs if often a copy of MPI_COMM_WORLD) +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 communicator for multi FESOM group loop parallelization + integer :: MPI_COMM_FESOM_WORLD + +! kh 17.11.21 communicator for multi FESOM group loop parallelization + integer :: MPI_COMM_FESOM_SAME_RANK_IN_GROUPS +#endif + ! MPI Datatypes for interface exchange ! Element fields (2D; 2D integer; 3D with nl-1 or nl levels, 1 - 4 values) ! small halo and / or full halo @@ -87,6 +95,11 @@ module MOD_PARTIT integer, allocatable :: s_mpitype_nod3D(:,:,:), r_mpitype_nod3D(:,:,:) integer :: MPIERR + +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 multi FESOM group loop parallelization + integer :: my_fesom_group +#endif !!! remPtr_* are constructed during the runtime and shall not be dumped!!! integer, allocatable :: remPtr_nod2D(:), remList_nod2D(:) diff --git a/src/MOD_TRACER.F90 b/src/MOD_TRACER.F90 index efeeed623..d1b872a18 100644 --- a/src/MOD_TRACER.F90 +++ b/src/MOD_TRACER.F90 @@ -20,6 +20,9 @@ MODULE MOD_TRACER real(kind=WP) :: tra_adv_pv = 1. ! a parameter to be used in horizontal advection (for QR4C it is the fraction of fourth-order contribution in the solution) integer :: AB_order=2 integer :: ID +!___________________________________________________________________________ +! TODO: Make it as a part of namelist.tra +logical :: ltra_diag = .true. ! OG - tra_diag contains procedure WRITE_T_TRACER_DATA @@ -41,7 +44,14 @@ MODULE MOD_TRACER ! compute Tstar = 0.5*( T^(n+1) + T^n) real(kind=WP), allocatable, dimension(:,:,:) :: dvd_trflx_hor, dvd_trflx_ver -!_______________________________________________________________________________ +! in case ltra_diag=.true. --> calculate tracer diags ! OG - tra_diag +real(kind=WP), allocatable :: tra_advhoriz(:,:,:), tra_advvert(:,:,:) +real(kind=WP), allocatable :: tra_diff_part_hor_redi(:,:,:) +real(kind=WP), allocatable :: tra_diff_part_ver_expl(:,:,:) +real(kind=WP), allocatable :: tra_diff_part_ver_redi_expl(:,:,:) +real(kind=WP), allocatable :: tra_diff_part_ver_impl(:,:,:) +real(kind=WP), allocatable :: tra_recom_sms(:,:,:) + ! The fct part real(kind=WP),allocatable,dimension(:,:) :: fct_LO ! Low-order solution real(kind=WP),allocatable,dimension(:,:) :: adv_flux_hor ! Antidif. horiz. contrib. from edges / backup for iterafive fct scheme diff --git a/src/associate_part_ass.h b/src/associate_part_ass.h index d2ee2010d..9d4c7a3c9 100644 --- a/src/associate_part_ass.h +++ b/src/associate_part_ass.h @@ -1,3 +1,7 @@ +!DIR$ if defined(__recom) .AND. defined(__usetp) +MPI_COMM_FESOM_WORLD => partit%MPI_COMM_FESOM_WORLD +MPI_COMM_FESOM_SAME_RANK_IN_GROUPS => partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS +!DIR$ endif MPI_COMM_FESOM => partit%MPI_COMM_FESOM MPI_COMM_FESOM_IB => partit%MPI_COMM_FESOM_IB com_nod2D => partit%com_nod2D @@ -13,9 +17,12 @@ eDim_edge2D => partit%eDim_edge2D pe_status => partit%pe_status elem_full_flag => partit%elem_full_flag MPIERR => partit%MPIERR -MPIERR_IB => partit%MPIERR_IB +MPIERR_IB => partit%MPIERR_IB npes => partit%npes mype => partit%mype +!DIR$ if defined(__recom) .AND. defined(__usetp) +my_fesom_group => my_fesom_group +!DIR$ endif maxPEnum => partit%maxPEnum part => partit%part diff --git a/src/associate_part_def.h b/src/associate_part_def.h index 262780a4a..8baef3a9d 100644 --- a/src/associate_part_def.h +++ b/src/associate_part_def.h @@ -1,3 +1,7 @@ +!DIR$ if defined(__recom) .AND. defined(__usetp) + integer, pointer :: MPI_COMM_FESOM_WORLD + integer, pointer :: MPI_COMM_FESOM_SAME_RANK_IN_GROUPS +!DIR$ endif integer, pointer :: MPI_COMM_FESOM ! FESOM communicator (for ocean only runs if often a copy of MPI_COMM_WORLD) integer, pointer :: MPI_COMM_FESOM_IB ! FESOM communicator copy for icebergs LA: 2023-05-22 type(com_struct), pointer :: com_nod2D @@ -20,6 +24,9 @@ integer, pointer :: MPIERR_IB ! copy for icebergs LA: 2023-05-22 integer, pointer :: npes integer, pointer :: mype +!DIR$ if defined(__recom) .AND. defined(__usetp) + integer, pointer :: my_fesom_group +!DIR$ endif integer, pointer :: maxPEnum integer, dimension(:), pointer :: part diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 0c2ad2bc1..3e402cd19 100644 --- a/src/cpl_driver.F90 +++ b/src/cpl_driver.F90 @@ -14,6 +14,9 @@ module cpl_driver ! use mod_oasis ! oasis module use g_config, only : dt, use_icebergs, lwiso +#if defined(__recom) && defined(__usetp) + use g_config, only : num_fesom_groups +#endif use o_param, only : rad USE MOD_PARTIT implicit none @@ -310,7 +313,12 @@ subroutine node_contours(my_x_corners, my_y_corners, partit, mesh) my_y_corners=my_y_corners/rad end subroutine node_contours - subroutine cpl_oasis3mct_init(partit, localCommunicator ) +#if defined(__recom) && defined(__usetp) + subroutine cpl_oasis3mct_init(partit, localCommunicator, num_fesom_groups) +#else + subroutine cpl_oasis3mct_init(partit, localCommunicator) +#endif + USE MOD_PARTIT implicit none save @@ -324,6 +332,9 @@ subroutine cpl_oasis3mct_init(partit, localCommunicator ) ! integer, intent(OUT) :: localCommunicator type(t_partit), intent(inout), target :: partit +#if defined(__recom) && defined(__usetp) + integer, intent(inout) :: num_fesom_groups +#endif ! ! Local declarations ! @@ -345,7 +356,11 @@ subroutine cpl_oasis3mct_init(partit, localCommunicator ) !------------------------------------------------------------------ ! 1st Initialize the OASIS3-MCT coupling system for the application !------------------------------------------------------------------ +#if defined(__recom) && defined(__usetp) + CALL oasis_init_comp(comp_id, comp_name, ierror, num_program_groups = num_fesom_groups) +#else CALL oasis_init_comp(comp_id, comp_name, ierror ) +#endif IF (ierror /= 0) THEN CALL oasis_abort(comp_id, 'cpl_oasis3mct_init', 'Init_comp failed.') ENDIF @@ -356,7 +371,11 @@ subroutine cpl_oasis3mct_init(partit, localCommunicator ) CALL oasis_abort(comp_id, 'cpl_oasis3mct_init', 'comm_rank failed.') ENDIF +#if defined(__recom) && defined(__usetp) + CALL oasis_get_localcomm_all_groups( localCommunicator, ierror ) +#else CALL oasis_get_localcomm( localCommunicator, ierror ) +#endif IF (ierror /= 0) THEN CALL oasis_abort(comp_id, 'cpl_oasis3mct_init', 'get_local_comm failed.') ENDIF @@ -606,6 +625,10 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh) print *, 'FESOM after Barrier' endif +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif + if (mype .eq. localroot) then print *, 'FESOM before grid writing to oasis grid files' CALL oasis_start_grids_writing(il_flag) @@ -632,6 +655,9 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh) print *, 'FESOM after terminate_grids_writing' endif !localroot +#if defined(__recom) && defined(__usetp) + end if !(partit%my_fesom_group == 0) then +#endif DEALLOCATE(all_x_coords, all_y_coords, my_x_coords, my_y_coords, displs_from_all_pes, counts_from_all_pes) @@ -900,15 +926,49 @@ subroutine cpl_oasis3mct_recv(ind, data_array, action, partit) endif #endif +#if defined(__recom) && defined(__usetp) +! the coupling is in principle as it was before, i.e. the fesom processes - in group 0 - receive their data from echam + if(partit%my_fesom_group == 0) then +#endif + call oasis_get(recv_id(ind), seconds_til_now, exfld,info) + +#if defined(__recom) && defined(__usetp) + else + +! defensive: assignment statement "action=(info==3 ..." below is "don't care" in this case, because the actual value for action +! is received via MPI_Bcast anyway + info = 0 + + end if +#endif + t2=MPI_Wtime() ! ! FESOM's interpolation routine interpolates structured ! VarStrLoc coming from OASIS3MCT to local unstructured data_array ! and delivered back to FESOM. action=(info==3 .OR. info==10 .OR. info==11 .OR. info==12 .OR. info==13) + +#if defined(__recom) && defined(__usetp) + if(num_fesom_groups > 1) then + call MPI_Bcast(action, 1, MPI_LOGICAL, 0, partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, partit%MPIerr) + end if +#endif + if (action) then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif data_array(1:partit%myDim_nod2d) = exfld +#if defined(__recom) && defined(__usetp) + end if + + if(num_fesom_groups > 1) then + call MPI_Bcast(data_array, partit%myDim_nod2d, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, partit%MPIerr) + end if +#endif + call exchange_nod(data_array, partit) end if t3=MPI_Wtime() diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 3302f61f4..097b9e5f9 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -36,6 +36,7 @@ module fesom_main_storage_module use age_tracer_init_interface use iceberg_params use iceberg_step + use mod_transit ! Define icepack module #if defined (__icepack) @@ -60,6 +61,9 @@ module fesom_main_storage_module integer :: which_readr ! read which restart files (0=netcdf, 1=core dump,2=dtype) integer :: total_nsteps integer, pointer :: mype, npes, MPIerr, MPI_COMM_FESOM, MPI_COMM_WORLD, MPI_COMM_FESOM_IB +#if defined(__recom) && defined(__usetp) + integer, pointer :: my_fesom_group, MPI_COMM_FESOM_WORLD, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS +#endif real(kind=WP) :: t0, t1, t2, t3, t4, t5, t6, t7, t8, t0_ice, t1_ice, t0_frc, t1_frc real(kind=WP) :: rtime_fullice, rtime_write_restart, rtime_write_means, rtime_compute_diag, rtime_read_forcing real(kind=real32) :: rtime_setup_mesh, rtime_setup_ocean, rtime_setup_forcing @@ -112,10 +116,25 @@ subroutine fesom_init(fesom_total_nsteps) #if defined(__MULTIO) use iom #endif + use cpl_driver integer, intent(out) :: fesom_total_nsteps ! EO parameters logical mpi_is_initialized integer :: tr_num + +#if defined(__recom) && defined(__usetp) +! multi FESOM group loop parallelization +! moved from fvom_main.F90 + integer :: npes_fesom_world + integer :: mype_fesom_world + integer :: processes_per_group + integer :: npes_check + integer :: mype_check + +! get current value for num_fesom_groups + call read_namelist_run_config +#endif + #if !defined __ifsinterface if(command_argument_count() > 0) then call command_line_options%parse() @@ -148,7 +167,13 @@ subroutine fesom_init(fesom_total_nsteps) #if defined (__oasis) - call cpl_oasis3mct_init(f%partit,f%partit%MPI_COMM_FESOM) +! pass num_fesom_groups to coupler +#if defined(__recom) && defined(__usetp) + call cpl_oasis3mct_init(f%partit, f%partit%MPI_COMM_FESOM, num_fesom_groups) +#else + call cpl_oasis3mct_init(f%partit, f%partit%MPI_COMM_FESOM) +#endif + #endif f%t1 = MPI_Wtime() @@ -162,6 +187,106 @@ subroutine fesom_init(fesom_total_nsteps) f%npes =>f%partit%npes +#if defined(__recom) && defined(__usetp) +! prepare communicator splitting for multi FESOM group loop parallelization + f%my_fesom_group=>f%partit%my_fesom_group + + f%MPI_COMM_FESOM_WORLD=> f%partit%MPI_COMM_FESOM_WORLD + f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS=> f%partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS + + f%MPI_COMM_FESOM_WORLD = f%MPI_COMM_FESOM + npes_fesom_world = f%npes + mype_fesom_world = f%mype + if(mype_fesom_world == 0) then + write(*,*) 'npes_fesom_world, num_fesom_groups', npes_fesom_world, num_fesom_groups + end if + if(mod(npes_fesom_world, num_fesom_groups) /= 0) then + if(mype_fesom_world == 0) then + write(*,*) 'MPI_comm_split mismatch npes_fesom_world, num_fesom_groups', npes_fesom_world, num_fesom_groups + end if + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + processes_per_group = npes_fesom_world / num_fesom_groups + if(mype_fesom_world == 0) then + write(*,*) 'processes_per_group', processes_per_group + end if + f%npes = processes_per_group + f%my_fesom_group = mype_fesom_world / processes_per_group + f%mype = mod(mype_fesom_world, processes_per_group) + +! split to num_fesom_groups + call MPI_comm_split(f%MPI_COMM_FESOM_WORLD, f%my_fesom_group, 0, f%MPI_COMM_FESOM, f%MPIerr) + if (f%MPIerr /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_split(MPI_COMM_FESOM_WORLD, my_fesom_group, 0, MPI_COMM_FESOM, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + call MPI_comm_size(f%MPI_COMM_FESOM, npes_check, f%MPIerr) + if(f%MPIerr /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_size(MPI_COMM_FESOM, npes_check, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + call MPI_comm_rank(f%MPI_COMM_FESOM, mype_check, f%MPIerr) + if(f%MPIerr /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_rank(MPI_COMM_FESOM, mype_check, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + if(npes_check /= f%npes) then + write(*,*) 'npes mismatch, npes, npes_check', f%npes, npes_check + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + if(mype_check /= f%mype) then + write(*,*) 'mype mismatch, mype, mype_check', f%mype, mype_check + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + +! group same ranks in each group for broadcasting + + call MPI_comm_split(f%MPI_COMM_FESOM_WORLD, f%mype, f%my_fesom_group, f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, f%MPIERR) + if (f%MPIERR /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_split(MPI_COMM_FESOM_WORLD, mype, my_fesom_group, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + call MPI_comm_size(f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, npes_check, f%MPIERR) + if(f%MPIERR /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_size(MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, npes_check, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + call MPI_comm_rank(f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, mype_check, f%MPIERR) + if(f%MPIERR /= MPI_SUCCESS) then + write(*,*) 'MPI_comm_rank(MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, mype_check, MPIERR) failed' + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + if(npes_check /= num_fesom_groups) then + write(*,*) 'npes mismatch, num_fesom_groups, npes_check', num_fesom_groups, npes_check + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + if(mype_check /= f%my_fesom_group) then + write(*,*) 'mype mismatch, my_fesom_group, mype_check', f%my_fesom_group, mype_check + call par_ex(f%MPI_COMM_FESOM, f%mype) + stop + end if + + if(f%my_fesom_group==0) then +#endif if(f%mype==0) then write(*,*) @@ -171,20 +296,53 @@ subroutine fesom_init(fesom_total_nsteps) print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' print *, achar(27)//'[7;32m'//' --> FESOM BUILDS UP MODEL CONFIGURATION '//achar(27)//'[0m' end if + +#if defined(__recom) && defined(__usetp) + end if +#endif + !===================== ! Read configuration data, ! load the mesh and fill in ! auxiliary mesh arrays !===================== + if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call setup_model'//achar(27)//'[0m' call setup_model(f%partit) ! Read Namelists, always before clock_init + + if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call clock_init'//achar(27)//'[0m' call clock_init(f%partit) ! read the clock file + + if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call get_run_steps'//achar(27)//'[0m' call get_run_steps(fesom_total_nsteps, f%partit) f%total_nsteps=fesom_total_nsteps + if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call mesh_setup'//achar(27)//'[0m' call mesh_setup(f%partit, f%mesh) +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif + if (f%mype==0) write(*,*) 'FESOM mesh_setup... complete' +#if defined(__recom) && defined(__usetp) + end if +#endif + +! Transient tracers: control output of initial input values + if(use_transit .and. anthro_transit .and. f%mype==0) then + write (*,*) + write (*,*) "*** Transient tracers: Initial atmospheric input values >>>" + write (*,*) "Year CE, xCO2, D14C_NH, D14C_TZ, D14C_SH, xCFC-11_NH, xCFC-11_SH, xCFC-12_NH, xCFC-12_SH, xSF6_NH, xSF6_SH" + write (*, fmt="(2x,i4,10(2x,f6.2))") & + year_ce(ti_transit), xCO2_ti(ti_transit) * 1.e6, & + (r14c_nh(ti_transit) - 1.) * 1000., (r14c_tz(ti_transit) - 1.) * 1000., (r14c_sh(ti_transit) - 1.) * 1000., & + xf11_nh(ti_transit) * 1.e12, xf11_sh(ti_transit) * 1.e12, & + xf12_nh(ti_transit) * 1.e12, xf12_sh(ti_transit) * 1.e12, & + xsf6_nh(ti_transit) * 1.e12, xsf6_sh(ti_transit) * 1.e12 + write (*,*) + end if + !===================== ! Allocate field variables ! and additional arrays needed for @@ -220,15 +378,36 @@ subroutine fesom_init(fesom_total_nsteps) ! recom setup #if defined (__recom) +#if defined (__usetp) + if(f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call recom_init'//achar(27)//'[0m' +#if defined (__usetp) + end if +#endif + f%t0_recom=MPI_Wtime() call recom_init(f%tracers, f%partit, f%mesh) ! adjust values for recom tracers (derived type "t_tracer") f%t1_recom=MPI_Wtime() + +#if defined (__usetp) + if(f%my_fesom_group==0) then +#endif if (f%mype==0) write(*,*) 'RECOM recom_init... complete' +#if defined (__usetp) + end if +#endif #endif if (f%mype==0) then +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif write(*,*) 'FESOM ocean_setup... complete' +#if defined(__recom) && defined(__usetp) + end if +#endif + f%t3=MPI_Wtime() endif call forcing_setup(f%partit, f%mesh) @@ -239,7 +418,15 @@ subroutine fesom_init(fesom_total_nsteps) call ice_setup(f%ice, f%tracers, f%partit, f%mesh) f%ice%ice_steps_since_upd = f%ice%ice_ave_steps-1 f%ice%ice_update=.true. + +#if defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (f%mype==0) write(*,*) 'EVP scheme option=', f%ice%whichEVP +#if defined(__usetp) + end if +#endif + else ! create a dummy ice derived type with only a_ice, m_ice, m_snow and ! uvice since oce_timesteps still needs in moment @@ -264,10 +451,29 @@ subroutine fesom_init(fesom_total_nsteps) !---age-code-end #if defined (__oasis) +! only mype == 0 in my_fesom_group == 0 handles coupling with extern models +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif call cpl_oasis3mct_define_unstr(f%partit, f%mesh) - if(f%mype==0) write(*,*) 'FESOM ----> cpl_oasis3mct_define_unstr nsend, nrecv:',nsend, nrecv +#if defined(__recom) && defined(__usetp) + end if #endif + +#if defined(__recom) && defined(__usetp) + call MPI_Barrier(f%MPI_COMM_FESOM_WORLD, f%MPIERR) + if(num_fesom_groups > 1) then + + call MPI_Bcast(cpl_send, sizeof(cpl_send), MPI_CHARACTER, 0, f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, f%MPIerr) + call MPI_Bcast(cpl_recv, sizeof(cpl_recv), MPI_CHARACTER, 0, f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, f%MPIerr) + +! needed in SUBROUTINE net_rec_from_atm(action) + call MPI_Bcast(target_root, 1, MPI_INTEGER, 0, f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, f%MPIerr) + end if +#endif + +#endif ! defined (__oasis) ! -------------- ! LA icebergs: 2023-05-17 @@ -286,13 +492,21 @@ subroutine fesom_init(fesom_total_nsteps) call init_icepack(f%ice, f%tracers%data(1), f%mesh) if (f%mype==0) write(*,*) 'Icepack: setup complete' #endif + call clock_newyear ! check if it is a new year if (f%mype==0) f%t6=MPI_Wtime() !___CREATE NEW RESTART FILE IF APPLICABLE___________________________________ call restart(0, 0, 0, r_restart, f%which_readr, f%ice, f%dynamics, f%tracers, f%partit, f%mesh) if (f%mype==0) f%t7=MPI_Wtime() ! store grid information into netcdf file + +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (.not. r_restart) call write_mesh_info(f%partit, f%mesh) +#if defined(__recom) && defined(__usetp) + end if +#endif !___IF RESTART WITH ZLEVEL OR ZSTAR IS DONE, ALSO THE ACTUAL LEVELS AND ____ !___MIDDEPTH LEVELS NEEDS TO BE CALCULATET AT RESTART_______________________ @@ -313,6 +527,10 @@ subroutine fesom_init(fesom_total_nsteps) f%rtime_setup_recom = real( f%t1_recom - f%t0_recom ,real32) #endif +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif + write(*,*) '==========================================' write(*,*) 'MODEL SETUP took on mype=0 [seconds] ' write(*,*) 'runtime setup total ',real(f%t8-f%t1,real32) @@ -326,6 +544,11 @@ subroutine fesom_init(fesom_total_nsteps) write(*,*) ' > runtime setup recom ',f%rtime_setup_recom #endif write(*,*) '============================================' + +#if defined(__recom) && defined(__usetp) + end if +#endif + endif #if defined(__MULTIO) @@ -443,24 +666,46 @@ subroutine fesom_runloop(current_nsteps) ! -------------- ! LA icebergs: 2023-05-17 +!YY: only when using icebergs? + if (use_icebergs) then f%MPI_COMM_FESOM_IB = f%MPI_COMM_FESOM if (f%mype==0) then ! write (*,*) 'ib_async_mode, initial omp_num_threads ', ib_async_mode, omp_get_num_threads() write (*,*) 'current_nsteps, steps_per_ib_step, icb_outfreq :', current_nsteps, steps_per_ib_step, icb_outfreq end if + end if ! -------------- - +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (f%mype==0) write(*,*) 'FESOM start iteration before the barrier...' +#if defined(__recom) && defined(__usetp) + end if +#endif call MPI_Barrier(f%MPI_COMM_FESOM, f%MPIERR) if (f%mype==0) then +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif write(*,*) 'FESOM start iteration after the barrier...' +#if defined(__recom) && defined(__usetp) + end if +#endif f%t0 = MPI_Wtime() endif + +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if(f%mype==0) then write(*,*) print *, achar(27)//'[32m' //'____________________________________________________________'//achar(27)//'[0m' print *, achar(27)//'[7;32m'//' --> FESOM STARTS TIME LOOP '//achar(27)//'[0m' end if +#if defined(__recom) && defined(__usetp) + end if +#endif + !___MODEL TIME STEPPING LOOP________________________________________________ if (use_global_tides) then call foreph_ini(yearnew, month, f%partit) @@ -526,18 +771,32 @@ subroutine fesom_runloop(current_nsteps) call foreph(f%partit, f%mesh) end if mstep = n + +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (mod(n,logfile_outfreq)==0 .and. f%mype==0) then write(*,*) 'FESOM =======================================================' ! write(*,*) 'FESOM step:',n,' day:', n*dt/24./3600., write(*,*) 'FESOM step:',n,' day:', daynew,' year:',yearnew write(*,*) end if +#if defined(__recom) && defined(__usetp) + end if +#endif + #if defined (__oifs) || defined (__oasis) seconds_til_now=INT(dt)*(n-1) #endif call clock !___compute horizontal velocity on nodes (originaly on elements)________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call compute_vel_nodes'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif call compute_vel_nodes(f%dynamics, f%partit, f%mesh) ! -------------- ! LA icebergs: 2023-05-17 @@ -551,11 +810,23 @@ subroutine fesom_runloop(current_nsteps) f%t1 = MPI_Wtime() if(use_ice) then !___compute fluxes from ocean to ice________________________________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call ocean2ice(n)'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif call ocean2ice(f%ice, f%dynamics, f%tracers, f%partit, f%mesh) !___compute update of atmospheric forcing____________________________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call update_atm_forcing(n)'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif f%t0_frc = MPI_Wtime() call update_atm_forcing(n, f%ice, f%tracers, f%dynamics, f%partit, f%mesh) f%t1_frc = MPI_Wtime() @@ -567,7 +838,13 @@ subroutine fesom_runloop(current_nsteps) f%ice%ice_update=.false. f%ice%ice_steps_since_upd=f%ice%ice_steps_since_upd+1 endif +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call ice_timestep(n)'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif if (f%ice%ice_update) call ice_timestep(n, f%ice, f%partit, f%mesh) @@ -581,7 +858,13 @@ subroutine fesom_runloop(current_nsteps) !___compute fluxes to the ocean: heat, freshwater, momentum_________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call oce_fluxes_mom...'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif call oce_fluxes_mom(f%ice, f%dynamics, f%partit, f%mesh) ! momentum only call oce_fluxes(f%ice, f%dynamics, f%tracers, f%partit, f%mesh) end if @@ -590,26 +873,51 @@ subroutine fesom_runloop(current_nsteps) !___now recom____________________________________________________ #if defined (__recom) +#if defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (f%mype==0 .and. n==1) print *, achar(27)//'[46' //'_____________________________________________________________'//achar(27)//'[0m' if (f%mype==0 .and. n==1) print *, achar(27)//'[46;1m'//' --> call REcoM '//achar(27)//'[0m' +#if defined(__usetp) + end if +#endif + f%t0_recom = MPI_Wtime() call recom(f%ice, f%dynamics, f%tracers, f%partit, f%mesh) f%t1_recom = MPI_Wtime() #endif !___model ocean step____________________________________________________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call oce_timestep_ale'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif call oce_timestep_ale(n, f%ice, f%dynamics, f%tracers, f%partit, f%mesh) f%t3 = MPI_Wtime() !___compute energy diagnostics..._______________________________________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call compute_diagnostics(1)'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + end if +#endif call compute_diagnostics(1, f%dynamics, f%tracers, f%ice, f%partit, f%mesh) f%t4 = MPI_Wtime() !___prepare output______________________________________________________ +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif if (flag_debug .and. f%mype==0) print *, achar(27)//'[34m'//' --> call output (n)'//achar(27)//'[0m' call output (n, f%ice, f%dynamics, f%tracers, f%partit, f%mesh) +#if defined(__recom) && defined(__usetp) + end if +#endif ! LA icebergs: 2023-05-17 if (use_icebergs .and. mod(n, steps_per_ib_step)==0.0) then @@ -629,6 +937,24 @@ subroutine fesom_runloop(current_nsteps) #if defined (__recom) f%rtime_compute_recom = f%rtime_compute_recom + f%t1_recom - f%t0_recom #endif + +! Transient tracers: update of input values between restarts + if(use_transit .and. anthro_transit .and. (daynew == ndpyr) .and. (timenew==86400.)) then + ti_transit = ti_transit + 1 + if (f%mype==0) then + write (*,*) + write (*,*) "*** Transient tracers: Updated atmospheric input values >>>" + write (*,*) "Year CE, xCO2, D14C_NH, D14C_TZ, D14C_SH, xCFC-11_NH, xCFC-11_SH, xCFC-12_NH, xCFC-12_SH, xSF6_NH, xSF6_SH" + write (*, fmt="(2x,i4,10(2x,f6.2))") & + year_ce(ti_transit), xCO2_ti(ti_transit) * 1.e6, & + (r14c_nh(ti_transit) - 1.) * 1000., (r14c_tz(ti_transit) - 1.) * 1000., (r14c_sh(ti_transit) - 1.) * 1000., & + xf11_nh(ti_transit) * 1.e12, xf11_sh(ti_transit) * 1.e12, & + xf12_nh(ti_transit) * 1.e12, xf12_sh(ti_transit) * 1.e12, & + xsf6_nh(ti_transit) * 1.e12, xsf6_sh(ti_transit) * 1.e12 + write (*,*) + end if + endif + end do !call cray_acc_set_debug_global_level(3) f%from_nstep = f%from_nstep+current_nsteps @@ -636,7 +962,6 @@ subroutine fesom_runloop(current_nsteps) ! write(0,*) 'f%from_nstep after the loop:', f%from_nstep end subroutine - subroutine fesom_finalize() use fesom_main_storage_module #if defined(__MULTIO) @@ -646,18 +971,41 @@ subroutine fesom_finalize() ! EO parameters real(kind=real32) :: mean_rtime(15), max_rtime(15), min_rtime(15) integer :: tr_num + integer :: i !YY ! -------------- ! LA icebergs: 2023-05-17 if (use_icebergs) then call iceberg_out(f%partit) end if ! -------------- +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif call finalize_output() call finalize_restart() +#if defined(__recom) && defined(__usetp) + end if +#endif !___FINISH MODEL RUN________________________________________________________ +#if !defined (__usetp) +! multi FESOM group loop parallelization call MPI_Barrier(f%MPI_COMM_FESOM, f%MPIERR) +#endif +#if defined(__recom) && defined (__usetp) +! list statistics for all fesom_groups +! fesom groups are listed backwards, so info for the main fesom group 0 is at the end in the log + do i = num_fesom_groups - 1, 0, -1 + +! use a barrier to "sort" the output but the mpi output can still get a bit mixed up, +! because MPI does not define the handling of the order of the output lines + call MPI_Barrier(f%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, f%MPIERR) + +! for the sake of output clarity produce output only for my_fesom_group == 0 for now + if(i == f%my_fesom_group .and. f%my_fesom_group == 0) then +#endif + !$ACC EXIT DATA DELETE (f%ice%delta_min, f%ice%Tevp_inv, f%ice%cd_oce_ice) !$ACC EXIT DATA DELETE (f%ice%work%fct_tmax, f%ice%work%fct_tmin) !$ACC EXIT DATA DELETE (f%ice%work%fct_fluxes, f%ice%work%fct_plus, f%ice%work%fct_minus) @@ -745,10 +1093,20 @@ subroutine fesom_finalize() call par_ex(f%partit%MPI_COMM_FESOM, f%partit%mype) #endif +#if defined(__recom) && defined (__usetp) + end if + end do ! i = num_fesom_groups - 1, 0, -1 +#endif + #if defined(__MULTIO) && !defined(__ifsinterface) && !defined(__oasis) call mpp_stop #endif if(f%fesom_did_mpi_init) call par_ex(f%partit%MPI_COMM_FESOM, f%partit%mype) ! finalize MPI before FESOM prints its stats block, otherwise there is sometimes output from other processes from an earlier time in the programm AFTER the starts block (with parastationMPI) + +#if defined(__recom) && defined(__usetp) + if (f%my_fesom_group==0) then +#endif + if (f%mype==0) then 41 format (a35,a10,2a15) !Format for table heading 42 format (a30,3f15.4) !Format for table content @@ -787,6 +1145,11 @@ subroutine fesom_finalize() write(*,*) '======================================================' write(*,*) end if + +#if defined(__recom) && defined(__usetp) + end if +#endif + ! call clock_finish end subroutine diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index be12855ce..0a0bdcbc5 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -269,7 +269,15 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) print *, 'not installed yet or error in cpl_oasis3mct_send', mype #endif ! oifs endif + +! kh 30.11.21 +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif call cpl_oasis3mct_send(i, exchange, action, partit) +#if defined(__recom) && defined(__usetp) + endif +#endif end do #ifdef VERBOSE do i=1, nsend @@ -834,11 +842,20 @@ SUBROUTINE net_rec_from_atm(action, partit) use o_PARAM, only: WP USE MOD_PARTIT USE MOD_PARSUP + +#if defined(__recom) && defined(__usetp) +! kh 10.21.21 + use g_config, only: num_fesom_groups +#endif IMPLICIT NONE LOGICAL, INTENT (IN) :: action type(t_partit), intent(inout), target :: partit INTEGER :: my_global_rank, ierror +#if defined(__recom) && defined(__usetp) +! kh 10.12.21 + INTEGER :: my_global_rank_test +#endif INTEGER :: n INTEGER :: status(MPI_STATUS_SIZE,partit%npes) INTEGER :: request(2) @@ -851,11 +868,30 @@ SUBROUTINE net_rec_from_atm(action, partit) CALL MPI_COMM_RANK(MPI_COMM_WORLD, my_global_rank, ierror) atm_net_fluxes_north=0. atm_net_fluxes_south=0. +#if defined(__recom) && defined(__usetp) +! kh 10.12.21 + my_global_rank_test = my_global_rank - (partit%my_fesom_group * partit%npes) +#endif + +#if defined(__recom) && defined(__usetp) +! kh 10.12.21 check for is root in group + if (my_global_rank_test==target_root) then + if(partit%my_fesom_group == 0) then +#else if (my_global_rank==target_root) then - CALL MPI_IRecv(atm_net_fluxes_north(1), nrecv, MPI_DOUBLE_PRECISION, source_root, 111, MPI_COMM_WORLD, request(1), partit%MPIerr) +#endif + CALL MPI_IRecv(atm_net_fluxes_north(1), nrecv, MPI_DOUBLE_PRECISION, source_root, 111, MPI_COMM_WORLD, request(1), partit%MPIerr) CALL MPI_IRecv(atm_net_fluxes_south(1), nrecv, MPI_DOUBLE_PRECISION, source_root, 112, MPI_COMM_WORLD, request(2), partit%MPIerr) CALL MPI_Waitall(2, request, status, partit%MPIerr) end if + +#if defined(__recom) && defined(__usetp) + if(num_fesom_groups > 1) then + call MPI_Bcast(atm_net_fluxes_north(1), nrecv, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, partit%MPIerr) + call MPI_Bcast(atm_net_fluxes_south(1), nrecv, MPI_DOUBLE_PRECISION, 0, partit%MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, partit%MPIerr) + end if + end if ! (my_global_rank_test==target_root) then +#endif call MPI_Barrier(partit%MPI_COMM_FESOM, partit%MPIerr) call MPI_AllREDUCE(atm_net_fluxes_north(1), aux, nrecv, MPI_DOUBLE_PRECISION, MPI_SUM, partit%MPI_COMM_FESOM, partit%MPIerr) atm_net_fluxes_north=aux diff --git a/src/gen_ic3d.F90 b/src/gen_ic3d.F90 index 7513fec4b..50681b76c 100644 --- a/src/gen_ic3d.F90 +++ b/src/gen_ic3d.F90 @@ -613,7 +613,9 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) call MPI_AllREDUCE(locSmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' `-> gobal min init. salt. =', glo #if defined(__recom) - +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (partit%mype==0) write(*,*) "Sanity check for REcoM variables" call MPI_AllREDUCE(locDINmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' |-> gobal max init. DIN. =', glo @@ -640,6 +642,9 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) if (partit%mype==0) write(*,*) ' |-> gobal max init. O2. =', glo call MPI_AllREDUCE(locO2min , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, partit%MPI_COMM_FESOM, partit%MPIerr) if (partit%mype==0) write(*,*) ' `-> gobal min init. O2. =', glo +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif #endif END SUBROUTINE do_ic3d diff --git a/src/gen_model_setup.F90 b/src/gen_model_setup.F90 index 13a669a3a..deaa3a789 100755 --- a/src/gen_model_setup.F90 +++ b/src/gen_model_setup.F90 @@ -174,6 +174,34 @@ subroutine setup_model(partit) endif ! if ((output_length_unit=='s').or.(int(real(step_per_day)/24.0)<=1)) use_means=.false. end subroutine setup_model + + +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 read num_fesom_groups for multi FESOM group loop parallelization +! ================================================================= +subroutine read_namelist_run_config(partit) + ! Reads run_config namelist and overwrite default parameters. + ! + ! kh 11.11.21 Copied by Kai Himstedt (based on read_namelist) + !-------------------------------------------------------------- + USE MOD_PARTIT + USE MOD_PARSUP + use g_config + implicit none + type(t_partit), intent(inout), target :: partit + + character(len=100) :: nmlfile + integer fileunit + + nmlfile ='namelist.config' ! name of general configuration namelist file + open (newunit=fileunit, file=nmlfile) + + open (fileunit,file=nmlfile) +! read (fileunit,NML=run_config) + read (fileunit,NML=run_config_tp) + close (fileunit) +end subroutine read_namelist_run_config +#endif ! ================================================================= subroutine get_run_steps(nsteps, partit) ! Coded by Qiang Wang diff --git a/src/gen_modules_config.F90 b/src/gen_modules_config.F90 index a12ff8c9a..d557ed771 100755 --- a/src/gen_modules_config.F90 +++ b/src/gen_modules_config.F90 @@ -151,7 +151,13 @@ module g_config use_cavity_partial_cell, cavity_partial_cell_thresh, & use_cavity_fw2press, toy_ocean, which_toy, flag_debug, flag_warn_cflz, lwiso, & use_transit - + +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 number of groups for multi FESOM group loop parallelization + integer :: num_fesom_groups=1 + namelist /run_config_tp/ num_fesom_groups +#endif + !_____________________________________________________________________________ ! *** others *** real(kind=WP) :: dt diff --git a/src/gen_surface_forcing.F90 b/src/gen_surface_forcing.F90 index bcc4de87b..182ec3b77 100644 --- a/src/gen_surface_forcing.F90 +++ b/src/gen_surface_forcing.F90 @@ -1220,9 +1220,22 @@ SUBROUTINE sbc_ini(partit, mesh) ! OPEN and read namelist for SBC REcoM open( unit=nm_sbc_unit+1, file='namelist.recom', form='formatted', access='sequential', status='old', iostat=iost ) if (iost == 0) then +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) WRITE(*,*) ' file : ', 'namelist.recom for sbc',' open ok' +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif else +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) WRITE(*,*) 'ERROR: --> bad opening file : ', 'namelist.recom for sbc',' ; iostat=',iost +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif call par_ex(partit%MPI_COMM_FESOM, partit%mype) stop endif @@ -1246,6 +1259,7 @@ SUBROUTINE sbc_do(partit, mesh) #if defined (__recom) use recom_config use recom_glovar + use REcoM_ciso #endif IMPLICIT NONE @@ -1263,6 +1277,11 @@ SUBROUTINE sbc_do(partit, mesh) real(kind=8), allocatable :: ncdata(:) integer :: CO2start, CO2count integer :: status, ncid, varid + character(300) :: sedfilename + logical :: do_read=.false. + integer :: n_lb + integer, dimension(2) :: istart, icount + real(kind=8) :: total_runoff #endif type(t_partit), intent(inout), target :: partit type(t_mesh), intent(in), target :: mesh @@ -1414,20 +1433,61 @@ SUBROUTINE sbc_do(partit, mesh) end if -#if defined (__recom) +#if defined(__recom) !< read surface atmospheric deposition for Fe, N, CO2 if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> Atm_input'//achar(27)//'[0m' ! ******** Atmospheric CO2 ********* if (mstep == 1) then ! The year has changed + if (use_atbox) then +! Atmospheric box model CO2 values + AtmCO2(:) = x_co2atm(1) + if (ciso) then + AtmCO2_13(:) = x_co2atm_13(1) + if (ciso_14) AtmCO2_14(:,1) = x_co2atm_14(1) + end if + else +! Prescribed atmospheric CO2 values if (constant_CO2) then AtmCO2(:) = CO2_for_spinup +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Constant_CO2 = ', CO2_for_spinup - if (mype==0) write(*,*),'Atm CO2=', AtmCO2 + if (mype==0) write(*,*),'Atm CO2=', AtmCO2 +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + if (ciso) then + AtmCO2_13 = CO2_for_spinup * (1. + 0.001 * delta_co2_13) + if (ciso_14) then +! Atmospheric 14C varies with latitude + do i=1, myDim_nod2D +! Latitude of atmospheric input data + lat_val = geo_coord_nod2D(2,i) / rad +! Binning to latitude zones + if (ciso_organic_14) then +! Convert Delta_14C to delta_14C + delta_co2_14 = (big_delta_co2_14(lat_zone(lat_val)) + 2. * delta_co2_13 + 50.) / & + (0.95 - 0.002 * delta_co2_13) + else +! "Inorganic" 14C approximation: delta_14C := Delta_14C + delta_co2_14 = big_delta_co2_14(lat_zone(lat_val)) + end if + AtmCO2_14(lat_zone(lat_val),:) = CO2_for_spinup * (1. + 0.001 * delta_co2_14) + end do + end if + end if else filename=trim(nm_co2_data_file) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Updating CO2 climatology for month ', i,' from ', trim(filename) +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif totnumyear = lastyearoffesomcycle-firstyearoffesomcycle+1 firstyearofcurrentCO2cycle = lastyearoffesomcycle-numofCO2cycles*totnumyear+(currentCO2cycle-1)*totnumyear @@ -1454,12 +1514,317 @@ SUBROUTINE sbc_do(partit, mesh) status=nf_get_vara_double(ncid,varid,CO2start,CO2count,ncdata) AtmCO2(:)=ncdata(:) deallocate(ncdata) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*),'Current carbon year=',currentCO2year if (mype==0) write(*,*),'Atm CO2=', AtmCO2 +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif status=nf_close(ncid) end if + end if ! atmospheric box model or prescribed CO2 values + +! Control output of atmospheric CO2 values + if (mype==0) then !OG +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + print *, "In atm_input: AtmCO2 = ", AtmCO2(1) + if (ciso) then + print *, " AtmCO2_13 = ", AtmCO2_13(1) + if (ciso_14) print *, " AtmCO2_14 = ", AtmCO2_14(:,1) + end if + if (use_atbox) print *, " use_atbox = .true." +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + end if + + + +! ******** Sediment input ********* +!-Checking if files need to be opened--------------------------------------------- + if(use_MEDUSA .and. (sedflx_num .ne. 0)) then + allocate(ncdata(9)) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> Sed_input'//achar(27)//'[0m' +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + ! MEDUSA input needs to be renamed via jobscript + sedfilename = trim(ResultPath)//'medusa_flux2fesom.nc' +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) write(*,*) 'Updating sedimentary input first time from', sedfilename +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + +!-Opening files-------------------------------------------------------------------- + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_din', 1, GloSed(:,1), partit,mesh) +! if (mype==0) write(*,*) mype, 'sediment DIN flux:', maxval(GloSed(:,1)), minval(GloSed(:,1)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic', 1, GloSed(:,2), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC flux:', maxval(GloSed(:,2)), minval(GloSed(:,2)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_alk', 1, GloSed(:,3), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment Alk flux:', maxval(GloSed(:,3)), minval(GloSed(:,3)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dsi', 1, GloSed(:,4), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DSi flux:', maxval(GloSed(:,4)), minval(GloSed(:,4)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_o2', 1, GloSed(:,5), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment O2 flux:', maxval(GloSed(:,5)), minval(GloSed(:,5)) + + if(ciso) then + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic13', 1, GloSed(:,6), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC13 flux:', maxval(GloSed(:,6)), minval(GloSed(:,6)) + if(ciso_14) then + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic14', 1, GloSed(:,7), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC14 flux:', maxval(GloSed(:,7)), minval(GloSed(:,7)) + end if ! ciso_14 + end if ! ciso + +! unit conversion + GloSed(:,:)=GloSed(:,:)/86400 + +! read loopback fluxes from the same file + if(add_loopback) then +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) write(*,*) 'adding loopback fluxes through runoff for the first time' !OG +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + + istart = (/1,1/) + icount = (/1,1/) + ncdata = 0.d0 + + total_runoff = 8.76d5*86400 + + status=nf_open(sedfilename, nf_nowrite, ncid) + if(status.ne.nf_noerr) call handle_err(status) + + status=nf_inq_varid(ncid, 'loopback_orgm_din', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(1)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_din (mmolN/day):', ncdata(1) !OG + + status=nf_inq_varid(ncid, 'loopback_orgm_dic', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(2)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_dic (mmolC/day):', ncdata(2) !OG + + status=nf_inq_varid(ncid, 'loopback_orgm_alk', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(3)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_alk (mmolAlk/day):', ncdata(3) !OG + + status=nf_inq_varid(ncid, 'loopback_opal', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(4)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_opal (mmolSi/day):', ncdata(4) !OG + + status=nf_inq_varid(ncid, 'loopback_caco3', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(5)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco3 (mmolC/day):', ncdata(5) !OG + + if(ciso) then + status=nf_inq_varid(ncid, 'loopback_orgm_dic13', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(6)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_dic13:', ncdata(6) !OG + + status=nf_inq_varid(ncid, 'loopback_caco313', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(7)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco313:', ncdata(7)!OG + + if(ciso_14 .and. ciso_organic_14) then + status=nf_inq_varid(ncid, 'loopback_orgm_dic14', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(8)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_dic14:', ncdata(8) !OG + + status=nf_inq_varid(ncid, 'loopback_caco314', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(9)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco314:', ncdata(9) !OG + + end if ! ciso_14 .and. ciso_organic_14 + end if ! ciso + deallocate(ncdata) + status=nf_close(ncid) + +! calculating fluxes back to ocean surface through rivers (mmol/m2/s) +! converting from fluxes out of sediment to fluxes into the ocean + do n_lb = 1,9 + lb_flux(:,n_lb) = -runoff*ncdata(n_lb)/total_runoff*lb_tscale + end do + + end if ! add_loopback + + else + +!-Checking if files need to be opened--------------------------------------------- + call monthly_event(do_read) + if(do_read) then ! file is opened and read every year + i=month + if (i > 12) i=1 +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) write(*,*) 'Updating sedimentary input for month', i, 'from', sedfilename !OG +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_din', 1, GloSed(:,1), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIN flux:', maxval(GloSed(:,1)), minval(GloSed(:,1)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic', 1, GloSed(:,2), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC flux:', maxval(GloSed(:,2)), minval(GloSed(:,2)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_alk', 1, GloSed(:,3), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment Alk flux:', maxval(GloSed(:,3)), minval(GloSed(:,3)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dsi', 1, GloSed(:,4), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DSi flux:', maxval(GloSed(:,4)), minval(GloSed(:,4)) + + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_o2', 1, GloSed(:,5), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment O2 flux:', maxval(GloSed(:,5)), minval(GloSed(:,5)) + + if(ciso) then + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic13', 1, GloSed(:,6), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC13 flux:', maxval(GloSed(:,6)), minval(GloSed(:,6)) + if(ciso_14) then + call read_2ddata_on_grid_NetCDF(sedfilename, 'df_dic14', 1, GloSed(:,7), partit, mesh) +! if (mype==0) write(*,*) mype, 'sediment DIC14 flux:', maxval(GloSed(:,7)), minval(GloSed(:,7)) + end if ! ciso_14 + end if ! ciso + +!to mmol/m2/s + GloSed(:,:)=GloSed(:,:)/86400 + +! read loopback fluxes from the same file + if(add_loopback) then +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) write(*,*) 'adding loopback fluxes into the ocean monthly' !OG +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + istart = (/1,1/) + icount = (/1,1/) + ncdata = 0.d0 + + total_runoff = 8.76d5*86400 + + status=nf_open(sedfilename, nf_nowrite, ncid) + if(status.ne.nf_noerr) call handle_err(status) + + status=nf_inq_varid(ncid, 'loopback_orgm_din', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(1)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_din (mmolN/day):', ncdata(1) !OG + + status=nf_inq_varid(ncid, 'loopback_orgm_dic', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(2)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_dic (mmolC/day):', ncdata(2) !OG + + status=nf_inq_varid(ncid, 'loopback_orgm_alk', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(3)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_orgm_alk (mmolAlk/day):', ncdata(3) !OG + + status=nf_inq_varid(ncid, 'loopback_opal', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(4)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_opal (mmolSi/day):', ncdata(4) !OG + + status=nf_inq_varid(ncid, 'loopback_caco3', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(5)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco3 (mmolC/day):', ncdata(5) !OG + + if(ciso) then + status=nf_inq_varid(ncid, 'loopback_orgm_dic13', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(6)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_dic13:', ncdata(6) !OG + + status=nf_inq_varid(ncid, 'loopback_caco313', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(7)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco313:', ncdata(7) !OG + + if(ciso_14 .and. ciso_organic_14) then + status=nf_inq_varid(ncid, 'loopback_orgm_dic14', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(8)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_dic14:', ncdata(8) !OG + + status=nf_inq_varid(ncid, 'loopback_caco314', varid) + if(status.ne.nf_noerr) call handle_err(status) + status=nf_get_vara_double(ncid,varid,istart,icount,ncdata(9)) + if(status.ne.nf_noerr) call handle_err(status) +! if (mype==0) write(*,*) mype, 'loopback_caco314:', ncdata(9) !OG + + end if ! ciso_14 .and. ciso_organic_14 + end if ! ciso + status=nf_close(ncid) + +! calculating fluxes back to ocean surface through rivers (mmol/m2/s) +! converting from fluxes out of sediment to fluxes into the ocean + do n_lb = 1,9 + lb_flux(:,n_lb) = -runoff*ncdata(n_lb)/total_runoff*lb_tscale + end do + + end if ! add_loopback + + end if ! do_read + +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) write(*,*) 'sedimentary input from MEDUSA not used!' !OG +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + end if ! use_MEDUSA and sedflx_num not 0 + end if + ! ******** Fe deposition ********* if (fe_data_source=='Albani') then if (update_monthly_flag) then @@ -1467,11 +1832,23 @@ SUBROUTINE sbc_do(partit, mesh) if (mstep > 1) i=i+1 if (i > 12) i=1 filename=trim(nm_fe_data_file) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Updating iron climatology for month ', i,' from ', trim(filename) +#if defined(__usetp) + endif +#endif call read_2ddata_on_grid_NetCDF(filename,'DustClim', i, GloFeDust, partit, mesh) end if else +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Albani is switched off --> Check namelist.recom' +#if defined(__usetp) + endif +#endif end if ! ******** N deposition ********* @@ -1483,7 +1860,13 @@ SUBROUTINE sbc_do(partit, mesh) ! if (i > 12) i=1 ! if (mype==0) write(*,*) 'Updating iron climatology for month ', i filename=trim(nm_aen_data_file) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Updating nitrogen climatology for month ', i,' from ', trim(filename) +#if defined(__usetp) + endif +#endif if (yearnew .gt. 2009) then Nvari = 'NDep2009' else if (yearnew .lt. 1850) then @@ -1496,7 +1879,13 @@ SUBROUTINE sbc_do(partit, mesh) end if else GloNDust = 0.0_WP +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mstep==1 .and. mype==0) write(*,*) 'useAeolianN is switched off' +#if defined(__usetp) + endif +#endif end if ! ******** Riverine input (Nutrients) ********* @@ -1513,7 +1902,13 @@ SUBROUTINE sbc_do(partit, mesh) if (mstep > 1) i=i+1 if (i > 12) i=1 filename=trim(nm_river_data_file) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Updating riverine restoring data for month', i,' from ', trim(filename) +#if defined(__usetp) + endif +#endif call read_2ddata_on_grid_NetCDF(filename,'Alkalinity', i, RiverAlk2D, partit, mesh) ! write(*,*) mype, 'RiverAlk2D', maxval(RiverAlk2D(:)), minval(RiverAlk2D(:)) ! molar convertion of [CaCo3] * 2 -> [total Alkalinity] @@ -1535,7 +1930,13 @@ SUBROUTINE sbc_do(partit, mesh) end if else is_riverinput = 0.0d0 +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0 .and. mstep==1) write(*,*) 'No riverine input' +#if defined(__usetp) + endif +#endif end if ! ******** Riverine input of iron ********* @@ -1553,7 +1954,13 @@ SUBROUTINE sbc_do(partit, mesh) !< read erosion input ! *** River inputs are in mmol/m2/s *** ! add erosion nutrients as surface boundary condition (surface_bc function in oce_ale_tracers) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> Erosion_input'//achar(27)//'[0m' +#if defined(__usetp) + endif +#endif is_erosioninput = 1.0d0 @@ -1562,7 +1969,13 @@ SUBROUTINE sbc_do(partit, mesh) if (mstep > 1) i=i+1 if (i > 12) i=1 filename=trim(nm_erosion_data_file) +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) 'Updating erosion restoring data for month ', i,' from ', trim(filename) +#if defined(__usetp) + endif +#endif call read_2ddata_on_grid_NetCDF(filename,'POC', i, ErosionTOC2D, partit, mesh) ! write(*,*) mype, 'ErosionTOC2D', maxval(ErosionTOC2D(:)), minval(ErosionTOC2D(:)) @@ -1574,7 +1987,13 @@ SUBROUTINE sbc_do(partit, mesh) end if else is_erosioninput = 0.0d0 +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0 .and. mstep==1) write(*,*) 'No erosion input' +#if defined(__usetp) + endif +#endif end if #endif diff --git a/src/int_recom/recom_atbox.F90 b/src/int_recom/recom_atbox.F90 new file mode 100644 index 000000000..4e579c51e --- /dev/null +++ b/src/int_recom/recom_atbox.F90 @@ -0,0 +1,92 @@ + subroutine recom_atbox(partit, mesh) +! Simple 0-d box model to calculate the temporal evolution of atmospheric CO2. +! Initially the box model was part of module recom_ciso. Now it can be run also +! without carbon isotopes (ciso==.false.) +! mbutzin, 2021-07-08 + +! Settings are copied from subroutine bio_fluxes, +! some of the following modules may be unnecessary here +! use REcoM_declarations +! use REcoM_LocVar + use REcoM_GloVar + use recom_config + use recom_ciso + + use mod_mesh + USE MOD_PARTIT + USE MOD_PARSUP + + use g_config + use o_arrays + use g_comm_auto + use g_forcing_arrays + use g_support + + + implicit none + integer :: n, elem, elnodes(3),n1 + real(kind=WP) :: total_co2flux, & ! (mol / s) + total_co2flux_13, & ! (mol / s) carbon-13 + total_co2flux_14 ! (mol / s) radiocarbon + real(kind=WP), parameter :: mol_allatm = 1.7726e20 ! atmospheric inventory of all compounds (mol) + type(t_partit), intent(inout), target :: partit + type(t_mesh) , intent(inout), target :: mesh + +#include "../associate_part_def.h" +#include "../associate_mesh_def.h" +#include "../associate_part_ass.h" +#include "../associate_mesh_ass.h" + +! Globally integrated air-sea CO2 flux (mol / s) + total_co2flux = 0. + call integrate_nod(0.001 * GloCO2flux_seaicemask , total_CO2flux, partit, mesh) + +! Atmospheric carbon budget (mol) +! mass of the dry atmosphere = 5.1352e18 kg (Trenberth & Smith 2005, doi:10.1175/JCLI-3299.1) +! mean density of air = 0.02897 kg / mol (https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html) +! => total molecular inventory of the dry atmosphere: moles_atm = 1.7726e20 mol == constant. +! mol_co2atm = mol_co2atm - total_co2flux * dt +! Atmospheric mixing ratios in ppm +! x_co2atm(1) = mol_co2atm / mol_allatm * 1.e6 ! ppm + x_co2atm(1) = x_co2atm(1) - total_co2flux / mol_allatm * dt * 1.e6 + x_co2atm = x_co2atm(1) + + if (ciso) then +! Consider 13CO2 (and maybe also 14CO2) + +! Globally integrated air-sea 13CO2 flux (mol / s) + total_co2flux_13 = 0. + call integrate_nod(0.001 * GloCO2flux_seaicemask_13, total_co2flux_13, partit, mesh) + +! Atmospheric carbon-13 budget (mol) +! mol_co2atm_13 = mol_co2atm_13 - total_co2flux_13 * dt +! Budget in terms of the 13C / 12C volume mixing ratio +! x_co2atm_13(1) = mol_co2atm_13 / mol_allatm * 1.e6 + x_co2atm_13(1) = x_co2atm_13(1) - total_co2flux_13 / mol_allatm * dt * 1.e6 + x_co2atm_13 = x_co2atm_13(1) + + if (ciso_14) then + total_co2flux_14 = 0. ! globally integrated air-sea 14CO2 flux (mol / s) + call integrate_nod(0.001 * GloCO2flux_seaicemask_14, total_co2flux_14, partit, mesh) +! Atmospheric radiocarbon budget in mol: +! mol_co2atm_14 = mol_co2atm_14 + dt * (cosmic_14(1) - mol_co2atm_14 * lambda_14 - total_co2flux_14) +! = (mol_co2atm_14 + dt * (cosmic_14(1) - total_co2flux_14)) / (1 + lambda_14 * dt) +! Budget in terms of the 14C / 12C volume mixing ratio + x_co2atm_14(1) = (x_co2atm_14(1) + dt * (cosmic_14(1) - total_co2flux_14) / mol_allatm * 1.e6) / & + (1 + lambda_14 * dt) + x_co2atm_14 = x_co2atm_14(1) + +! Adjust cosmogenic 14C production (mol / s) in spinup runs, + r_atm_14 = x_co2atm_14(1) / x_co2atm(1) +! r_atm_spinup_14 is calculated once-only in subroutine recom_init + if (atbox_spinup .and. abs(r_atm_14 - r_atm_spinup_14) > 0.001) then + cosmic_14(1) = cosmic_14(1) * (r_atm_spinup_14 / r_atm_14) +! cosmic_14(1) = cosmic_14(1) * (1 + 0.01 * (r_atm_14_spinup / r_atm_14)) + end if + cosmic_14 = cosmic_14(1) + endif + end if + + return + end subroutine recom_atbox + diff --git a/src/int_recom/recom_extra.F90 b/src/int_recom/recom_extra.F90 index 5495e0c8f..eb9011cfd 100644 --- a/src/int_recom/recom_extra.F90 +++ b/src/int_recom/recom_extra.F90 @@ -201,4 +201,3 @@ subroutine krill_resp(n, partit, mesh) end if endif end subroutine krill_resp - diff --git a/src/int_recom/recom_forcing.F90 b/src/int_recom/recom_forcing.F90 index 81650779d..161fa1625 100644 --- a/src/int_recom/recom_forcing.F90 +++ b/src/int_recom/recom_forcing.F90 @@ -272,6 +272,54 @@ subroutine REcoM_Forcing(zNodes, n, Nn, state, SurfSW, Loc_slp , Temp, Sali, Sal if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> ciso after REcoM_Forcing'//achar(27)//'[0m' + if (ciso) then +! Calculate carbon-isotopic fractionation, radioactive decay is calculated in oce_ale_tracer.F90 + +! Fractionation due to air-sea exchange and chemical speciation of CO2 + call recom_ciso_airsea(recom_t(1), co3(1), recom_dic(1)) ! -> alpha_aq, alpha_dic. CO3 is taken from mocsy + +! Isotopic ratios of dissolved CO2, also needed to calculate biogenic fractionation + r_dic_13 = max(tiny*1e-3,state(1,idic_13)*1e-3) / recom_dic(1) + r_co2s_13 = alpha_aq_13 / alpha_dic_13 * r_dic_13 +! Calculate air-sea fluxes of 13|14CO2 in mmol / m**2 / s + kwco2 = kw660(1) * (660/scco2(REcoM_T(1)))**0.5 ! Piston velocity (via mocsy) + co2sat = co2flux(1) / (kwco2 + tiny) + co2(1) ! Saturation concentration of CO2 (via mocsy) +! co2flux_13 = kwco2 * alpha_k_13 * (alpha_aq_13 * r_atm_13 * co2sat - r_co2s_13 * co2(1)) +! co2flux_13 = alpha_k_13 * alpha_aq_13 * kwco2 * (r_atm_13 * co2sat - r_dic_13 * co2(1) / alpha_dic_13) +! Fractionation factors were determined for freshwater, include a correction for enhanced fractionation in seawater + co2flux_13 = (alpha_k_13 * alpha_aq_13 - 0.0002) * kwco2 * (r_atm_13 * co2sat - r_dic_13 * co2(1) / alpha_dic_13) + co2flux_seaicemask_13 = co2flux_13 * 1.e3 + +! Biogenic fractionation due to photosynthesis of plankton +! phyc_13|14 and diac_13|14 are only used in REcoM_sms to calculate DIC_13|14, DOC_13|14 and DetC_13|14 + + call recom_ciso_photo(co2(1)) ! -> alpha_p + r_phyc_13 = r_co2s_13 / alpha_p_13 + r_diac_13 = r_co2s_13 / alpha_p_dia_13 + state(1:nn,iphyc_13) = max((tiny_C * r_phyc_13), (state(1:nn,iphyc) * r_phyc_13)) + state(1:nn,idiac_13) = max((tiny_C_d * r_diac_13), (state(1:nn,idiac) * r_diac_13)) + +! The same for radiocarbon, fractionation factors have been already derived above + if (ciso_14) then +! Air-sea exchange + r_dic_14 = max(tiny*1e-3,state(1,idic_14)*1e-3) / recom_dic(1) + r_co2s_14 = alpha_aq_14 / alpha_dic_14 * r_dic_14 +! co2flux_14 = kwco2 * alpha_k_14 * (alpha_aq_14 * r_atm_14 * co2sat - r_co2s_14 * co2(1)) +! Fractionation factors were determined for freshwater, include a correction for enhanced fractionation seawater + co2flux_14 = (alpha_k_14 * alpha_aq_14 - 0.0004) * kwco2 * (r_atm_14 * co2sat - r_dic_14 * co2(1) / alpha_dic_14) + co2flux_seaicemask_14 = co2flux_14 * 1.e3 +! Biogenic fractionation + if (ciso_organic_14) then + r_phyc_14 = r_co2s_14 / alpha_p_14 + r_diac_14 = r_co2s_14 / alpha_p_dia_14 + state(1:nn,iphyc_14) = max((tiny_C * r_phyc_14), (state(1:nn,iphyc) * r_phyc_14)) + state(1:nn,idiac_14) = max((tiny_C_d * r_diac_14), (state(1:nn,idiac) * r_diac_14)) + end if + end if +! Radiocarbon + end if +! ciso + !------------------------------------------------------------------------------- ! Diagnostics if (Diags) then @@ -291,10 +339,12 @@ subroutine REcoM_Forcing(zNodes, n, Nn, state, SurfSW, Loc_slp , Temp, Sali, Sal locNNAd = sum(vertNNAd(1:nn) * thick(1:nn)) locChldegd = sum(vertChldegd(1:nn) * thick(1:nn)) +#if defined (__coccos) locNPPc = sum(vertNPPc(1:nn) * thick(1:nn)) locGPPc = sum(vertGPPc(1:nn) * thick(1:nn)) locNNAc = sum(vertNNAc(1:nn) * thick(1:nn)) locChldegc = sum(vertChldegc(1:nn) * thick(1:nn)) +#endif end if end subroutine REcoM_Forcing diff --git a/src/int_recom/recom_init.F90 b/src/int_recom/recom_init.F90 index e9b5aebbf..0186199c5 100644 --- a/src/int_recom/recom_init.F90 +++ b/src/int_recom/recom_init.F90 @@ -246,6 +246,100 @@ subroutine recom_init(tracers, partit, mesh) Sinkingvel1(:,:) = 0.d0 Sinkingvel2(:,:) = 0.d0 + if (use_MEDUSA) then + allocate(GloSed(node_size,sedflx_num)) + allocate(SinkFlx(node_size,bottflx_num)) + allocate(SinkFlx_tr(node_size,bottflx_num,num_tracers)) ! kh 25.03.22 buffer sums per tracer index + + SinkFlx(:,:) = 0.d0 + SinkFlx_tr(:,:,:) = 0.0d0 ! kh 25.03.22 + GloSed(:,:) = 0.d0 + allocate(lb_flux(node_size,9)) + lb_flux(:,:) = 0.d0 + end if + + if (useRivFe) then + allocate(RiverFe ( node_size )) + RiverFe(:) = 0.d0 + end if + +! Atmospheric box model + if (use_atbox) then +! if (mype==0 .and. my_fesom_group == 0) print *, "Initializing the atmospheric isoCO2 box model ..." !OG + allocate(x_co2atm(node_size)) + x_co2atm = CO2_for_spinup + if (ciso) then + allocate(x_co2atm_13(node_size)) + r_atm_spinup_13 = 1. + 0.001 * delta_co2_13 + x_co2atm_13 = CO2_for_spinup * r_atm_spinup_13 + if (ciso_14) then + allocate(x_co2atm_14(node_size)) + allocate(cosmic_14(node_size)) + if (ciso_organic_14) then + delta_co2_14 = (big_delta_co2_14(1) + 2. * delta_co2_13 + 50.) / (0.95 - 0.002 * delta_co2_13) + else + delta_co2_14 = big_delta_co2_14(1) + end if + r_atm_spinup_14 = 1. + 0.001 * delta_co2_14 + x_co2atm_14 = CO2_for_spinup * r_atm_spinup_14 +! Conversion of initial cosmogenic 14C production rates (mol / s) to fluxes (atoms / s / cm**2) +! Since 14C values are scaled to 12C, we need to include the standard 14C / 12C ratio here: +! 1.176e-12 (Karlen et al., 1964) * 6.0221e23 (Avogadro constant) * 1.e-4 (cm**2 / m**2) +! = 7.0820e7 cm**2 / m**2 + production_rate_to_flux_14 = 7.0820e7 / ocean_area + cosmic_14 = cosmic_14_init / production_rate_to_flux_14 + end if + end if + end if ! use_atbox + + if (ciso) then +!! Define ciso variables assigning additional ciso tracer indices +! idic_13 = bgc_base_num + 1 +! iphyc_13 = bgc_base_num + 2 +! idetc_13 = bgc_base_num + 3 +! ihetc_13 = bgc_base_num + 4 +! idoc_13 = bgc_base_num + 5 +! idiac_13 = bgc_base_num + 6 +! iphycal_13 = bgc_base_num + 7 +! idetcal_13 = bgc_base_num + 8 +! idic_14 = bgc_base_num + 9 +! iphyc_14 = bgc_base_num + 10 +! idetc_14 = bgc_base_num + 11 +! ihetc_14 = bgc_base_num + 12 +! idoc_14 = bgc_base_num + 13 +! idiac_14 = bgc_base_num + 14 +! iphycal_14 = bgc_base_num + 15 +! idetcal_14 = bgc_base_num + 16 + + !< Allocate 13CO2 surface fields + allocate(GloPCO2surf_13 ( node_size )) + allocate(GloCO2flux_13 ( node_size )) + allocate(GloCO2flux_seaicemask_13 ( node_size )) + + GloPCO2surf_13 = 0.d0 + GloCO2flux_13 = 0.d0 + GloCO2flux_seaicemask_13 = 0.0d0 + + !< Allocate auxiliary inital delta13C_DIC field + allocate(delta_dic_13_init (nl-1, nod2D )) + + if (ciso_14) then + !< Allocate 14CO2 surface fields + allocate(GloPCO2surf_14 ( node_size )) + allocate(GloCO2flux_14 ( node_size )) + allocate(GloCO2flux_seaicemask_14 ( node_size )) + + GloPCO2surf_14 = 0.d0 + GloCO2flux_14 = 0.d0 + GloCO2flux_seaicemask_14 = 0.0d0 + + !< Allocate auxiliary inital d|Delta14C_DIC fields + allocate(delta_dic_14_init ( nl-1, nod2D )) + allocate(big_delta_dic_14_init ( nl-1, nod2D )) + end if ! ciso_14 + + end if ! ciso + DO i=num_tracers-bgc_num+1, num_tracers id=tracers%data(i)%ID @@ -332,20 +426,26 @@ subroutine recom_init(tracers, partit, mesh) #if defined (__coccos) & defined (__3Zoo2Det) CASE (1029) tracers%data(i)%values(:,:) = tiny_chl/chl2N_max ! CoccoN + CASE (1030) tracers%data(i)%values(:,:) = tiny_chl/chl2N_max/NCmax ! CoccoC + CASE (1031) tracers%data(i)%values(:,:) = tiny_chl ! CoccoChl + ! ******************* ! CASE 3phy 1zoo 1det ! ******************* #elif defined (__coccos) & !defined (__3Zoo2Det) CASE (1023) tracers%data(i)%values(:,:) = tiny_chl/chl2N_max ! CoccoN + CASE (1024) tracers%data(i)%values(:,:) = tiny_chl/chl2N_max/NCmax ! CoccoC + CASE (1025) tracers%data(i)%values(:,:) = tiny_chl ! CoccoChl + #endif ! ******************* @@ -354,16 +454,20 @@ subroutine recom_init(tracers, partit, mesh) #if defined (__coccos) & defined (__3Zoo2Det) CASE (1032) tracers%data(i)%values(:,:) = tiny ! Zoo3N + CASE (1033) tracers%data(i)%values(:,:) = tiny * Redfield ! Zoo3C + #elif !defined (__coccos) & defined (__3Zoo2Det) ! ******************* ! CASE 2phy 3zoo 2det ! ******************* CASE (1029) tracers%data(i)%values(:,:) = tiny ! Zoo3N + CASE (1030) tracers%data(i)%values(:,:) = tiny * Redfield ! Zoo3C + #endif END SELECT diff --git a/src/int_recom/recom_main.F90 b/src/int_recom/recom_main.F90 index c9c5575e5..c19041a81 100755 --- a/src/int_recom/recom_main.F90 +++ b/src/int_recom/recom_main.F90 @@ -94,10 +94,11 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) real(kind=8) :: SW, Loc_slp integer :: tr_num, num_tracers - integer :: nz, n, nzmin, nzmax + integer :: nz, n, nzmin, nzmax, nu1, nl1 integer :: idiags real(kind=8) :: Sali + logical :: do_update = .false. real(kind=8), allocatable :: Temp(:), Sali_depth(:), zr(:), PAR(:) real(kind=8), allocatable :: C(:,:) @@ -113,6 +114,7 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) real(kind=8), allocatable :: OmegaC_watercolumn(:) real(kind=8), allocatable :: kspc_watercolumn(:) real(kind=8), allocatable :: rhoSW_watercolumn(:) + real(kind=WP) :: ttf_rhs_bak (mesh%nl-1, tracers%num_tracers) ! local variable ! OG - tra_diag #include "../associate_part_def.h" #include "../associate_mesh_def.h" @@ -135,6 +137,17 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) if (restore_alkalinity) call bio_fluxes(tracers, partit, mesh) if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> bio_fluxes'//achar(27)//'[0m' + if (use_atbox) then ! MERGE +! Prognostic atmospheric isoCO2 + call recom_atbox(partit,mesh) +! optional I/O of isoCO2 and inferred cosmogenic 14C production; this may cost some CPU time + if (ciso .and. ciso_14) then + call annual_event(do_update) + if (do_update .and. mype==0) write (*, fmt = '(a50,2x,i6,4(2x,f6.2))') & + 'Year, xCO2 (ppm), cosmic 14C flux (at / cm² / s):', & + yearold, x_co2atm(1), x_co2atm_13(1), x_co2atm_14(1), cosmic_14(1) * production_rate_to_flux_14 + end if + end if ! ====================================================================================== !********************************* LOOP STARTS ***************************************** @@ -176,6 +189,31 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) !!---- Atmospheric CO2 in LocVar LocAtmCO2 = AtmCO2(month) +! Update of prognostic atmospheric CO2 values + if (use_atbox) then + LocAtmCO2 = x_co2atm(1) + if (ciso) then + LocAtmCO2_13 = x_co2atm_13(1) + if (ciso_14) LocAtmCO2_14 = x_co2atm_14(1) + end if + else +! Consider prescribed atmospheric CO2 values + if (ciso) then + LocAtmCO2_13 = AtmCO2_13(month) + if (ciso_14) then +! Latitude of nodal point n + lat_val = geo_coord_nod2D(2,n) / rad +! Zonally binned NH / SH / TZ 14CO2 input values + LocAtmCO2_14 = AtmCO2_14(lat_zone(lat_val), month) + end if + end if + end if ! use_atbox + + if (ciso) then + r_atm_13 = LocAtmCO2_13(1) / LocAtmCO2(1) + if (ciso_14) r_atm_14 = LocAtmCO2_14(1) / LocAtmCO2(1) + end if + !!---- Shortwave penetration SW = parFrac * shortwave(n) SW = SW * (1.d0 - a_ice(n)) @@ -206,6 +244,14 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) C(1:nzmax, tr_num-2) = tracers%data(tr_num)%values(1:nzmax, n) end do + ttf_rhs_bak = 0.0 ! OG - tra_diag + + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do tr_num=1, num_tracers + ttf_rhs_bak(1:nzmax,tr_num) = tracers%data(tr_num)%values(1:nzmax, n) + end do + end if + !!---- Depth of the nodes in the water column zr(1:nzmax) = Z_3d_n(1:nzmax, n) @@ -219,29 +265,36 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) if (Diags) then !! * Allocate 3D diagnostics * - allocate(vertrespmeso(nl-1), vertrespmacro(nl-1), vertrespmicro(nl-1)) + allocate(vertrespmeso(nl-1)) vertrespmeso = 0.d0 + +#if defined (__3Zoo2Det) + allocate(vertrespmacro(nl-1), vertrespmicro(nl-1)) vertrespmacro = 0.d0 vertrespmicro = 0.d0 - +#endif allocate(vertcalcdiss(nl-1), vertcalcif(nl-1)) vertcalcdiss = 0.d0 vertcalcif = 0.d0 - allocate(vertaggn(nl-1), vertaggd(nl-1), vertaggc(nl-1)) + allocate(vertaggn(nl-1), vertaggd(nl-1)) vertaggn = 0.d0 vertaggd = 0.d0 - vertaggc = 0.d0 - allocate(vertdocexn(nl-1), vertdocexd(nl-1), vertdocexc(nl-1)) + allocate(vertdocexn(nl-1), vertdocexd(nl-1)) vertdocexn = 0.d0 vertdocexd = 0.d0 - vertdocexc = 0.d0 - allocate(vertrespn(nl-1), vertrespd(nl-1), vertrespc(nl-1)) + allocate(vertrespn(nl-1), vertrespd(nl-1)) vertrespn = 0.d0 vertrespd = 0.d0 + +#if defined (__coccos) + allocate(vertaggc(nl-1), vertdocexc(nl-1), vertrespc(nl-1)) + vertaggc = 0.d0 + vertdocexc = 0.d0 vertrespc = 0.d0 +#endif !! * Allocate 2D diagnostics * allocate(vertNPPn(nl-1), vertGPPn(nl-1), vertNNAn(nl-1), vertChldegn(nl-1)) @@ -256,11 +309,13 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) vertNNAd = 0.d0 vertChldegd = 0.d0 +#if defined (__coccos) allocate(vertNPPc(nl-1), vertGPPc(nl-1), vertNNAc(nl-1), vertChldegc(nl-1)) vertNPPc = 0.d0 vertGPPc = 0.d0 vertNNAc = 0.d0 vertChldegc = 0.d0 +#endif end if if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> REcoM_Forcing'//achar(27)//'[0m' @@ -282,6 +337,14 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) tracers%data(tr_num)%values(1:nzmax, n) = C(1:nzmax, tr_num-2) end do + ! recom_sms + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do tr_num=1, num_tracers + tracers%work%tra_recom_sms(1:nzmax,n,tr_num) = tracers%data(tr_num)%values(1:nzmax, n) - ttf_rhs_bak(1:nzmax,tr_num) + !if (mype==0) print *, tra_recom_sms(:,:,tr_num) + end do + end if + !!---- Local variables that have been changed during the time-step are stored so they can be saved Benthos(n,1:benthos_num) = LocBenthos(1:benthos_num) GlodecayBenthos(n, 1:benthos_num) = decayBenthos(1:benthos_num)/SecondsPerDay ! convert from [mmol/m2/d] to [mmol/m2/s] @@ -304,35 +367,52 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) !! * Update 3D diagnostics * respmeso (1:nzmax,n) = vertrespmeso (1:nzmax) +#if defined (__3Zoo2Det) respmacro (1:nzmax,n) = vertrespmacro (1:nzmax) respmicro (1:nzmax,n) = vertrespmicro (1:nzmax) +#endif calcdiss (1:nzmax,n) = vertcalcdiss (1:nzmax) calcif (1:nzmax,n) = vertcalcif (1:nzmax) + aggn (1:nzmax,n) = vertaggn (1:nzmax) - aggd (1:nzmax,n) = vertaggd (1:nzmax) - aggc (1:nzmax,n) = vertaggc (1:nzmax) docexn (1:nzmax,n) = vertdocexn (1:nzmax) - docexd (1:nzmax,n) = vertdocexd (1:nzmax) - docexc (1:nzmax,n) = vertdocexc (1:nzmax) respn (1:nzmax,n) = vertrespn (1:nzmax) - respd (1:nzmax,n) = vertrespd (1:nzmax) - respc (1:nzmax,n) = vertrespc (1:nzmax) NPPn3D (1:nzmax,n) = vertNPPn (1:nzmax) + + aggd (1:nzmax,n) = vertaggd (1:nzmax) + docexd (1:nzmax,n) = vertdocexd (1:nzmax) + respd (1:nzmax,n) = vertrespd (1:nzmax) NPPd3D (1:nzmax,n) = vertNPPd (1:nzmax) + +#if defined (__coccos) + aggc (1:nzmax,n) = vertaggc (1:nzmax) + docexc (1:nzmax,n) = vertdocexc (1:nzmax) + respc (1:nzmax,n) = vertrespc (1:nzmax) NPPc3D (1:nzmax,n) = vertNPPc (1:nzmax) +#endif + +!YY: why printing this? if (recom_debug .and. mype==0) print *, achar(27)//'[36m'//' --> ciso after REcoM_Forcing'//achar(27)//'[0m' !! * Deallocating 2D diagnostics * deallocate(vertNPPn, vertGPPn, vertNNAn, vertChldegn) deallocate(vertNPPd, vertGPPd, vertNNAd, vertChldegd) +#if defined (__coccos) deallocate(vertNPPc, vertGPPc, vertNNAc, vertChldegc) +#endif !! * Deallocating 3D Diagnostics * - deallocate(vertrespmeso, vertrespmacro, vertrespmicro ) - deallocate(vertcalcdiss, vertcalcif ) - deallocate(vertaggn, vertaggd, vertaggc ) - deallocate(vertdocexn, vertdocexd, vertdocexc ) - deallocate(vertrespn, vertrespd, vertrespc ) + deallocate(vertrespmeso) +#if defined (__3Zoo2Det) + deallocate(vertrespmacro, vertrespmicro) +#endif + deallocate(vertcalcdiss, vertcalcif) + deallocate(vertaggn, vertdocexn, vertrespn) + deallocate(vertaggd, vertdocexd, vertrespd) +#if defined (__coccos) +! deallocate(vertgrazmeso_c) + deallocate(vertaggc, vertdocexc, vertrespc) +#endif end if @@ -345,6 +425,12 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) GloCO2flux(n) = dflux(1) ! [mmol/m2/d] GloCO2flux_seaicemask(n) = co2flux_seaicemask(1) ! [mmol/m2/s] GloO2flux_seaicemask(n) = o2flux_seaicemask(1) ! [mmol/m2/s] + if (ciso) then + GloCO2flux_seaicemask_13(n) = co2flux_seaicemask_13(1) ! [mmol/m2/s] + if (ciso_14) then + GloCO2flux_seaicemask_14(n) = co2flux_seaicemask_14(1) ! [mmol/m2/s] + end if + end if GloO2flux(n) = oflux(1) ! [mmol/m2/d] PAR3D(1:nzmax,n) = PAR(1:nzmax) @@ -375,6 +461,16 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) call exchange_nod(GloCO2flux_seaicemask, partit) call exchange_nod(GloO2flux_seaicemask, partit) + if (ciso) then + call exchange_nod(GloPCO2surf_13, partit) + call exchange_nod(GloCO2flux_13, partit) + call exchange_nod(GloCO2flux_seaicemask_13, partit) + if (ciso_14) then + call exchange_nod(GloPCO2surf_14, partit) + call exchange_nod(GloCO2flux_14, partit) + call exchange_nod(GloCO2flux_seaicemask_14, partit) + end if + end if do n=1, benthos_num call exchange_nod(Benthos(:,n), partit) end do @@ -388,10 +484,12 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) call exchange_nod(NNAd, partit) call exchange_nod(Chldegn, partit) call exchange_nod(Chldegd, partit) +#if defined (__coccos) call exchange_nod(NPPc, partit) call exchange_nod(GPPc, partit) call exchange_nod(NNAc, partit) call exchange_nod(Chldegc, partit) +#endif endif do n=1, benthos_num @@ -412,6 +510,7 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) call exchange_nod(OmegaC3D, partit) call exchange_nod(kspc3D, partit) call exchange_nod(rhoSW3D, partit) + end subroutine recom ! ====================================================================================== @@ -496,4 +595,4 @@ subroutine bio_fluxes(tracers, partit, mesh) relax_alk=relax_alk-net/ocean_area ! at ocean surface layer -end subroutine bio_fluxes \ No newline at end of file +end subroutine bio_fluxes diff --git a/src/int_recom/recom_modules.F90 b/src/int_recom/recom_modules.F90 index c4113377c..ca4918da8 100644 --- a/src/int_recom/recom_modules.F90 +++ b/src/int_recom/recom_modules.F90 @@ -110,7 +110,7 @@ module recom_config Logical :: REcoM_restart = .false. Integer :: bgc_num = 33 ! NEW increased the number from 28 to 34 (added coccos and respiration) ! NEW 3Zoo changed from 31 to 33 - integer :: bgc_base_num = 22 ! standard tracers + integer :: bgc_base_num = 22 ! tracer number for case 2phy 1zoo 1det Integer :: diags3d_num = 28 ! Number of diagnostic 3d tracers to be saved Real(kind=8) :: VDet = 20.d0 ! Sinking velocity, constant through the water column and positive downwards Real(kind=8) :: VDet_zoo2 = 200.d0 ! Sinking velocity, constant through the water column diff --git a/src/int_recom/recom_sinking.F90 b/src/int_recom/recom_sinking.F90 index 0971e5bd0..2cbbe48c4 100644 --- a/src/int_recom/recom_sinking.F90 +++ b/src/int_recom/recom_sinking.F90 @@ -41,11 +41,6 @@ subroutine ver_sinking_recom_benthos(tr_num, tracer, partit, mesh) end interface end module !=============================================================================== -! YY: sinking of second detritus adapted from Ozgur's code -! but not using recom_det_tracer_id, since -! second detritus has a different sinking speed than the first -! define recom_det2_tracer_id to make it consistent??? -!=============================================================================== subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) use MOD_MESH @@ -138,34 +133,134 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) !! * Particulate Organic Nitrogen * if( tracers%data(tr_num)%ID==1004 .or. & !iphyn tracers%data(tr_num)%ID==1007 .or. & !idetn - tracers%data(tr_num)%ID==1013 .or. & !idian - tracers%data(tr_num)%ID==1025 ) then !idetz2n + tracers%data(tr_num)%ID==1013 ) then !idian Benthos(n,1)= Benthos(n,1) + add_benthos_2d(n) ![mmol] + +#if defined(__usetp) +! kh 25.03.22 buffer sums per tracer index to avoid non bit identical results regarding global sums when running the tracer loop in parallel + Benthos_tr(n,1,tr_num)= Benthos_tr(n,1,tr_num) + add_benthos_2d(n) ![mmol] + + if (use_MEDUSA) then + SinkFlx_tr(n,1,tr_num) = SinkFlx_tr(n,1,tr_num) + add_benthos_2d(n) / area(1,n)/dt ![mmol/m2] + ! now SinkFlx hat the unit mmol/time step + ! but mmol/m2/time is needed for MEDUSA: thus /area + endif +#endif endif !! * Particulate Organic Carbon * if( tracers%data(tr_num)%ID==1005 .or. & !iphyc tracers%data(tr_num)%ID==1008 .or. & !idetc - tracers%data(tr_num)%ID==1014 .or. & !idiac - tracers%data(tr_num)%ID==1026 ) then !idetz2c + tracers%data(tr_num)%ID==1014 ) then Benthos(n,2)= Benthos(n,2) + add_benthos_2d(n) + +#if defined(__usetp) + Benthos_tr(n,2,tr_num)= Benthos_tr(n,2,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,2,tr_num) = SinkFlx_tr(n,2,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif endif !! *Particulate Organic Silicon * if( tracers%data(tr_num)%ID==1016 .or. & !idiasi - tracers%data(tr_num)%ID==1017 .or. & !idetsi - tracers%data(tr_num)%ID==1027 ) then !idetz2si + tracers%data(tr_num)%ID==1017 ) then Benthos(n,3)= Benthos(n,3) + add_benthos_2d(n) + +#if defined(__usetp) + Benthos_tr(n,3,tr_num)= Benthos_tr(n,3,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,3,tr_num) = SinkFlx_tr(n,3,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif endif !! * Cal * if( tracers%data(tr_num)%ID==1020 .or. & !iphycal - tracers%data(tr_num)%ID==1021 .or. & !idetcal - tracers%data(tr_num)%ID==1028 ) then !idetz2cal + tracers%data(tr_num)%ID==1021 ) then !idetcal Benthos(n,4)= Benthos(n,4) + add_benthos_2d(n) + +#if defined(__usetp) + Benthos_tr(n,4,tr_num)= Benthos_tr(n,4,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,4,tr_num) = SinkFlx_tr(n,4,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif + endif + + ! flux of 13C into the sediment + if (ciso) then + if( tracers%data(tr_num)%ID==1305 .or. & !iphyc_13 + tracers%data(tr_num)%ID==1308 .or. & !idetc_13 + tracers%data(tr_num)%ID==1314 ) then !idiac_13 + Benthos(n,5)= Benthos(n,5) + add_benthos_2d(n) +#if defined(__usetp) + Benthos_tr(n,5,tr_num)= Benthos_tr(n,5,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,5,tr_num) = SinkFlx_tr(n,5,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif + endif + + if( tracers%data(tr_num)%ID==1320 .or. & !iphycal_13 + tracers%data(tr_num)%ID==1321 ) then !idetcal_13 + Benthos(n,6)= Benthos(n,6) + add_benthos_2d(n) +#if defined(__usetp) + Benthos_tr(n,6,tr_num)= Benthos_tr(n,6,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,6,tr_num) = SinkFlx_tr(n,6,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif + endif + endif + + ! flux of 14C into the sediment + if (ciso .and. ciso_organic_14) then + if( tracers%data(tr_num)%ID==1405 .or. & !iphyc_14 + tracers%data(tr_num)%ID==1408 .or. & !idetc_14 + tracers%data(tr_num)%ID==1414 ) then !idiac_14 + Benthos(n,7)= Benthos(n,7) + add_benthos_2d(n) +#if defined(__usetp) + Benthos_tr(n,7,tr_num)= Benthos_tr(n,7,tr_num) + add_benthos_2d(n) + + if (use_MEDUSA) then + SinkFlx_tr(n,7,tr_num) = SinkFlx_tr(n,7,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif + endif + + if( tracers%data(tr_num)%ID==1420 .or. & !iphycal_14 + tracers%data(tr_num)%ID==1421 ) then !idetcal_14 + Benthos(n,8)= Benthos(n,8) + add_benthos_2d(n) +#if defined(__usetp) + Benthos_tr(n,8,tr_num)= Benthos_tr(n,8,tr_num) + add_benthos_2d(n) + if (use_MEDUSA) then + SinkFlx_tr(n,8,tr_num) = SinkFlx_tr(n,8,tr_num) + add_benthos_2d(n) / area(1,n)/dt + endif +#endif + endif + endif + end do + +#if defined(__usetp) + if(use_MEDUSA) then + do n=1, bottflx_num + call exchange_nod(SinkFlx_tr(:,n,tr_num), partit) + end do + end if ! use_MEDUSA +#endif + do n=1, benthos_num +#if defined(__usetp) + call exchange_nod(Benthos_tr(:,n,tr_num), partit) +#endif call exchange_nod(Benthos(:,n), partit) end do @@ -223,6 +318,43 @@ subroutine diff_ver_recom_expl(tr_num, tracers, partit, mesh) bottom_flux = 0._WP id = tracers%data(tr_num)%ID +#if defined(__recom) +if (use_MEDUSA .and. (sedflx_num .ne. 0)) then + !CV update: the calculation later has been changed by Ozgur in such + !a way that now the variable bottom_flux is in (mol/time) units, + !rather than a flux in (mol/time/area). I therefore multiply the + !Medusa fluxes by the area to get the same unit. + + SELECT CASE (id) + CASE (1001) + bottom_flux = GloSed(:,1) * area(1,:) ! DIN + CASE (1002) + bottom_flux = GloSed(:,2) * area(1,:) ! DIC + CASE (1003) + bottom_flux = GloSed(:,3) * area(1,:) ! Alk + CASE (1018) + bottom_flux = GloSed(:,4) * area(1,:) ! Si + CASE (1019) + bottom_flux = GloSed(:,1) * Fe2N_benthos * area(1,:) + CASE (1022) + bottom_flux = GloSed(:,5) * area(1,:) ! Oxy + CASE (1302) + if (ciso) then + bottom_flux = GloSed(:,6) * area(1,:) ! DIC_13 and Calc: DIC_13 + end if + CASE (1402) + if (ciso) then + bottom_flux = GloSed(:,7) * area(1,:) ! DIC_14 and Calc: DIC_14 + end if + CASE DEFAULT + if (partit%mype==0) then + write(*,*) 'check specified in boundary conditions' + write(*,*) 'the model will stop!' + end if + call par_ex(partit%MPI_COMM_FESOM, partit%mype) + stop + END SELECT +else SELECT CASE (id) CASE (1001) bottom_flux = GlodecayBenthos(:,1) !*** DIN [mmolN/m^2/s] *** @@ -236,6 +368,14 @@ subroutine diff_ver_recom_expl(tr_num, tracers, partit, mesh) bottom_flux = GlodecayBenthos(:,1) * Fe2N_benthos !*** DFe *** CASE (1022) bottom_flux = -GlodecayBenthos(:,2) * redO2C !*** O2 *** + CASE (1302) + if (ciso) then + bottom_flux = GlodecayBenthos(:,5) + GlodecayBenthos(:,6) !*** DIC_13 and Calc: DIC_13 *** + end if + CASE (1402) + if (ciso) then + bottom_flux = GlodecayBenthos(:,7) + GlodecayBenthos(:,8) !*** DIC_14 and Calc: DIC_14 *** + end if CASE DEFAULT if (partit%mype==0) then write(*,*) 'check specified in boundary conditions' @@ -244,6 +384,8 @@ subroutine diff_ver_recom_expl(tr_num, tracers, partit, mesh) call par_ex(partit%MPI_COMM_FESOM, partit%mype) stop END SELECT +endif ! (use_MEDUSA .and. (sedflux_num .gt. 0)) +#endif do n=1, myDim_nod2D @@ -438,7 +580,7 @@ subroutine ver_sinking_recom(tr_num, tracers, partit, mesh) endif #endif - end do + end do !nz=nzmin,nzmax+1 dt_sink = dt vd_flux = 0.0d0 @@ -479,7 +621,7 @@ subroutine ver_sinking_recom(tr_num, tracers, partit, mesh) tv= (0.5 * wPs * (trarr(nz,n) + psiM * Rj)+ & 0.5 * wM * (trarr(max(nzmin,nz-1),n) + psiP * Rj)) vd_flux(nz)= - tv*area(nz,n) - end do + end do !nz=nzmax, nzmin+1,-1 end if ! 3rd Order DST Sceheme with flux limiting if (.FALSE.) then ! simple upwind @@ -503,12 +645,12 @@ subroutine ver_sinking_recom(tr_num, tracers, partit, mesh) trarr(nz ,n)*(Wvel_flux(nz)+abs(Wvel_flux(nz)))) vd_flux(nz)= tv*area(nz,n) - end do + end do !nz=nzmin+1,nzmax end if ! simple upwind do nz=nzmin,nzmax vert_sink(nz,n) = vert_sink(nz,n) + (vd_flux(nz)-vd_flux(nz+1))*dt/areasvol(nz,n)/hnode_new(nz,n) !/(zbar_3d_n(nz,n)-zbar_3d_n(nz+1,n)) end do - end do + end do !n = 1,myDim_nod2D end if ! Vsink .gt. 0.1 end subroutine ver_sinking_recom diff --git a/src/int_recom/recom_sms.F90 b/src/int_recom/recom_sms.F90 index 044d59c1e..7259fd893 100644 --- a/src/int_recom/recom_sms.F90 +++ b/src/int_recom/recom_sms.F90 @@ -258,6 +258,50 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & if (Grazing_detritus) recipDet2 = DetZ2C / DetZ2N #endif + if (ciso) then +!< additional variables are declared in module REcoM_ciso + DIC_13 = max(tiny,state(k,idic_13) + sms(k,idic_13 )) + PhyC_13 = max(tiny_C,state(k,iphyc_13) + sms(k,iphyc_13 )) + DetC_13 = max(tiny,state(k,idetc_13) + sms(k,idetc_13 )) + HetC_13 = max(tiny,state(k,ihetc_13) + sms(k,ihetc_13 )) + EOC_13 = max(tiny,state(k,idoc_13) + sms(k,idoc_13 )) + DiaC_13 = max(tiny_C,state(k,idiac_13) + sms(k,idiac_13 )) + PhyCalc_13 = max(tiny,state(k,iphycal_13) + sms(k,iphycal_13)) + DetCalc_13 = max(tiny,state(k,idetcal_13) + sms(k,idetcal_13)) + + calc_diss_13 = alpha_dcal_13 * calc_diss + + quota_13 = PhyN / PhyC_13 + recipQuota_13 = real(one) / quota_13 + + quota_dia_13 = DiaN / DiaC_13 + recipQuota_dia_13 = real(one) / quota_dia_13 + + recipQZoo_13 = HetC_13 / HetN + + if (ciso_14) then + DIC_14 = max(tiny,state(k,idic_14) + sms(k,idic_14 )) + if (ciso_organic_14) then + PhyC_14 = max(tiny_C,state(k,iphyc_14) + sms(k,iphyc_14 )) + DetC_14 = max(tiny,state(k,idetc_14) + sms(k,idetc_14 )) + HetC_14 = max(tiny,state(k,ihetc_14) + sms(k,ihetc_14 )) + EOC_14 = max(tiny,state(k,idoc_14) + sms(k,idoc_14 )) + DiaC_14 = max(tiny_C,state(k,idiac_14) + sms(k,idiac_14 )) + PhyCalc_14 = max(tiny,state(k,iphycal_14) + sms(k,iphycal_14)) + DetCalc_14 = max(tiny,state(k,idetcal_14) + sms(k,idetcal_14)) + + calc_diss_14 = alpha_dcal_14 * calc_diss + + quota_14 = PhyN / PhyC_14 + recipQuota_14 = real(one) / quota_14 + + quota_dia_14 = DiaN / DiaC_14 + recipQuota_dia_14 = real(one) / quota_dia_14 + recipQZoo_14 = HetC_14 / HetN + end if ! ciso_organic_14 + end if ! ciso_14 + end if ! ciso + !------------------------------------------------------------------------------- !> Temperature dependence of rates !------------------------------------------------------------------------------- @@ -939,6 +983,18 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & HetRespFlux = max(zero, HetRespFlux) !!!!!!!! CHECK Judith Valid for het_resp_noredfield case as well ???????? Then move it below endif + if (ciso) then +!MB set HetRespFlux_plus = .true. in namelist.recom +! HetRespFlux_13 = max(zero, recip_res_het * arrFunc * (hetC_13 * recip_hetN_plus - redfield) * HetC_13) +! Numerically safer parametrization avoiding instable results which may result from different cutoff values -- CHECK + HetRespFlux_13 = HetRespFlux * HetC_13 / HetC +!! HetRespFlux_13 = HetRespFlux * (HetC_13 / HetC) **2 + if (ciso_14 .and. ciso_organic_14) then +! HetRespFlux_14 = max(zero, recip_res_het * arrFunc * (hetC_14 * recip_hetN_plus - redfield) * HetC_14) + HetRespFlux_14 = HetRespFlux * HetC_14 / HetC +!! HetRespFlux_14 = HetRespFlux * (HetC_14 / HetC) **2 + end if + end if !------------------------------------------------------------------------------- !< Zooplanton mortality (Quadratic) @@ -1059,6 +1115,19 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & calc_loss_gra3 = grazingFlux_phy3 * aux ! 3Zoo #endif #endif + + if (ciso) then + calcification_13 = calc_prod_ratio * Cphot * PhyC_13 * alpha_calc_13 + calcification_13 = calcification * alpha_calc_13 + calc_loss_agg_13 = aggregationRate * PhyCalc_13 + calc_loss_gra_13 = grazingFlux_phy * recipQuota_13/(PhyC_13 + tiny) * PhyCalc_13 + if (ciso_14 .and. ciso_organic_14) then + calcification_14 = calc_prod_ratio * Cphot * PhyC_14 * alpha_calc_14 + calc_loss_agg_14 = aggregationRate * PhyCalc_14 + calc_loss_gra_14 = grazingFlux_phy * recipQuota_14/(PhyC_14 + tiny) * PhyCalc_14 + end if + end if + !------------------------------------------------------------------------------- ! Sources minus sinks (SMS) !------------------------------------------------------------------------------- @@ -1219,25 +1288,25 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & if (Grazing_detritus) then #if defined (__3Zoo2Det) sms(k,idetn) = ( & - + grazingFlux_phy3 & - - grazingFlux_phy3 * grazEff3 & - + grazingFlux_dia3 & - - grazingFlux_dia3 * grazEff3 & + + grazingFlux_phy3 & ! --> grazing on small phytoplankton by third zooplankton + - grazingFlux_phy3 * grazEff3 & ! --> fraction of grazingFlux_phy3 into microzooplankton pool + + grazingFlux_dia3 & ! --> grazing on diatoms by third zooplankton + - grazingFlux_dia3 * grazEff3 & ! --> fraction of grazingFlux_dia3 into microzooplankton pool #if defined (__coccos) - + grazingFlux_Cocco3 & - - grazingFlux_Cocco3 * grazEff3 & + + grazingFlux_Cocco3 & ! --> grazing on coccolithophores by third zooplankton + - grazingFlux_Cocco3 * grazEff3 & ! --> fraction of grazingFlux_Cocco3 into microzooplankton pool + aggregationRate * CoccoN & #endif - - grazingFlux_Det * grazEff & - - grazingFlux_Det2 * grazEff2 & ! --> okay, grazing of second zoo on first detritus + - grazingFlux_Det * grazEff & ! --> grazing of first zoo (meso) on first detritus class + - grazingFlux_Det2 * grazEff2 & ! --> grazing of second zoo (macro) on first detritus class + aggregationRate * PhyN & + aggregationRate * DiaN & - + miczooLossFlux & - - reminN * arrFunc * O2Func * DetN & ! O2remin + + miczooLossFlux & ! --> microzooplankton, mortality + - reminN * arrFunc * O2Func * DetN & ! --> O2remin ) * dt_b + sms(k,idetn) #else sms(k,idetn) = ( & - + grazingFlux_phy & + + grazingFlux_phy & ! Technically it is mesooooooooooooooooo - grazingFlux_phy * grazEff & + grazingFlux_dia & - grazingFlux_dia * grazEff & @@ -1301,7 +1370,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & + aggregationRate * CoccoC & #endif - grazingFlux_Det * recipDet * grazEff & - - grazingFlux_Det2 * recipDet2 * grazEff2 & + - grazingFlux_Det2 * recipDet * grazEff2 & ! corrected recipDet2 -> recipDet + aggregationRate * PhyC & + aggregationRate * DiaC & + miczooLossFlux * recipQZoo3 & @@ -1319,7 +1388,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & + aggregationRate * CoccoC & #endif - grazingFlux_Det * recipDet * grazEff & - - grazingFlux_Det2 * recipDet2 * grazEff & !!!!!! CHECK + ! - grazingFlux_Det2 * recipDet2 * grazEff & !!!!!! CHECK + aggregationRate * phyC & + aggregationRate * DiaC & + hetLossFlux * recipQZoo & @@ -2001,6 +2070,180 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & #endif ) * redO2C * dt_b + sms(k,ioxy) ! + if (ciso) then +!------------------------------------------------------------------------------- +! DIC_13 + sms(k,idic_13) = ( & + - Cphot * PhyC_13 & + + phyRespRate * PhyC_13 & + - Cphot_Dia * DiaC_13 & + + phyRespRate_Dia * DiaC_13 & + + rho_C1 * arrFunc * EOC_13 & + + HetRespFlux_13 & + + calc_diss_13 * DetCalc_13 & + + calc_loss_gra_13 * calc_diss_guts & + - calcification_13 & + ) * dt_b + sms(k,idic_13) +!------------------------------------------------------------------------------- +! Phytoplankton C_13 + sms(k,iphyc_13) = ( & + + Cphot * PhyC_13 & + - lossC * limitFacN * PhyC_13 & + - phyRespRate * PhyC_13 & + - aggregationRate * PhyC_13 & + - grazingFlux_phy * recipQuota_13 & + ) * dt_b + sms(k,iphyc_13) +!------------------------------------------------------------------------------- +! Detritus C_13 + sms(k,idetc_13) = ( & + + grazingFlux_phy * recipQuota_13 & + - grazingFlux_phy * recipQuota_13 * grazEff & + + grazingFlux_Dia * recipQuota_dia_13 & + - grazingFlux_Dia * recipQuota_dia_13 * grazEff & + + aggregationRate * phyC_13 & + + aggregationRate * DiaC_13 & + + hetLossFlux * recipQZoo_13 & + - reminC * arrFunc * DetC_13 & + ) * dt_b + sms(k,idetc_13) +!------------------------------------------------------------------------------- +! Heterotrophic C_13 + sms(k,ihetc_13) = ( & + + grazingFlux_phy * recipQuota_13 * grazEff & + + grazingFlux_Dia * recipQuota_dia_13 * grazEff & + - hetLossFlux * recipQZoo_13 & + - lossC_z * HetC_13 & + - hetRespFlux_13 & + ) * dt_b + sms(k,ihetc_13) +!------------------------------------------------------------------------------- +! EOC_13 + sms(k,idoc_13) = ( & + + lossC * limitFacN * phyC_13 & + + lossC_d * limitFacN_dia * DiaC_13 & + + reminC * arrFunc * DetC_13 & + + lossC_z * HetC_13 & + - rho_c1 * arrFunc * EOC_13 & + + LocRiverDOC * r_iorg_13 & + ) * dt_b + sms(k,idoc_13) +!------------------------------------------------------------------------------- +! Diatom C_13 + sms(k,idiac_13) = ( & + + Cphot_dia * DiaC_13 & + - lossC_d * limitFacN_dia * DiaC_13 & + - phyRespRate_dia * DiaC_13 & + - aggregationRate * DiaC_13 & + - grazingFlux_dia * recipQuota_dia_13 & + ) * dt_b + sms(k,idiac_13) +!------------------------------------------------------------------------------- +! Small phytoplankton calcite_13 + sms(k,iphycal_13) = ( & + + calcification_13 & + - lossC * limitFacN * phyCalc_13 & + - phyRespRate * phyCalc_13 & + - calc_loss_agg_13 & + - calc_loss_gra_13 & + ) * dt_b + sms(k,iphycal_13) +!------------------------------------------------------------------------------- +! Detritus calcite_13 + sms(k,idetcal_13) = ( & + + lossC * limitFacN * phyCalc_13 & + + phyRespRate * phyCalc_13 & + + calc_loss_agg_13 & + + calc_loss_gra_13 & + - calc_loss_gra_13 * calc_diss_guts & + - calc_diss_13 * DetCalc_13 & + ) * dt_b + sms(k,idetcal_13) +!------------------------------------------------------------------------------- + if (ciso_14) then +!------------------------------------------------------------------------------- + if (ciso_organic_14) then +! DIC_14 + sms(k,idic_14) = ( & + - Cphot * PhyC_14 & + + phyRespRate * PhyC_14 & + - Cphot_Dia * DiaC_14 & + + phyRespRate_Dia * DiaC_14 & + + rho_C1 * arrFunc * EOC_14 & + + HetRespFlux_14 & + + calc_diss_14 * DetCalc_14 & + + calc_loss_gra_14 * calc_diss_guts & + - calcification_14 & + ) * dt_b + sms(k,idic_14) +!------------------------------------------------------------------------------- +! Phytoplankton C_14 + sms(k,iphyc_14) = ( & + + Cphot * PhyC_14 & + - lossC * limitFacN * PhyC_14 & + - phyRespRate * PhyC_14 & + - aggregationRate * PhyC_14 & + - grazingFlux_phy * recipQuota_14 & + ) * dt_b + sms(k,iphyc_14) +!------------------------------------------------------------------------------- +! Detritus C_14 + sms(k,idetc_14) = ( & + + grazingFlux_phy * recipQuota_14 & + - grazingFlux_phy * recipQuota_14 * grazEff & + + grazingFlux_Dia * recipQuota_dia_14 & + - grazingFlux_Dia * recipQuota_dia_14 * grazEff & + + aggregationRate * phyC_14 & + + aggregationRate * DiaC_14 & + + hetLossFlux * recipQZoo_14 & + - reminC * arrFunc * DetC_14 & + ) * dt_b + sms(k,idetc_14) +!------------------------------------------------------------------------------- +! Heterotrophic C_14 + sms(k,ihetc_14) = ( & + + grazingFlux_phy * recipQuota_14 * grazEff & + + grazingFlux_Dia * recipQuota_dia_14 * grazEff & + - hetLossFlux * recipQZoo_14 & + - lossC_z * HetC_14 & + - hetRespFlux_14 & + ) * dt_b + sms(k,ihetc_14) +!------------------------------------------------------------------------------- +! EOC_14 + sms(k,idoc_14) = ( & + + lossC * limitFacN * phyC_14 & + + lossC_d * limitFacN_dia * DiaC_14 & + + reminC * arrFunc * DetC_14 & + + lossC_z * HetC_14 & + - rho_c1 * arrFunc * EOC_14 & + + LocRiverDOC * r_iorg_14 & + ) * dt_b + sms(k,idoc_14) +!------------------------------------------------------------------------------- +! Diatom C_14 + sms(k,idiac_14) = ( & + + Cphot_dia * DiaC_14 & + - lossC_d * limitFacN_dia * DiaC_14 & + - phyRespRate_dia * DiaC_14 & + - aggregationRate * DiaC_14 & + - grazingFlux_dia * recipQuota_dia_14 & + ) * dt_b + sms(k,idiac_14) +!------------------------------------------------------------------------------- +! Small phytoplankton calcite_14 + sms(k,iphycal_14) = ( & + + calcification_14 & + - lossC * limitFacN * phyCalc_14 & + - phyRespRate * phyCalc_14 & + - calc_loss_agg_14 & + - calc_loss_gra_14 & + ) * dt_b + sms(k,iphycal_14) +!------------------------------------------------------------------------------- +! Detritus calcite_14 + sms(k,idetcal_14) = ( & + + lossC * limitFacN * phyCalc_14 & + + phyRespRate * phyCalc_14 & + + calc_loss_agg_14 & + + calc_loss_gra_14 & + - calc_loss_gra_14 * calc_diss_guts & + - calc_diss_14 * DetCalc_14 & + ) * dt_b + sms(k,idetcal_14) +!------------------------------------------------------------------------------- + else +! "Abiotic" DIC_14, identical to DIC except for radioactive decay (-> +! recom_forcing) + sms(k,idic_14) = sms(k,idic) + end if ! ciso_organic_14 + end if ! ciso_14 + end if ! ciso !------------------------------------------------------------------------------- ! Diagnostics: Averaged rates @@ -2153,6 +2396,13 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !------------------------------------------------------------------------------- ! Remineralization from the sediments into the bottom layer + + if (use_MEDUSA .and. (sedflx_num .ne. 0)) then + if (mype==0) then !OG + write(*,*) ' --> Sedimentary input of nutrients through MEDUSA' + endif + + else ! not use_MEDUSA or sedflx_num = 0 !*** DIN *** !< decayRateBenN: Remineralization rate for benthic N [day^-1] !< LocBenthos(1): Vertically integrated N concentration in benthos (1 layer) [mmolN/m^2] @@ -2175,6 +2425,29 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & decayBenthos(4) = calc_diss_ben * LocBenthos(4) ! NEW DISS changed calc_diss to calc_diss_ben to not make the dissolution omega dependent when using the switch OmegaC_diss LocBenthos(4) = LocBenthos(4) - decayBenthos(4) * dt_b + if (ciso) then +!*** DIC_13 *** We ignore isotopic fractionation during remineralization. + decayBenthos(5) = alpha_dcal_13 * decayRateBenC * LocBenthos(5) + LocBenthos(5) = LocBenthos(5) - decayBenthos(5) * dt_b +!*** Calc: DIC_13 *** + decayBenthos(6) = calc_diss_13 * LocBenthos(6) + LocBenthos(6) = LocBenthos(6) - decayBenthos(6) * dt_b ! / depth of benthos + if (ciso_14) then + if (ciso_organic_14) then +!*** DIC_14 *** We ignore isotopic fractionation during remineralization. + decayBenthos(7) = alpha_dcal_14 * decayRateBenC * LocBenthos(7) + LocBenthos(7) = LocBenthos(7) - decayBenthos(7) * dt_b +!*** Calc: DIC_14 *** + decayBenthos(8) = calc_diss_14 * LocBenthos(8) + LocBenthos(8) = LocBenthos(8) - decayBenthos(8) * dt_b ! / depth of benthos + else +! Do nothing here because sms(idic_14) is defined as sms(idic) further +! above + end if ! ciso_organic_14 + end if ! ciso_14 + end if ! ciso + endif ! use_MEDUSA + end do ! Main time loop ends diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index afc58f87f..44e702bd7 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -523,6 +523,32 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) if (use_REcoM) then call def_stream(nod2D, myDim_nod2D, 'benCalc','Benthos calcite','mmol', Benthos(:,4), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) end if +! ciso +CASE ('benC_13 ') + if (use_REcoM) then + if (ciso) then + call def_stream(nod2D, myDim_nod2D, 'benC_13','Benthos Carbon-13','mmol/m2', Benthos(:,5), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + end if +CASE ('benC_14 ') + if (use_REcoM) then + if (ciso) then + call def_stream(nod2D, myDim_nod2D, 'benC_14','Benthos Carbon-14','mmol/m2', Benthos(:,6), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + end if +CASE ('benCalc_13') + if (use_REcoM) then + if (ciso) then + call def_stream(nod2D, myDim_nod2D, 'benCalc_13','Benthos calcite-13','mmol/m2', Benthos(:,7), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + end if +CASE ('benCalc_14') + if (use_REcoM) then + if (ciso) then + call def_stream(nod2D, myDim_nod2D, 'benCalc_14','Benthos calcite-14','mmol/m2', Benthos(:,8), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + end if +!ciso CASE ('NPPn ') if (use_REcoM) then call def_stream(nod2D, myDim_nod2D, 'NPPn','Mean NPP nanophytoplankton','mmolC/m2/d', NPPn, io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) @@ -707,6 +733,32 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) else if (tracers%data(j)%ID==1002) then if (use_REcoM) then call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC', 'Dissolved Inorganic C', '[mmol/m3]', tracers%data(j)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + if (tracers%data(j)%ltra_diag) then ! OG - tra_diag + ! horizontal advection + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_hor_adv', 'Horizontal advection part of dissolved Inorganic C', '[mmol/m3]', tracers%work%tra_advhoriz(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + ! vertical advection + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_ver_adv', 'Vertical advection part of dissolved Inorganic C', '[mmol/m3]', tracers%work%tra_advvert(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + ! horizontal diffusion + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_tra_diff_part_hor_redi', 'Horizontal diffusion of dissolved Inorganic C (includes Redi diffusivity if Redi=.true.)', '[mmol/m3]', tracers%work%tra_diff_part_hor_redi(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + if (.not. tracers%data(j)%i_vert_diff) then + ! vertical diffusion (Explicit) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_tra_diff_part_ver_expl', 'Vertical diffusion of dissolved Inorganic C (Explicit)', '[mmol/m3]', tracers%work%tra_diff_part_ver_expl(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + + ! projection of horizontal Redi diffussivity onto vertical + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_tra_diff_part_ver_redi_expl', 'Projection of horizontal Redi diffussivity onto vertical for dissolved Inorganic C (Explicit)', '[mmol/m3]', tracers%work%tra_diff_part_ver_redi_expl(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + ! vertical diffusion (Implicit) + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_tra_diff_part_ver_impl', 'Vertical diffusion of dissolved Inorganic C (Implicit)', '[mmol/m3]', tracers%work%tra_diff_part_ver_impl(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + + ! recom_sms + call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'DIC_recom_sms', 'Recom SMS', '[mmol/m3]', tracers%work%tra_recom_sms(:,:,j), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + endif else if (tracers%data(j)%ID==1003) then diff --git a/src/io_restart.F90 b/src/io_restart.F90 index 9077e27af..03c2c087a 100644 --- a/src/io_restart.F90 +++ b/src/io_restart.F90 @@ -224,6 +224,24 @@ subroutine ini_bio_io(year, tracers, partit, mesh) call bio_files%def_node_var('BenSi', 'Benthos Silicate', 'mmol/m3', Benthos(:,3), mesh, partit); call bio_files%def_node_var('BenCalc', 'Benthos Calcite', 'mmol/m3', Benthos(:,4), mesh, partit); call bio_files%def_node_var('HPlus', 'Conc. of H-plus ions in the surface water', 'mol/kg', GloHplus, mesh, partit); + if (ciso) then + call bio_files%def_node_var('BenC_13', 'Benthos Carbon-13', 'mmol/m3', Benthos(:,5), mesh, partit); + call bio_files%def_node_var('BenCalc_13', 'Benthos Calcite-13', 'mmol/m3', Benthos(:,6), mesh, partit); + if (ciso_14 .and. ciso_organic_14) then + call bio_files%def_node_var('BenC_14', 'Benthos Carbon-14', 'mmol/m3', Benthos(:,7), mesh, partit); + call bio_files%def_node_var('BenCalc_14', 'Benthos Calcite-14', 'mmol/m3', Benthos(:,8), mesh, partit); + end if ! ciso_14 + end if ! ciso + if (use_atbox) then + call bio_files%def_node_var('xCO2', 'Atm. CO2 mixing ratio', 'mol / mol', x_co2atm(:), mesh, partit); + if (ciso) then + call bio_files%def_node_var('xCO2_13', 'Atm. 13CO2 mixing ratio', 'mol / mol', x_co2atm_13(:), mesh, partit); + if (ciso_14) then + call bio_files%def_node_var('xCO2_14', 'Atm. 14CO2 mixing ratio', 'mol / mol', x_co2atm_14(:), mesh, partit); + call bio_files%def_node_var('cosmic_14', 'Cosmic 14C production', 'mol / s', cosmic_14(:), mesh, partit); + end if + end if + end if ! use_atbox end subroutine ini_bio_io #endif @@ -252,8 +270,19 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr logical rawfiles_exist, binfiles_exist logical, save :: initialized_raw = .false. logical, save :: initialized_bin = .false. - integer mpierr - +! integer mpierr + +#if defined(__recom) && defined(__usetp) +! kh 31.03.22 + integer :: tr_arr_slice_count_fix_1 + integer :: group_i + integer :: tr_num_start + integer :: tr_num_end + integer :: tr_num_in_group + logical :: has_one_added_tracer + integer :: num_tracers +#endif + !which_readr = ... ! 0 ... read netcdf restart ! 1 ... read dump file restart (binary) @@ -261,6 +290,16 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr integer, intent(out):: which_readr integer :: cstep + +#if defined(__recom) && defined(__usetp) +! kh 31.03.22 nl is required +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + + num_tracers = tracers%num_tracers +#endif !_____________________________________________________________________________ ! initialize directory for core dump restart if(.not. initialized_raw) then @@ -268,12 +307,19 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr raw_restart_dirpath = trim(ResultPath)//"fesom_raw_restart/np"//int_to_txt(partit%npes) raw_restart_infopath = trim(ResultPath)//"fesom_raw_restart/np"//int_to_txt(partit%npes)//".info" if(raw_restart_length_unit /= "off") then - if(partit%mype == RAW_RESTART_METADATA_RANK) then + +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG master rank creates the folder +#endif + if(partit%mype == RAW_RESTART_METADATA_RANK) then ! execute_command_line with mkdir sometimes fails, use a custom implementation around mkdir from C instead - call mkdir(trim(ResultPath)//"fesom_raw_restart") ! we have no mkdir -p, create the intermediate dirs separately - call mkdir(raw_restart_dirpath) - end if - call MPI_Barrier(partit%MPI_COMM_FESOM, mpierr) ! make sure the dir has been created before we continue... + call mkdir(trim(ResultPath)//"fesom_raw_restart") ! we have no mkdir -p, create the intermediate dirs separately + call mkdir(raw_restart_dirpath) + end if +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif + call MPI_Barrier(partit%MPI_COMM_FESOM, partit%mpierr) ! make sure the dir has been created before we continue... end if end if @@ -284,12 +330,18 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr bin_restart_dirpath = trim(ResultPath)//"fesom_bin_restart/np"//int_to_txt(partit%npes) bin_restart_infopath = trim(ResultPath)//"fesom_bin_restart/np"//int_to_txt(partit%npes)//".info" if(bin_restart_length_unit /= "off") then - if(partit%mype == RAW_RESTART_METADATA_RANK) then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG +#endif + if(partit%mype == RAW_RESTART_METADATA_RANK) then ! execute_command_line with mkdir sometimes fails, use a custom implementation around mkdir from C instead - call mkdir(trim(ResultPath)//"fesom_bin_restart") ! we have no mkdir -p, create the intermediate dirs separately - call mkdir(bin_restart_dirpath) - end if - call MPI_Barrier(partit%MPI_COMM_FESOM, mpierr) ! make sure the dir has been created before we continue... + call mkdir(trim(ResultPath)//"fesom_bin_restart") ! we have no mkdir -p, create the intermediate dirs separately + call mkdir(bin_restart_dirpath) + end if +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif + call MPI_Barrier(partit%MPI_COMM_FESOM, partit%mpierr) ! make sure the dir has been created before we continue... end if end if @@ -323,13 +375,13 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr if(partit%mype == RAW_RESTART_METADATA_RANK) then inquire(file=raw_restart_infopath, exist=rawfiles_exist) end if - call MPI_Bcast(rawfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, mpierr) + call MPI_Bcast(rawfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, partit%mpierr) ! check if folder for derived type binary restarts exist if(partit%mype == RAW_RESTART_METADATA_RANK) then inquire(file=bin_restart_infopath, exist=binfiles_exist) end if - call MPI_Bcast(binfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, mpierr) + call MPI_Bcast(binfiles_exist, 1, MPI_LOGICAL, RAW_RESTART_METADATA_RANK, partit%MPI_COMM_FESOM, partit%mpierr) !___________________________________________________________________________ ! read core dump file restart @@ -359,10 +411,22 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ! read netcdf file restart else which_readr = 0 +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then !OG +#endif if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ocean'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group == 0) then ! OG +#endif call read_restart(oce_path, oce_files, partit%MPI_COMM_FESOM, partit%mype) if (use_ice) then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG +#endif if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: ice'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group == 0) then ! OG +#endif call read_restart(ice_path, ice_files, partit%MPI_COMM_FESOM, partit%mype) end if @@ -370,19 +434,34 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr !RECOM restart !read here if (REcoM_restart) then +#if defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> read restarts from netcdf file: bio'//achar(27)//'[0m' +#if defined(__usetp) + endif !(partit%my_fesom_group == 0) then +#endif call read_restart(bio_path, bio_files, partit%MPI_COMM_FESOM, partit%mype) end if #endif ! immediately create a raw core dump restart if(raw_restart_length_unit /= "off") then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG master rank reads +#endif call write_all_raw_restarts(istep, partit%MPI_COMM_FESOM, partit%mype) +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif end if ! immediately create a derived type binary restart if(bin_restart_length_unit /= "off") then ! current (total) model step --> cstep = globalstep+istep +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG +#endif call write_all_bin_restarts((/globalstep+istep, int(ctime), yearnew/), & bin_restart_dirpath, & bin_restart_infopath, & @@ -391,6 +470,9 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ice, & dynamics, & tracers ) +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif end if end if end if @@ -420,15 +502,40 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr is_bin_restart_write = is_due(trim(bin_restart_length_unit), bin_restart_length, istep) end if +#if defined(__recom) && defined(__usetp) + if(num_fesom_groups > 1) then + tr_arr_slice_count_fix_1 = 1 * (nl - 1) * (myDim_nod2D + eDim_nod2D) + + do group_i = 0, num_fesom_groups - 1 + call calc_slice(num_tracers, num_fesom_groups, group_i, tr_num_start, tr_num_end, tr_num_in_group, has_one_added_tracer) + + call MPI_Bcast(tracers%data(tr_num_start)%valuesAB(:, :), tr_arr_slice_count_fix_1 * tr_num_in_group, MPI_DOUBLE_PRECISION, group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, mpierr) + end do + end if +#endif + !_____________________________________________________________________________ ! finally write restart for netcdf, core dump and derived type binary ! write netcdf restart if(is_portable_restart_write) then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif ! if(partit%mype==0) write(*,*)'Do output (netCDF, restart) ...' if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ocean'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group == 0) then +#endif + call write_restart(oce_path, oce_files, istep) if(use_ice) then +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: ice'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group == 0) then +#endif call write_restart(ice_path, ice_files, istep) end if @@ -436,8 +543,14 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr !RECOM restart !write here if (REcoM_restart .or. use_REcoM) then - +#if defined(__usetp) + if(partit%my_fesom_group == 0) then +#endif if (partit%mype==RAW_RESTART_METADATA_RANK) print *, achar(27)//'[1;33m'//' --> write restarts to netcdf file: bio'//achar(27)//'[0m' +#if defined(__usetp) + endif !(partit%my_fesom_group == 0) then +#endif + call write_restart(bio_path, bio_files, istep) end if #endif @@ -445,12 +558,27 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr ! write core dump if(is_raw_restart_write) then + +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG master rank reads +#endif + call write_all_raw_restarts(istep, partit%MPI_COMM_FESOM, partit%mype) + +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif + end if ! write derived type binary if(is_bin_restart_write) then ! current (total) model step --> cstep = globalstep+istep + +#if defined(__recom) && defined(__usetp) + if(partit%my_fesom_group == 0) then ! OG master rank reads +#endif + call write_all_bin_restarts((/globalstep+istep, int(ctime), yearnew/), & bin_restart_dirpath, & bin_restart_infopath, & @@ -469,6 +597,10 @@ subroutine restart(istep, nstart, ntotal, l_read, which_readr, ice, dynamics, tr end if end if +#if defined(__recom) && defined(__usetp) + end if ! (my_fesom_group == 0) then +#endif + end subroutine restart ! ! diff --git a/src/io_restart_file_group.F90 b/src/io_restart_file_group.F90 index 3c7a327ca..b10aa4b0a 100644 --- a/src/io_restart_file_group.F90 +++ b/src/io_restart_file_group.F90 @@ -17,7 +17,7 @@ module restart_file_group_module type restart_file_group private - type(restart_file_type), public :: files(112) + type(restart_file_type), public :: files(200) ! .OG. 112 Before integer, public :: nfiles = 0 ! todo: allow dynamically allocated size without messing with shallow copied pointers contains diff --git a/src/mod_transit.F90 b/src/mod_transit.F90 index 6713032eb..4daffd2f4 100644 --- a/src/mod_transit.F90 +++ b/src/mod_transit.F90 @@ -16,6 +16,7 @@ MODULE mod_transit r39ar_a = 1.0, & ! 39Ar / 40 Ar (homogeneous) xarg_a = 9.34e-3, & ! Argon (homogeneous) xCO2_a = 284.32e-6, & ! CO2 (CMIP6 & OMIP-BGC: 284.32e-6 for 1700-1850, PMIP4: 190.00e-6 for 21 ka BP) + xf11_a = 0.0, & ! CFC-11 (latitude dependent) xf12_a = 0.0, & ! CFC-12 (latitude dependent) xsf6_a = 0.0 ! SF6 (latitude dependent) @@ -23,6 +24,7 @@ MODULE mod_transit real(kind=8), allocatable, dimension(:) :: r14c_nh, r14c_tz, r14c_sh, & ! 14CO2 / 12CO2, latitude-dependent (e.g., bomb 14C) r14c_ti, & ! 14CO2 / 12CO2, homogenous (e.g., IntCal) xCO2_ti, & ! CO2 + xf11_nh, xf11_sh, & ! CFC-11, latitude-dependent xf12_nh, xf12_sh, & ! CFC-12, latitude-dependent xsf6_nh, xsf6_sh ! SF6, latitude-dependent integer, allocatable, dimension(:) :: year_ce ! current year in anthropenic runs (control output) diff --git a/src/oce_adv_tra_driver.F90 b/src/oce_adv_tra_driver.F90 index 8ef6e9a4c..5dff19d7e 100644 --- a/src/oce_adv_tra_driver.F90 +++ b/src/oce_adv_tra_driver.F90 @@ -209,6 +209,27 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, #endif #endif + +! O:G - tra_diag +! LO solution +! fct_LO is zero before adv_flux_hor +! Up to now only horizontal +! contribution + + +!#if defined (__recom) + if (tracers%data(tr_num)%ltra_diag) then + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! Horizontal advection part for LO (FCT is .TRUE.) + tracers%work%tra_advhoriz(nz,n,tr_num) = fct_LO(nz,n) * dt/areasvol(nz,n)/hnode_new(nz,n) + end do + end do + end if +!#endif + ! compute the low order upwind vertical flux (explicit part only) ! zero the input/output flux before computation call adv_tra_ver_upw1(we, ttf, partit, mesh, adv_flux_ver, o_init_zero=.true.) @@ -273,6 +294,24 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, end if +! O:G - tra_diag +! LO solution +! fct_LO is zero before adv_flux_ver +! vertical contribution + +!#if defined (__recom) + if (tracers%data(tr_num)%ltra_diag) then + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! Vertical advection part for LO (FCT is .TRUE.) + tracers%work%tra_advvert (nz,n,tr_num) = (adv_flux_ver(nz, n)-adv_flux_ver(nz+1, n))*dt/areasvol(nz,n)/hnode_new(nz,n) + end do + end do + end if +!#endif + !_______________________________________________________________________ if (dynamics%use_wsplit) then !wvel/=wvel_e @@ -297,7 +336,7 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, ! do horizontal tracer advection, in case of FCT high order solution SELECT CASE(trim(tracers%data(tr_num)%tra_adv_hor)) CASE('MUSCL') - ! compute the untidiffusive horizontal flux (o_init_zero=.false.: input is the LO horizontal flux computed above) + ! compute the antidiffusive horizontal flux (o_init_zero=.false.: input is the LO horizontal flux computed above) call adv_tra_hor_muscl(vel, ttfAB, partit, mesh, opth, adv_flux_hor, edge_up_dn_grad, nboundary_lay, o_init_zero=do_zero_flux) CASE('MFCT') call adv_tra_hor_mfct(vel, ttfAB, partit, mesh, opth, adv_flux_hor, edge_up_dn_grad, o_init_zero=do_zero_flux) @@ -412,6 +451,36 @@ subroutine do_oce_adv_tra(dt, vel, w, wi, we, tr_num, dynamics, tracers, partit, end if end if !-->if ((ldiag_DVD) .and. (tr_num<=2)) then + + +! O:G - tra_diag +!#if defined (__recom) + if (tracers%data(tr_num)%ltra_diag) then + !_______________________________________________________________________ + if (trim(tracers%data(tr_num)%tra_adv_lim)=='FCT') then + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! part for LO + antidiffusive (FCT is .TRUE.) + tracers%work%tra_advhoriz(nz,n,tr_num) = tracers%work%tra_advhoriz(nz,n,tr_num) + dttf_h(nz,n)/hnode_new(nz,n) + tracers%work%tra_advvert(nz,n,tr_num) = tracers%work%tra_advvert(nz,n,tr_num) + dttf_v(nz,n)/hnode_new(nz,n) + end do + end do + else + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! (FCT .FALSE.) + tracers%work%tra_advhoriz(nz,n,tr_num) = dttf_h(nz,n)/hnode_new(nz,n) + tracers%work%tra_advvert (nz,n,tr_num) = dttf_v(nz,n)/hnode_new(nz,n) + end do + end do + end if + end if +!#endif + end subroutine do_oce_adv_tra ! ! diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 681bfa51e..30dc7d793 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -3757,6 +3757,10 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) !___________________________________________________________________________ ! write out global fields for debugging if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call write_step_info'//achar(27)//'[0m' +#if defined(__recom) && defined(__usetp) +! kh 19.11.21 + if(partit%my_fesom_group == 0) then +#endif call write_step_info(n,logfile_outfreq, ice, dynamics, tracers, partit, mesh) !___________________________________________________________________________ @@ -3769,6 +3773,10 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) ! togeather around 2.5% of model runtime if (flag_debug .and. mype==0) print *, achar(27)//'[36m'//' --> call check_blowup'//achar(27)//'[0m' call check_blowup(n, ice, dynamics, tracers, partit, mesh) +#if defined(__recom) && defined(__usetp) + endif +#endif + t10=MPI_Wtime() !___________________________________________________________________________ @@ -3781,6 +3789,11 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) rtime_oce_GMRedi = rtime_oce_GMRedi + (t6-t5) rtime_oce_solvetra = rtime_oce_solvetra + (t8-t7) rtime_tot = rtime_tot + (t10-t0)-(t10-t9) + +#if defined(__recom) && defined(__usetp) +! kh 19.11.21 + if(partit%my_fesom_group == 0) then +#endif if(mod(n,logfile_outfreq)==0 .and. mype==0) then write(*,*) '___ALE OCEAN STEP EXECUTION TIMES______________________' write(*,"(A, ES10.3)") ' Oce. Mix,Press.. :', t1-t0 @@ -3799,6 +3812,10 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) write(*,"(A, ES10.3)") ' Oce. TOTAL :', t10-t0 write(*,*) write(*,*) - end if + end if +#if defined(__recom) && defined(__usetp) + endif +#endif + end subroutine oce_timestep_ale diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 3e64502ad..82e1ed93b 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -163,6 +163,8 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) #if defined(__recom) use recom_glovar use recom_config + use recom_ciso + use o_arrays #endif use diagnostics, only: ldiag_DVD use g_forcing_param, only: use_age_tracer !---age-code @@ -173,8 +175,42 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) type(t_tracer), intent(inout), target :: tracers type(t_partit), intent(inout), target :: partit type(t_mesh) , intent(in) , target :: mesh + +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 multi FESOM group loop parallelization + integer :: num_tracers + integer :: tr_num_start_memo + +! kh 15.11.21 + integer :: group_i + integer :: tr_num_start + +! kh 19.11.21 + logical :: has_one_added_tracer + logical :: has_one_added_tracer_local_dummy + logical :: tr_num_end_local_dummy + logical :: tr_num_in_group_local_dummy + integer :: tr_num_end + logical :: tr_num_in_group_dummy + integer :: tr_arr_slice_count_fix_1 + +! kh 28.03.22 + integer :: Sinkflx_tr_slice_count_fix_1 + integer :: Benthos_tr_slice_count_fix_1 + + integer :: tr_num_start_local + integer :: tr_num_to_send + +! kh 22.11.21 + logical :: completed + + logical :: bBreak +#endif + !___________________________________________________________________________ integer :: i, tr_num, node, elem, nzmax, nzmin + real(kind=WP) :: ttf_rhs_bak (mesh%nl-1, partit%myDim_nod2D+partit%eDim_elem2D) ! local variable ! OG - tra_diag + integer :: nz, n, nu1, nl1 ! OG - tra_diag !___________________________________________________________________________ ! pointer on necessary derived types real(kind=WP), dimension(:,:,:), pointer :: UV, fer_UV @@ -194,6 +230,10 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) end if del_ttf => tracers%work%del_ttf +#if defined(__recom) && defined(__usetp) + num_tracers=tracers%num_tracers +#endif + !___________________________________________________________________________ if (SPP) call cal_rejected_salt(ice, partit, mesh) if (SPP) call app_rejected_salt(tracers%data(2)%values, partit, mesh) @@ -216,13 +256,55 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) !$OMP END PARALLEL DO end if + ! Set advective and diffusive components of total tracer fluxes to zero + ! Before tr_num loop +!#if defined (__recom) ! not necessarily should belong to recom case +! tracers%work%tra_advhoriz = 0.0 ! O:G - tra_diag +! tracers%work%tra_advvert = 0.0 + ttf_rhs_bak = 0.0 +!#endif + !___________________________________________________________________________ ! loop over all tracers !$ACC UPDATE DEVICE(dynamics%w, dynamics%w_e, dynamics%uv) !!! async(1) !!! !$ACC UPDATE DEVICE(tracers%work%fct_ttf_min, tracers%work%fct_ttf_max, tracers%work%fct_plus, tracers%work%fct_minus) !$ACC UPDATE DEVICE (mesh%helem, mesh%hnode, mesh%hnode_new, mesh%zbar_3d_n, mesh%z_3d_n) + +#if defined(__recom) && defined(__usetp) +! kh 11.11.21 multi FESOM group loop parallelization + call calc_slice(num_tracers, num_fesom_groups, partit%my_fesom_group, tr_num_start, tr_num_end, tr_num_in_group_dummy, has_one_added_tracer) + +! kh 19.11.21 + tr_arr_slice_count_fix_1 = 1 * (nl - 1) * (myDim_nod2D + eDim_nod2D) + +! kh 28.03.22 + Sinkflx_tr_slice_count_fix_1 = 1 * (myDim_nod2D + eDim_nod2D) * bottflx_num + Benthos_tr_slice_count_fix_1 = 1 * (myDim_nod2D + eDim_nod2D) * benthos_num + + tr_num_start_memo = tr_num_start + +! kh 22.11.21 + request_count = 0 +#endif + +#if defined(__recom) && defined(__usetp) + do tr_num = tr_num_start, tr_num_end +#else do tr_num=1, tracers%num_tracers - +#endif + +#if defined(__recom) + if(use_MEDUSA) then + SinkFlx = 0.0d0 +#if defined(__usetp) + SinkFlx_tr(:, :, tr_num) = 0.0d0 +#endif + endif +#if defined(__usetp) + Benthos_tr(:, :, tr_num) = 0.0d0 +#endif +#endif !__recom + ! do tracer AB (Adams-Bashfort) interpolation only for advectiv part ! needed if (flag_debug .and. mype==0) print *, achar(27)//'[37m'//' --> call init_tracers_AB'//achar(27)//'[0m' @@ -250,15 +332,52 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) end do !$OMP END PARALLEL DO +! O:G +! Save horizontal and vertical advective fluxes. +! We have the values on the nodes +! We do not know how much each edge contributes +! to the nodes it connects +! Notes from Patrick: del_ttf includes +! Low-order solution. But, del_ttf_advhoriz and +! del_ttf_advvert contain antidiffusive fluxes +! from the FCT scheme + +!if (.FALSE.) then +! O:G - tra_diag +!#if defined (__recom) +! if (tracers%data(tr_num)%ltra_diag) then +! do n=1, myDim_nod2D+eDim_nod2D +! nu1 = ulevels_nod2D(n) +! nl1 = nlevels_nod2D(n) +! do nz = nu1, nl1-1 + ! Horizontal advection part +! tracers%work%tra_advhoriz(nz,n,tr_num) = tracers%work%del_ttf_advhoriz(nz,n) + ! Vertical advection part +! tracers%work%tra_advvert (nz,n,tr_num) = tracers%work%del_ttf_advvert(nz,n) +! end do +! end do +! end if +!#endif +!endif !___________________________________________________________________________ ! diffuse tracers if (flag_debug .and. mype==0) print *, achar(27)//'[37m'//' --> call diff_tracers_ale'//achar(27)//'[0m' call diff_tracers_ale(tr_num, dynamics, tracers, ice, partit, mesh) -! Radioactive decay of 14C and 39Ar + !___________________________________________________________________________ + ! Radioactive decay of 14C and 39Ar if (tracers%data(tr_num)%ID == 14) tracers%data(tr_num)%values(:,:) = tracers%data(tr_num)%values(:,:) * exp(-decay14 * dt) if (tracers%data(tr_num)%ID == 39) tracers%data(tr_num)%values(:,:) = tracers%data(tr_num)%values(:,:) * exp(-decay39 * dt) +!YY: C14 seems to be calculated both in fesom and recom +!YY: decay differently calculated??? +#if defined(__recom) + ! radioactive decay of 14C + if (ciso_14 .and. any(c14_tracer_id == tracers%data(tr_num)%ID)) then + tracers%data(tr_num)%values(:,:) = tracers%data(tr_num)%values(:,:) * (1 - lambda_14 * dt) + end if ! ciso & ciso_14 +#endif + !___________________________________________________________________________ ! relax to salt and temp climatology if (flag_debug .and. mype==0) print *, achar(27)//'[37m'//' --> call relax_to_clim'//achar(27)//'[0m' @@ -270,10 +389,93 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) end if call exchange_nod(tracers%data(tr_num)%values(:,:), partit) !$OMP BARRIER - end do +! end do !!! !$ACC UPDATE HOST (tracers%work%fct_ttf_min, tracers%work%fct_ttf_max, tracers%work%fct_plus, tracers%work%fct_minus) & !!! !$ACC HOST (tracers%work%edge_up_dn_grad) +#if defined(__recom) && defined(__usetp) +! kh 19.11.21 broadcast tracer results to fesom groups + if(num_fesom_groups > 1) then + + do group_i = 0, num_fesom_groups - 1 + call calc_slice(num_tracers, num_fesom_groups, group_i, tr_num_start_local, tr_num_end_local_dummy, tr_num_in_group_local_dummy, has_one_added_tracer_local_dummy) + + tr_num_to_send = tr_num_start_local + (tr_num - tr_num_start_memo) + + if((tr_num == tr_num_end) .and. has_one_added_tracer) then + ! skip: if last tracer in group was added to compensate for fragementation it is skipped here and handled after the loop + else + request_count = request_count + 1 + +! kh 22.11.21 non-blocking communication overlapped with computation in loop + call MPI_IBcast(tracers%data(tr_num_to_send)%values(:, :), tr_arr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, tr_arr_requests(request_count), MPIerr) + + if(use_MEDUSA) then + call MPI_IBcast(Sinkflx_tr (:, :, tr_num_to_send), Sinkflx_tr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, SinkFlx_tr_requests(request_count), MPIerr) + endif + call MPI_IBcast(Benthos_tr (:, :, tr_num_to_send), Benthos_tr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, Benthos_tr_requests(request_count), MPIerr) + end if + end do + end if ! (num_fesom_groups > 1) then +#endif + end do ! tr_num = tr_num_start, tr_num_end + +#if defined(__recom) && defined(__usetp) +! kh 19.11.21 if tracer in group was added to compensate for fragmentation its broadcast of the last index is handled here + if(num_fesom_groups > 1) then + do group_i = 0, num_fesom_groups - 1 + call calc_slice(num_tracers, num_fesom_groups, group_i, tr_num_start, tr_num_end, tr_num_in_group_dummy, has_one_added_tracer) + + if(has_one_added_tracer) then + + request_count = request_count + 1 + + call MPI_IBcast(tracers%data(tr_num_end)%values(:, :), tr_arr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, tr_arr_requests(request_count), MPIerr) + if(use_MEDUSA) then + call MPI_IBcast(Sinkflx_tr (:, :, tr_num_end), Sinkflx_tr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, SinkFlx_tr_requests(request_count), MPIerr) + endif + call MPI_IBcast(Benthos_tr (:, :, tr_num_end), Benthos_tr_slice_count_fix_1, MPI_DOUBLE_PRECISION, & + group_i, MPI_COMM_FESOM_SAME_RANK_IN_GROUPS, Benthos_tr_requests(request_count), MPIerr) + end if + end do + end if !(num_fesom_groups > 1) then + + if(num_fesom_groups > 1) then + completed = .false. + do while (.not. completed) + call MPI_TESTALL(request_count, tr_arr_requests(:), completed, MPI_STATUSES_IGNORE, MPIerr) + end do + + if(use_MEDUSA) then + completed = .false. + do while (.not. completed) + call MPI_TESTALL(request_count, SinkFlx_tr_requests(:), completed, MPI_STATUSES_IGNORE, MPIerr) + end do + endif ! (use_MEDUSA) then + + completed = .false. + do while (.not. completed) + call MPI_TESTALL(request_count, Benthos_tr_requests(:), completed, MPI_STATUSES_IGNORE, MPIerr) + end do + end if ! (num_fesom_groups > 1) then +#endif + +#if defined(__recom) && defined(__usetp) +! kh 25.03.22 SinkFlx and Benthos values are buffered per tracer index in the loop above and now summed up to +! avoid non bit identical results regarding global sums when running the tracer loop in parallel + do tr_num = 1, num_tracers + if(use_MEDUSA) then + SinkFlx = SinkFlx + SinkFlx_tr(:, :, tr_num) + endif + Benthos = Benthos + Benthos_tr(:, :, tr_num) + end do +#endif + !___________________________________________________________________________ ! 3D restoring for "passive" tracers !!!$OMPTODO: add OpenMP later, not needed right now! @@ -366,6 +568,8 @@ subroutine diff_tracers_ale(tr_num, dynamics, tracers, ice, partit, mesh) type(t_mesh) , intent(in) , target :: mesh !___________________________________________________________________________ integer :: n, nzmax, nzmin + real(kind=WP) :: ttf_rhs_bak (mesh%nl-1, partit%myDim_nod2D+partit%eDim_nod2D) ! OG - tra_diag + integer :: nz, nu1, nl1 ! OG - tra_diag !___________________________________________________________________________ ! pointer on necessary derived types real(kind=WP), pointer :: del_ttf(:,:) @@ -379,20 +583,91 @@ subroutine diff_tracers_ale(tr_num, dynamics, tracers, ice, partit, mesh) str_bf = 0.0_WP vert_sink = 0.0_WP #endif + + ttf_rhs_bak = 0.0 ! OG - tra_diag + + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ttf_rhs_bak(nz,n) = del_ttf(nz,n) + end do + end do + end if !___________________________________________________________________________ - ! do horizontal diffusiion + ! do horizontal diffusion ! write there also horizontal diffusion rhs to del_ttf which is equal the R_T^n ! in danilovs srcipt ! includes Redi diffusivity if Redi=.true. call diff_part_hor_redi(tracers, partit, mesh) ! seems to be ~9% faster than diff_part_hor + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! horizontal diffusion (w/out Redi) + tracers%work%tra_diff_part_hor_redi(nz,n,tr_num) = (del_ttf(nz,n) - ttf_rhs_bak(nz,n)) / hnode_new(nz,n) ! Unit [Conc] + !if (mype==0) print *, tracers%work%tra_diff_part_hor_redi(nz,n,tr_num) + end do + end do + end if + + if ((.not. tracers%data(tr_num)%i_vert_diff) .and. tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ttf_rhs_bak(nz,n) = del_ttf(nz,n) + end do + end do + end if !___________________________________________________________________________ ! do vertical diffusion: explicit if (.not. tracers%data(tr_num)%i_vert_diff) call diff_ver_part_expl_ale(tr_num, tracers, partit, mesh) + + ! OG i_vert_diff = TRUE so, we dont call explicit scheme + ! If we use this, check surface forcing for recom variables (They are not updated) + if ((.not. tracers%data(tr_num)%i_vert_diff) .and. tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! vertical diffusion: explicit + tracers%work%tra_diff_part_ver_expl(nz,n,tr_num) = (del_ttf(nz,n) - ttf_rhs_bak(nz,n)) / hnode_new(nz,n) ! Unit [Conc] + !if (mype==0) print *, tra_diff_part_ver_expl(:,:,tr_num) + end do + end do + end if + ! A projection of horizontal Redi diffussivity onto vertical. This par contains horizontal ! derivatives and has to be computed explicitly! + + if (tracers%data(tr_num)%ltra_diag .and. Redi) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ttf_rhs_bak(nz,n) = del_ttf(nz,n) + end do + end do + end if + if (Redi) call diff_ver_part_redi_expl(tracers, partit, mesh) + if (tracers%data(tr_num)%ltra_diag .and. Redi) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ! Redi diffussivity onto vertical: explicit + tracers%work%tra_diff_part_ver_redi_expl(nz,n,tr_num) = (del_ttf(nz,n) - ttf_rhs_bak(nz,n)) / hnode_new(nz,n) ! Unit [Conc] + !if (mype==0) print *, tra_diff_part_ver_redi_expl(:,:,tr_num) + end do + end do + end if + ! if (recom_debug .and. mype==0) print *, tracers%data(tr_num)%ID #if defined(__recom) @@ -478,7 +753,32 @@ subroutine diff_tracers_ale(tr_num, dynamics, tracers, ice, partit, mesh) !___________________________________________________________________________ if (tracers%data(tr_num)%i_vert_diff) then ! do vertical diffusion: implicite + + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + ttf_rhs_bak(nz,n) = tracers%data(tr_num)%values(nz,n) + end do + end do + end if + + ! (w/out Redi) call diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh) + + ! vertical diffusion: implicit + if (tracers%data(tr_num)%ltra_diag) then ! OG - tra_diag + do n=1, myDim_nod2D+eDim_nod2D + nu1 = ulevels_nod2D(n) + nl1 = nlevels_nod2D(n) + do nz = nu1, nl1-1 + tracers%work%tra_diff_part_ver_impl(nz,n,tr_num) = tracers%data(tr_num)%values(nz,n) - ttf_rhs_bak(nz,n) + !if (mype==0) print *, tra_diff_part_ver_impl(:,:,tr_num) + end do + end do + end if + end if !We DO not set del_ttf to zero because it will not be used in this timestep anymore !init_tracers_AB will set it to zero for the next timestep @@ -1458,6 +1758,8 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) #if defined (__recom) use recoM_declarations use recom_glovar + use recom_config + use recom_ciso #endif use mod_transit implicit none @@ -1479,16 +1781,35 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) + relax_salt(n) - real_salt_flux(n)*is_nonlinfs) #if defined(__recom) CASE (1001) ! DIN + if (use_MEDUSA .and. add_loopback) then ! OG: add is_MEDUSA_loopback flag is_MEDUSA_loopback flag * lb_flux(n,1) + bc_surface= dt*(AtmNInput(n) + RiverDIN2D(n) * is_riverinput & + + ErosionTON2D(n) * is_erosioninput + lb_flux(n,1)) + else bc_surface= dt*(AtmNInput(n) + RiverDIN2D(n) * is_riverinput & + ErosionTON2D(n) * is_erosioninput) + end if + CASE (1002) ! DIC + if (use_MEDUSA .and. add_loopback) then + bc_surface= dt*(GloCO2flux_seaicemask(n) & + + RiverDIC2D(n) * is_riverinput & + + ErosionTOC2D(n) * is_erosioninput & + + lb_flux(n,2) + lb_flux(n,5)) + else bc_surface= dt*(GloCO2flux_seaicemask(n) & + RiverDIC2D(n) * is_riverinput & + ErosionTOC2D(n) * is_erosioninput) + end if + CASE (1003) ! Alk + if (use_MEDUSA .and. add_loopback) then + bc_surface= dt*(virtual_alk(n) + relax_alk(n) & + + RiverAlk2D(n) * is_riverinput & + + lb_flux(n,3) + lb_flux(n,5)*2) !CaCO3:Alk burial=1:2 + else bc_surface= dt*(virtual_alk(n) + relax_alk(n) & + RiverAlk2D(n) * is_riverinput) - !bc_surface=0.0_WP + end if CASE (1004:1010) bc_surface=0.0_WP CASE (1011) ! DON @@ -1498,16 +1819,53 @@ FUNCTION bc_surface(n, id, sval, nzmin, partit) CASE (1013:1017) bc_surface=0.0_WP CASE (1018) ! DSi - bc_surface=dt*(RiverDSi2D(n) * is_riverinput + ErosionTSi2D(n) * is_erosioninput) + if (use_MEDUSA .and. add_loopback) then + bc_surface=dt*(RiverDSi2D(n) * is_riverinput & + + ErosionTSi2D(n) * is_erosioninput & + + lb_flux(n,4)) + else + bc_surface=dt*(RiverDSi2D(n) * is_riverinput + ErosionTSi2D(n) * is_erosioninput) + end if + CASE (1019) ! Fe + if (useRivFe) then + bc_surface= dt*(AtmFeInput(n) + RiverFe(n)) + else bc_surface= dt*AtmFeInput(n) + end if CASE (1020:1021) ! Cal bc_surface=0.0_WP CASE (1022) ! OXY bc_surface= dt*GloO2flux_seaicemask(n) ! bc_surface=0.0_WP - CASE (1023:1035) + CASE (1023:1033) bc_surface=0.0_WP ! OG added bc for recom fields + CASE (1302) + if (ciso) then + if (use_MEDUSA .and. add_loopback) then + bc_surface= dt*(GloCO2flux_seaicemask_13(n) & + + lb_flux(n,6) + lb_flux(n,7)) + else + bc_surface= dt*(GloCO2flux_seaicemask_13(n)) + end if + else + bc_surface=0.0_WP + end if + CASE (1305:1321) + bc_surface=0.0_WP ! organic 13C + CASE (1402) + if (ciso .and. ciso_14) then + if (use_MEDUSA .and. add_loopback .and. ciso_organic_14) then + bc_surface= dt*(GloCO2flux_seaicemask_14(n) & + + lb_flux(n,8) + lb_flux(n,9)) + else + bc_surface= dt*GloCO2flux_seaicemask_14(n) + end if + else + bc_surface=0.0_WP + end if + CASE (1405:1421) + bc_surface=0.0_WP ! organic 14C #endif CASE (101) ! apply boundary conditions to tracer ID=101 bc_surface= dt*(prec_rain(n))! - real_salt_flux(n)*is_nonlinfs) @@ -1687,3 +2045,38 @@ FUNCTION transit_bc_surface(n, id, sst, sss, aice, sval, nzmin, partit, mesh) END FUNCTION +!=============================================================================== +! kh 11.11.21 divide the range specified by indexcount into fesom_group_count equal slices and calculate +! the start_index and end_index for the given fesom_group_id. +! if necessary to compensate for fragmentation, the end index of the first n slices +! might be one higher than for the remaining slices. this is indicated by end_index_is_one_higher +subroutine calc_slice(index_count, fesom_group_count, fesom_group_id, start_index, end_index, index_count_in_group, end_index_is_one_higher) +! use g_config + + implicit none + integer, intent(in) :: index_count + integer, intent(in) :: fesom_group_count + integer, intent(in) :: fesom_group_id + integer, intent(out) :: start_index + integer, intent(out) :: end_index + integer, intent(out) :: index_count_in_group + logical, intent(out) :: end_index_is_one_higher + + integer :: group_id_limit_to_adjust_end_index + + index_count_in_group = index_count / fesom_group_count + group_id_limit_to_adjust_end_index = mod(index_count, fesom_group_count) + start_index = (fesom_group_id * index_count_in_group) + 1 + +! kh 11.11.21 adjust loop start and number of loop iterations by 1 if necessary + if(fesom_group_id < group_id_limit_to_adjust_end_index) then + start_index = start_index + fesom_group_id + index_count_in_group = index_count_in_group + 1 + end_index_is_one_higher = .true. + else + start_index = start_index + group_id_limit_to_adjust_end_index + end_index_is_one_higher = .false. + end if + + end_index = start_index + index_count_in_group - 1 +end subroutine calc_slice diff --git a/src/oce_mesh.F90 b/src/oce_mesh.F90 index 1843e345b..f3eab41f1 100755 --- a/src/oce_mesh.F90 +++ b/src/oce_mesh.F90 @@ -307,7 +307,13 @@ SUBROUTINE read_mesh(partit, mesh) read(fileID,*) n ! nod2D, we know it already error_status=0 if (n/=mesh%nod2D) error_status=1 !set the error status for consistency between rpart and nod2D +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif write(*,*) 'reading '// trim(file_name) +#if defined(__recom) && defined(__usetp) + end if +#endif end if ! check the error status call MPI_BCast(error_status, 1, MPI_INTEGER, 0, MPI_COMM_FESOM, ierror) diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 57f647e25..92d767fe0 100755 --- a/src/oce_modules.F90 +++ b/src/oce_modules.F90 @@ -220,6 +220,15 @@ MODULE o_ARRAYS #if defined(__recom) real(kind=WP), allocatable :: dtr_bf(:,:), str_bf(:,:) real(kind=WP), allocatable :: vert_sink(:,:) +#if defined(__usetp) +! kh 22.11.21 +integer :: request_count +integer, allocatable :: tr_arr_requests(:), tr_arr_old_requests(:) + +! kh 28.03.22 +integer, allocatable :: SinkFlx_tr_requests(:) +integer, allocatable :: Benthos_tr_requests(:) +#endif #endif !Viscosity and diff coefs diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index 4a1e96d9d..b7cab9719 100755 --- a/src/oce_setup_step.F90 +++ b/src/oce_setup_step.F90 @@ -291,12 +291,12 @@ SUBROUTINE tracer_init(tracers, partit, mesh) !___________________________________________________________________________ ! define tracer namelist parameter integer :: num_tracers - logical :: i_vert_diff, smooth_bh_tra + logical :: i_vert_diff, smooth_bh_tra , ltra_diag ! OG - tra_diag real(kind=WP) :: gamma0_tra, gamma1_tra, gamma2_tra integer :: AB_order = 2 namelist /tracer_listsize/ num_tracers namelist /tracer_list / nml_tracer_list - namelist /tracer_general / smooth_bh_tra, gamma0_tra, gamma1_tra, gamma2_tra, i_vert_diff, AB_order + namelist /tracer_general / smooth_bh_tra, gamma0_tra, gamma1_tra, gamma2_tra, i_vert_diff, AB_order, ltra_diag ! OG - tra_diag !___________________________________________________________________________ ! pointer on necessary derived types #include "associate_part_def.h" @@ -435,6 +435,7 @@ SUBROUTINE tracer_init(tracers, partit, mesh) tracers%data(n)%valuesAB = 0. tracers%data(n)%valuesold = 0. tracers%data(n)%i_vert_diff = i_vert_diff + tracers%data(n)%ltra_diag = ltra_diag ! OG - tra_diag end do allocate(tracers%work%del_ttf(nl-1,node_size)) allocate(tracers%work%del_ttf_advhoriz(nl-1,node_size),tracers%work%del_ttf_advvert(nl-1,node_size)) @@ -447,6 +448,20 @@ SUBROUTINE tracer_init(tracers, partit, mesh) tracers%work%dvd_trflx_hor = 0.0_WP tracers%work%dvd_trflx_ver = 0.0_WP end if + if (ltra_diag) then ! OG - tra_diag + allocate(tracers%work%tra_advhoriz(nl-1,node_size,num_tracers),tracers%work%tra_advvert(nl-1,node_size,num_tracers)) + tracers%work%tra_advhoriz = 0.0_WP + tracers%work%tra_advvert = 0.0_WP + allocate(tracers%work%tra_diff_part_hor_redi(nl-1,node_size,num_tracers),tracers%work%tra_diff_part_ver_expl(nl-1,node_size,num_tracers)) + allocate(tracers%work%tra_diff_part_ver_redi_expl(nl-1,node_size,num_tracers),tracers%work%tra_diff_part_ver_impl(nl-1,node_size,num_tracers)) + allocate(tracers%work%tra_recom_sms(nl-1,node_size,num_tracers)) + tracers%work%tra_diff_part_hor_redi = 0.0_WP + tracers%work%tra_diff_part_ver_expl = 0.0_WP + tracers%work%tra_diff_part_ver_redi_expl = 0.0_WP + tracers%work%tra_diff_part_ver_impl = 0.0_WP + tracers%work%tra_recom_sms = 0.0_WP + + end if END SUBROUTINE tracer_init ! ! @@ -796,6 +811,12 @@ SUBROUTINE arrays_init(num_tracers, partit, mesh) allocate(str_bf ( nl-1, node_size )) allocate(vert_sink ( nl-1, node_size )) allocate(Alk_surf ( node_size )) +#if defined(__usetp) +! kh 22.11.21 + allocate(tr_arr_requests(num_tracers), tr_arr_old_requests(num_tracers)) + allocate(SinkFlx_tr_requests(num_tracers)) + allocate(Benthos_tr_requests(num_tracers)) +#endif #endif ! ================= ! Visc and Diff coefs @@ -905,6 +926,13 @@ SUBROUTINE arrays_init(num_tracers, partit, mesh) str_bf = 0.0_WP vert_sink = 0.0_WP Alk_surf = 0.0_WP +#if defined(__usetp) +! kh 23.03.22 + tr_arr_requests = 0 + tr_arr_old_requests = 0 + SinkFlx_tr_requests = 0 + Benthos_tr_requests = 0 +#endif #endif ! init field for pressure force @@ -1000,6 +1028,9 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) end if if (mype==0) then +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif write(*,*) print *, achar(27)//'[36m'//'*************************'//achar(27)//'[0m' print *, achar(27)//'[36m'//' --> RECOM ON'//achar(27)//'[0m' @@ -1024,6 +1055,9 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) write(*,*) 'read Nitrate climatology from:', trim(filelist(6)) write(*,*) 'read Salt climatology from:', trim(filelist(7)) write(*,*) 'read Temperature climatology from:', trim(filelist(8)) +#if defined(__usetp) + end if ! (partit%my_fesom_group==0) then +#endif end if ! read ocean state ! this must be always done! First two tracers with IDs 0 and 1 are the temperature and salinity. @@ -1040,9 +1074,10 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) ! this must be always done! First two tracers with IDs 0 and 1 are the temperature and salinity. if(mype==0) write(*,*) 'read Temperature climatology from:', trim(filelist(1)) if(mype==0) write(*,*) 'read Salinity climatology from:', trim(filelist(2)) - #endif + if(any(idlist == 14) .and. mype==0) write(*,*) 'read radiocarbon climatology from:', trim(filelist(3)) + call do_ic3d(tracers, partit, mesh) Tclim=tracers%data(1)%values @@ -1060,9 +1095,17 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) #if defined(__recom) if (restore_alkalinity) then + +#if defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) write(*,*) if (mype==0) print *, achar(27)//'[46;1m'//' --> Set surface field for alkalinity restoring'//achar(27)//'[0m' - if (mype==0) write(*,*) + if (mype==0) write(*,*),'Alkalinity restoring = true. Field is read.' +#if defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + Alk_surf = tracers%data(5)%values(1,:) ! alkalinity is the 5th tracer endif @@ -1097,26 +1140,101 @@ SUBROUTINE oce_initial_state(tracers, partit, mesh) !_______________________________________________________________________ CASE (1004:1017) tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) then write (i_string, "(I4)") i write (id_string, "(I4)") id write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif CASE (1020:1021) tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) then write (i_string, "(I4)") i write (id_string, "(I4)") id write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif CASE (1023:1033) tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif if (mype==0) then write (i_string, "(I4)") i write (id_string, "(I4)") id write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) end if - !_______________________________________________________________________ +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif +!_______________________________________________________________________ +! Carbon isotopes +! Carbon-13 + CASE (1302) + tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) then + write (i_string, "(I4)") i + write (id_string, "(I4)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + CASE (1305:1321) + tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) then + write (i_string, "(I4)") i + write (id_string, "(I4)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif +! Radiocarbon + CASE (1402) + tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) then + write (i_string, "(I4)") i + write (id_string, "(I4)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif + CASE (1405:1421) + tracers%data(i)%values(:,:)=0.0_WP +#if defined(__recom) && defined(__usetp) + if (partit%my_fesom_group==0) then +#endif + if (mype==0) then + write (i_string, "(I4)") i + write (id_string, "(I4)") id + write(*,*) 'initializing '//trim(i_string)//'th tracer with ID='//trim(id_string) + end if +#if defined(__recom) && defined(__usetp) + endif !(partit%my_fesom_group==0) then +#endif +! End of carbon isotopes section +!_______________________________________________________________________ CASE (101) ! initialize tracer ID=101 tracers%data(i)%values(:,:)=0.0_WP if (mype==0) then