diff --git a/src/int_recom/recom_init.F90 b/src/int_recom/recom_init.F90 index 8c8d8eb6e..84f0e6eaa 100644 --- a/src/int_recom/recom_init.F90 +++ b/src/int_recom/recom_init.F90 @@ -613,31 +613,12 @@ subroutine recom_init(tracers, partit, mesh) end do if (mype==0) write(*,*) "Sanity check for REcoM variables after recom_init call" - call MPI_AllREDUCE(locDINmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. DIN. =', glo - call MPI_AllREDUCE(locDINmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal min init. DIN. =', glo - - call MPI_AllREDUCE(locDICmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. DIC. =', glo - call MPI_AllREDUCE(locDICmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal min init. DIC. =', glo - call MPI_AllREDUCE(locAlkmax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. Alk. =', glo - call MPI_AllREDUCE(locAlkmin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal min init. Alk. =', glo - call MPI_AllREDUCE(locDSimax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. DSi. =', glo - call MPI_AllREDUCE(locDSimin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal min init. DSi. =', glo - call MPI_AllREDUCE(locDFemax , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. DFe. =', glo - call MPI_AllREDUCE(locDFemin , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' `-> gobal min init. DFe. =', glo - call MPI_AllREDUCE(locO2max , glo , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' |-> gobal max init. O2. =', glo - call MPI_AllREDUCE(locO2min , glo , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) - if (mype==0) write(*,*) ' `-> gobal min init. O2. =', glo + call print_global_minmax('DIN', locDINmin, locDINmax, .true.) + call print_global_minmax('DIC', locDICmin, locDICmax, .true.) + call print_global_minmax('Alk', locAlkmin, locAlkmax, .true.) + call print_global_minmax('DSi', locDSimin, locDSimax, .true.) + call print_global_minmax('DFe', locDFemin, locDFemax, .false.) + call print_global_minmax('O2', locO2min, locO2max, .false.) if (enable_3zoo2det) then is_3zoo2det=1.0_WP @@ -650,5 +631,28 @@ subroutine recom_init(tracers, partit, mesh) else is_coccos=0.0_WP endif - end subroutine recom_init +contains + +subroutine print_global_minmax(label, loc_min, loc_max, use_pipe) + character(len=*), intent(in) :: label + real(kind=8) , intent(in) :: loc_min, loc_max + logical , intent(in) :: use_pipe + + real(kind=8) :: glo_min, glo_max + character(len=3) :: prefix_min, prefix_max + + if (use_pipe) then + prefix_max = ' |' + prefix_min = ' |' + else + prefix_max = ' |' + prefix_min = ' `' + end if + + call MPI_AllREDUCE(loc_max , glo_max , 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_FESOM, MPIerr) + if (mype==0) write(*,*) prefix_max//'-> gobal max init. '//trim(label)//' =', glo_max + call MPI_AllREDUCE(loc_min , glo_min , 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_FESOM, MPIerr) + if (mype==0) write(*,*) prefix_min//'-> gobal min init. '//trim(label)//' =', glo_min +end subroutine print_global_minmax + end subroutine recom_init diff --git a/src/int_recom/recom_main.F90 b/src/int_recom/recom_main.F90 index c6805ff43..80292eeb1 100755 --- a/src/int_recom/recom_main.F90 +++ b/src/int_recom/recom_main.F90 @@ -248,11 +248,7 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) ttf_rhs_bak = 0.0 - do tr_num=1, num_tracers - if (tracers%data(tr_num)%ltra_diag) then - ttf_rhs_bak(1:nzmax,tr_num) = tracers%data(tr_num)%values(1:nzmax, n) - end if - end do + call recom_diag_backup(tracers, n, nzmax, num_tracers, ttf_rhs_bak) !!---- Depth of the nodes in the water column zr(1:nzmax) = Z_3d_n(1:nzmax, n) @@ -395,15 +391,7 @@ 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 - - do tr_num=1, num_tracers - if (tracers%data(tr_num)%ltra_diag) then - 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 if - - end do + call recom_diag_store_delta(tracers, n, nzmax, num_tracers, ttf_rhs_bak) !!---- 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) @@ -626,6 +614,34 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) call exchange_nod(kspc3D, partit) call exchange_nod(rhoSW3D, partit) +contains + +subroutine recom_diag_backup(tracers, n, nzmax, num_tracers, ttf_rhs_bak) + type(t_tracer), intent(in) :: tracers + integer , intent(in) :: n, nzmax, num_tracers + real(kind=WP) , intent(inout) :: ttf_rhs_bak(:,:) + integer :: tr_num + + do tr_num=1, num_tracers + if (tracers%data(tr_num)%ltra_diag) then + ttf_rhs_bak(1:nzmax,tr_num) = tracers%data(tr_num)%values(1:nzmax, n) + end if + end do +end subroutine recom_diag_backup + +subroutine recom_diag_store_delta(tracers, n, nzmax, num_tracers, ttf_rhs_bak) + type(t_tracer), intent(inout) :: tracers + integer , intent(in) :: n, nzmax, num_tracers + real(kind=WP) , intent(in) :: ttf_rhs_bak(:,:) + integer :: tr_num + + do tr_num=1, num_tracers + if (tracers%data(tr_num)%ltra_diag) then + 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) + end if + end do +end subroutine recom_diag_store_delta + end subroutine recom ! ====================================================================================== diff --git a/src/int_recom/recom_sinking.F90 b/src/int_recom/recom_sinking.F90 index 1d8139e0f..fb1036fdf 100644 --- a/src/int_recom/recom_sinking.F90 +++ b/src/int_recom/recom_sinking.F90 @@ -140,19 +140,7 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) tracers%data(tr_num)%ID==1007 .or. & !idetn tracers%data(tr_num)%ID==1013 .or. & !idian tracers%data(tr_num)%ID==1025 ) then !idetz2n - Benthos(n,1)= Benthos(n,1) + add_benthos_2d(n) ![mmol] - - if (use_MEDUSA) then -! 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 - 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 - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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] - endif - + call accum_benthos(n, tr_num, 1, 1, add_benthos_2d(n), area(1,n), dt, .true.) endif !! * Particulate Organic Carbon * @@ -160,51 +148,21 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) tracers%data(tr_num)%ID==1008 .or. & !idetc tracers%data(tr_num)%ID==1014 .or. & !idiac tracers%data(tr_num)%ID==1026 ) then !idetz2c - Benthos(n,2)= Benthos(n,2) + add_benthos_2d(n) - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,2,tr_num) = SinkFlx_tr(n,2,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,2,tr_num)= Benthos_tr(n,2,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 2, 2, add_benthos_2d(n), area(1,n), dt, .true.) 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 - Benthos(n,3)= Benthos(n,3) + add_benthos_2d(n) - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,3,tr_num) = SinkFlx_tr(n,3,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,3,tr_num)= Benthos_tr(n,3,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 3, 3, add_benthos_2d(n), area(1,n), dt, .true.) 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 - Benthos(n,4)= Benthos(n,4) + add_benthos_2d(n) - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,4,tr_num) = SinkFlx_tr(n,4,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,4,tr_num)= Benthos_tr(n,4,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 4, 4, add_benthos_2d(n), area(1,n), dt, .true.) endif ! flux of 13C into the sediment @@ -212,30 +170,12 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) 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_14 - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,5,tr_num) = SinkFlx_tr(n,5,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,5,tr_num)= Benthos_tr(n,5,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 0, 5, add_benthos_2d(n), area(1,n), dt, .false.) endif if( tracers%data(tr_num)%ID==1320 .or. & !iphycal tracers%data(tr_num)%ID==1321 ) then !idetcal - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,6,tr_num) = SinkFlx_tr(n,6,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,6,tr_num)= Benthos_tr(n,6,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 0, 6, add_benthos_2d(n), area(1,n), dt, .false.) endif endif @@ -245,28 +185,12 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) if( tracers%data(tr_num)%ID==1405 .or. & !iphyc_13 tracers%data(tr_num)%ID==1408 .or. & !idetc_13 tracers%data(tr_num)%ID==1414 ) then !idiac_14 - - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,7,tr_num) = SinkFlx_tr(n,7,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,7,tr_num)= Benthos_tr(n,7,tr_num) + add_benthos_2d(n) - endif - + call accum_benthos(n, tr_num, 0, 7, add_benthos_2d(n), area(1,n), dt, .false.) endif if( tracers%data(tr_num)%ID==1420 .or. & !iphycal tracers%data(tr_num)%ID==1421 ) then !idetcal - if (use_MEDUSA) then -! 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 - SinkFlx_tr(n,8,tr_num) = SinkFlx_tr(n,8,tr_num) + add_benthos_2d(n) / area(1,n)/dt - endif - if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then -! 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,8,tr_num)= Benthos_tr(n,8,tr_num) + add_benthos_2d(n) - endif + call accum_benthos(n, tr_num, 0, 8, add_benthos_2d(n), area(1,n), dt, .false.) endif endif @@ -288,6 +212,26 @@ subroutine ver_sinking_recom_benthos(tr_num, tracers, partit, mesh) call exchange_nod(Benthos(:,n), partit) end do +contains + +subroutine accum_benthos(n, tr_num, benthos_idx, flx_idx, add_val, area_surf, dt_local, do_benthos) + integer , intent(in) :: n, tr_num + integer , intent(in) :: benthos_idx, flx_idx + real(kind=WP), intent(in) :: add_val, area_surf, dt_local + logical , intent(in) :: do_benthos + + if (do_benthos) then + Benthos(n,benthos_idx)= Benthos(n,benthos_idx) + add_val + end if + + if (use_MEDUSA) then + SinkFlx_tr(n,flx_idx,tr_num) = SinkFlx_tr(n,flx_idx,tr_num) + add_val / area_surf / dt_local + endif + if ((.not.use_MEDUSA).or.(sedflx_num.eq.0)) then + Benthos_tr(n,flx_idx,tr_num)= Benthos_tr(n,flx_idx,tr_num) + add_val + endif +end subroutine accum_benthos + end subroutine ver_sinking_recom_benthos ! ! diff --git a/src/int_recom/recom_sms.F90 b/src/int_recom/recom_sms.F90 index 88c62a473..30f33ca7f 100644 --- a/src/int_recom/recom_sms.F90 +++ b/src/int_recom/recom_sms.F90 @@ -6900,10 +6900,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !--------------------------------------------------------------------------- ! Calculate N remineralization flux [mmolN m-2 day-1] - decayBenthos(1) = decayRateBenN * LocBenthos(1) - - ! Update benthic N pool (remove remineralized N) - LocBenthos(1) = LocBenthos(1) - decayBenthos(1) * dt_b + call update_benthos_decay(decayRateBenN, 1) !=========================================================================== ! DISSOLVED INORGANIC CARBON (DIC) REMINERALIZATION @@ -6930,10 +6927,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !--------------------------------------------------------------------------- ! Calculate C remineralization flux [mmolC m-2 day-1] - decayBenthos(2) = decayRateBenC * LocBenthos(2) - - ! Update benthic C pool - LocBenthos(2) = LocBenthos(2) - decayBenthos(2) * dt_b + call update_benthos_decay(decayRateBenC, 2) !=========================================================================== ! SILICATE (Si) DISSOLUTION @@ -6962,10 +6956,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !--------------------------------------------------------------------------- ! Calculate Si dissolution flux [mmolSi m-2 day-1] - decayBenthos(3) = decayRateBenSi * LocBenthos(3) - - ! Update benthic Si pool - LocBenthos(3) = LocBenthos(3) - decayBenthos(3) * dt_b + call update_benthos_decay(decayRateBenSi, 3) !=========================================================================== ! CALCITE (CaCO3) DISSOLUTION @@ -6997,10 +6988,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !--------------------------------------------------------------------------- ! Calculate calcite dissolution flux [mmolC m-2 day-1] - decayBenthos(4) = calc_diss_ben * LocBenthos(4) - - ! Update benthic calcite pool - LocBenthos(4) = LocBenthos(4) - decayBenthos(4) * dt_b + call update_benthos_decay(calc_diss_ben, 4) if (ciso) then @@ -7027,10 +7015,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & ! Calculate 13C remineralization flux ! Ignores isotopic fractionation during remineralization (alpha ≈ 1) - decayBenthos(5) = alpha_dcal_13 * decayRateBenC * LocBenthos(5) - - ! Update benthic 13C pool - LocBenthos(5) = LocBenthos(5) - decayBenthos(5) * dt_b + call update_benthos_decay(alpha_dcal_13 * decayRateBenC, 5) !======================================================================= ! Carbon-13 Calcite Dissolution @@ -7044,10 +7029,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !----------------------------------------------------------------------- ! Calculate 13C calcite dissolution flux - decayBenthos(6) = calc_diss_13 * LocBenthos(6) - - ! Update benthic 13C-calcite pool - LocBenthos(6) = LocBenthos(6) - decayBenthos(6) * dt_b + call update_benthos_decay(calc_diss_13, 6) if (ciso_14) then @@ -7069,10 +7051,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & ! Calculate 14C remineralization flux ! Ignores isotopic fractionation during remineralization - decayBenthos(7) = alpha_dcal_14 * decayRateBenC * LocBenthos(7) - - ! Update benthic 14C pool - LocBenthos(7) = LocBenthos(7) - decayBenthos(7) * dt_b + call update_benthos_decay(alpha_dcal_14 * decayRateBenC, 7) !=============================================================== ! 3.4 Carbon-14 Calcite Dissolution @@ -7086,10 +7065,7 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & !--------------------------------------------------------------- ! Calculate 14C calcite dissolution flux - decayBenthos(8) = calc_diss_14 * LocBenthos(8) - - ! Update benthic 14C-calcite pool - LocBenthos(8) = LocBenthos(8) - decayBenthos(8) * dt_b + call update_benthos_decay(calc_diss_14, 8) else !--------------------------------------------------------------- @@ -7111,6 +7087,16 @@ subroutine REcoM_sms(n,Nn,state,thick,recipthick,SurfSR,sms,Temp, Sali_depth & end do ! Main time loop ends +contains + +subroutine update_benthos_decay(rate, idx) + real(kind=8), intent(in) :: rate + integer , intent(in) :: idx + + decayBenthos(idx) = rate * LocBenthos(idx) + LocBenthos(idx) = LocBenthos(idx) - decayBenthos(idx) * dt_b +end subroutine update_benthos_decay + end subroutine REcoM_sms !-------------------------------------------------------------------------------