Refactoring ab3

config/namelist.dyn

@@ -21,5 +21,6 @@ use_wsplit   = .false. ! Switch for implicite/explicte splitting of vert. veloci
 wsplit_maxcfl= 1.0     ! maximum allowed CFL criteria in vertical (0.5 < w_max_cfl < 1.) 
                        ! in older FESOM it used to be w_exp_max=1.e-3
 ldiag_KE=.false.       ! activates energy diagnostics
+AB_order=2
 /
 

config/namelist.tra

@@ -32,6 +32,7 @@ gamma0_tra    = 0.0005     ! gammaX_tra are analogous to those in the dynamical
 gamma1_tra    = 0.0125
 gamma2_tra    = 0.
 i_vert_diff   =.true.
+AB_order      = 2
 /
 
 &tracer_phys

src/MOD_DYN.F90

@@ -51,8 +51,9 @@ MODULE MOD_DYN
 TYPE T_DYN
 !___________________________________________________________________________
     ! instant zonal merdional velocity & Adams-Bashfort rhs
-    real(kind=WP), allocatable, dimension(:,:,:):: uv, uv_rhs, uv_rhsAB, fer_uv
-
+    real(kind=WP), allocatable, dimension(:,:,:)   :: uv, uv_rhs, fer_uv
+    real(kind=WP), allocatable, dimension(:,:,:,:) :: uv_rhsAB
+    integer                                        :: AB_order=2
     ! horizontal velocities at nodes
     real(kind=WP), allocatable, dimension(:,:,:):: uvnode
 
@@ -108,12 +109,12 @@ MODULE MOD_DYN
     real(kind=WP), allocatable, dimension(:,:,:) :: ke_adv, ke_cor, ke_pre, ke_hvis, ke_vvis, ke_du2, ke_umean, ke_u2mean
     real(kind=WP), allocatable, dimension(:,:)   :: ke_wind, ke_drag
     ! same as above but multiplied by velocity. we need both for later computation of turbulent fluxes
-    real(kind=WP), allocatable, dimension(:,:,:) :: ke_adv_xVEL, ke_cor_xVEL, ke_pre_xVEL, ke_hvis_xVEL, ke_vvis_xVEL
-    real(kind=WP), allocatable, dimension(:,:)   :: ke_wind_xVEL, ke_drag_xVEL
-    real(kind=WP), allocatable, dimension(:,:)   :: ke_wrho         !we use pressure to compute (W*dens) as it appeares much easier to compute (P*dW) instead of (dP*w)
-    real(kind=WP), allocatable, dimension(:,:)   :: ke_dW, ke_Pfull !for later computation of turbulent fluxes from the term above
-    real(kind=WP), allocatable, dimension(:,:,:) :: ke_adv_AB, ke_cor_AB
-    real(kind=WP), allocatable, dimension(:,:,:) :: ke_rhs_bak
+    real(kind=WP), allocatable, dimension(:,:,:)   :: ke_adv_xVEL, ke_cor_xVEL, ke_pre_xVEL, ke_hvis_xVEL, ke_vvis_xVEL
+    real(kind=WP), allocatable, dimension(:,:)     :: ke_wind_xVEL, ke_drag_xVEL
+    real(kind=WP), allocatable, dimension(:,:)     :: ke_wrho         !we use pressure to compute (W*dens) as it appeares much easier to compute (P*dW) instead of (dP*w)
+    real(kind=WP), allocatable, dimension(:,:)     :: ke_dW, ke_Pfull !for later computation of turbulent fluxes from the term above
+    real(kind=WP), allocatable, dimension(:,:,:,:) :: ke_adv_AB, ke_cor_AB
+    real(kind=WP), allocatable, dimension(:,:,:)   :: ke_rhs_bak
     ! surface fields to compute APE generation
     real(kind=WP), allocatable, dimension(:)     :: ke_J, ke_D, ke_G, ke_D2, ke_n0, ke_JD, ke_GD, ke_swA, ke_swB
 

src/MOD_READ_BINARY_ARRAYS.F90

@@ -7,7 +7,7 @@ MODULE MOD_READ_BINARY_ARRAYS
 private
 public :: read_bin_array, read1d_int_static
 INTERFACE read_bin_array
-          MODULE PROCEDURE read1d_real, read1d_int, read1d_char, read2d_real, read2d_int, read3d_real, read3d_int
+          MODULE PROCEDURE read1d_real, read1d_int, read1d_char, read2d_real, read2d_int, read3d_real, read3d_int, read4d_real, read4d_int
 END INTERFACE
 contains
 subroutine read1d_real(arr, unit, iostat, iomsg)
@@ -113,5 +113,31 @@ subroutine read3d_int(arr, unit, iostat, iomsg)
     if (.not. allocated(arr)) allocate(arr(s1,s2,s3))
     read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3)
 end subroutine read3d_int
+
+subroutine read4d_real(arr, unit, iostat, iomsg)
+    real(kind=WP), intent(inout), allocatable :: arr(:,:,:,:)
+    integer,       intent(in)                 :: unit
+    integer,       intent(out)                :: iostat
+    character(*),  intent(inout)              :: iomsg
+    integer                                   :: s1, s2, s3, s4
+
+    read(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+    if ((s1==0) .or. (s2==0) .or. (s3==0) .or. (s4==0)) return
+    if (.not. allocated(arr)) allocate(arr(s1,s2,s3,s4))
+    read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3, 1:s4)
+end subroutine read4d_real
+
+subroutine read4d_int(arr, unit, iostat, iomsg)
+    integer, intent(inout), allocatable :: arr(:,:,:,:)
+    integer,       intent(in)           :: unit
+    integer,       intent(out)          :: iostat
+    character(*),  intent(inout)        :: iomsg
+    integer                             :: s1, s2, s3, s4
+
+    read(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+    if ((s1==0) .or. (s2==0) .or. (s3==0) .or. (s4==0)) return
+    if (.not. allocated(arr)) allocate(arr(s1,s2,s3,s4))
+    read(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3, 1:s4)
+end subroutine read4d_int
 end module MOD_READ_BINARY_ARRAYS
 !==========================================================

src/MOD_TRACER.F90

@@ -8,14 +8,18 @@ MODULE MOD_TRACER
 SAVE
 
 TYPE T_TRACER_DATA
-real(kind=WP), allocatable, dimension(:,:)  :: values, valuesAB  ! instant values & Adams-Bashfort interpolation
-logical                                     :: smooth_bh_tra=.false.
-real(kind=WP)                               :: gamma0_tra, gamma1_tra, gamma2_tra
-logical                                     :: i_vert_diff =.false.
-character(20)                               :: tra_adv_hor, tra_adv_ver, tra_adv_lim ! type of the advection scheme for this tracer
-real(kind=WP)                               :: tra_adv_ph  = 1.  ! a parameter to be used in horizontal advection (for MUSCL it is the fraction of fourth-order contribution in the solution)
-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                                     :: ID
+real(kind=WP), allocatable, dimension(:,:)    :: values    ! instant values
+real(kind=WP), allocatable, dimension(:,:)    :: valuesAB  ! Adams-Bashfort interpolation
+real(kind=WP), allocatable, dimension(:,:,:)  :: valuesold ! previous timesteps
+
+logical                                       :: smooth_bh_tra=.false.
+real(kind=WP)                                 :: gamma0_tra, gamma1_tra, gamma2_tra
+logical                                       :: i_vert_diff =.false.
+character(20)                                 :: tra_adv_hor, tra_adv_ver, tra_adv_lim ! type of the advection scheme for this tracer
+real(kind=WP)                                 :: tra_adv_ph  = 1.  ! a parameter to be used in horizontal advection (for MUSCL it is the fraction of fourth-order contribution in the solution)
+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
 
 contains
 procedure WRITE_T_TRACER_DATA
@@ -95,6 +99,7 @@ subroutine WRITE_T_TRACER_DATA(tdata, unit)
     character(len=1024)                  :: iomsg
 
     call write_bin_array(tdata%values,   unit, iostat, iomsg)
+    call write_bin_array(tdata%valuesold,unit, iostat, iomsg)
     call write_bin_array(tdata%valuesAB, unit, iostat, iomsg)
     write(unit, iostat=iostat, iomsg=iomsg) tdata%smooth_bh_tra
     write(unit, iostat=iostat, iomsg=iomsg) tdata%gamma0_tra
@@ -117,8 +122,9 @@ subroutine READ_T_TRACER_DATA(tdata, unit)
     integer                              :: iostat
     character(len=1024)                  :: iomsg
 
-    call read_bin_array(tdata%values,   unit, iostat, iomsg)
-    call read_bin_array(tdata%valuesAB, unit, iostat, iomsg)
+    call read_bin_array(tdata%values,    unit, iostat, iomsg)
+    call read_bin_array(tdata%valuesold,unit, iostat, iomsg)
+    call read_bin_array(tdata%valuesAB,  unit, iostat, iomsg)
     read(unit, iostat=iostat, iomsg=iomsg) tdata%smooth_bh_tra
     read(unit, iostat=iostat, iomsg=iomsg) tdata%gamma0_tra
     read(unit, iostat=iostat, iomsg=iomsg) tdata%gamma1_tra

src/MOD_WRITE_BINARY_ARRAYS.F90

@@ -7,7 +7,7 @@ MODULE MOD_WRITE_BINARY_ARRAYS
 private
 public :: write_bin_array, write1d_int_static
 INTERFACE write_bin_array
-          MODULE PROCEDURE write1d_real, write1d_int, write1d_char, write2d_real, write2d_int, write3d_real, write3d_int
+          MODULE PROCEDURE write1d_real, write1d_int, write1d_char, write2d_real, write2d_int, write3d_real, write3d_int, write4d_real, write4d_int
 END INTERFACE
 contains
 
@@ -155,5 +155,51 @@ subroutine write3d_int(arr, unit, iostat, iomsg)
        write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3
     end if
 end subroutine write3d_int
+subroutine write4d_real(arr, unit, iostat, iomsg)
+    real(kind=WP), intent(in), allocatable :: arr(:,:,:,:)
+    integer,       intent(in)              :: unit
+    integer,       intent(out)             :: iostat
+    character(*),  intent(inout)           :: iomsg
+    integer                                :: s1, s2, s3, s4
+
+    if (allocated(arr)) then
+       s1=size(arr, 1)
+       s2=size(arr, 2)
+       s3=size(arr, 3)
+       s4=size(arr, 4)
+       write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+       write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3, 1:s4)
+    else
+       s1=0
+       s2=0
+       s3=0
+       s4=0
+       write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+    end if
+end subroutine write4d_real
+
+subroutine write4d_int(arr, unit, iostat, iomsg)
+    integer,       intent(in), allocatable :: arr(:,:,:,:)
+    integer,       intent(in)              :: unit
+    integer,       intent(out)             :: iostat
+    character(*),  intent(inout)           :: iomsg
+    integer                                :: s1, s2, s3, s4
+
+    if (allocated(arr)) then
+       s1=size(arr, 1)
+       s2=size(arr, 2)
+       s3=size(arr, 3)
+       s4=size(arr, 4)
+       write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+       write(unit, iostat=iostat, iomsg=iomsg) arr(1:s1, 1:s2, 1:s3, 1:s4)
+    else
+       s1=0
+       s2=0
+       s3=0
+       s4=0
+       write(unit, iostat=iostat, iomsg=iomsg) s1, s2, s3, s4
+    end if
+end subroutine write4d_int
+
 end module MOD_WRITE_BINARY_ARRAYS
 !==========================================================

src/io_blowup.F90

@@ -115,8 +115,8 @@ subroutine ini_blowup_io(year, ice, dynamics, tracers, partit, mesh)
 		call def_variable(bid, 'v'			, (/nl-1, elem2D/)	, 'meridional velocity', 'm/s', dynamics%uv(2,:,:));
 		call def_variable(bid, 'u_rhs'			, (/nl-1, elem2D/)	, 'zonal velocity', 'm/s', dynamics%uv_rhs(1,:,:));
 		call def_variable(bid, 'v_rhs'			, (/nl-1, elem2D/)	, 'meridional velocity', 'm/s', dynamics%uv_rhs(2,:,:));
-		call def_variable(bid, 'urhs_AB'	, (/nl-1, elem2D/)	, 'Adams–Bashforth for u', 'm/s', dynamics%uv_rhsAB(1,:,:));
-		call def_variable(bid, 'vrhs_AB'	, (/nl-1, elem2D/)	, 'Adams–Bashforth for v', 'm/s', dynamics%uv_rhsAB(2,:,:));
+		call def_variable(bid, 'urhs_AB'	, (/nl-1, elem2D/)	, 'Adams-Bashforth for u', 'm/s', dynamics%uv_rhsAB(1,1,:,:));
+		call def_variable(bid, 'vrhs_AB'	, (/nl-1, elem2D/)	, 'Adams-Bashforth for v', 'm/s', dynamics%uv_rhsAB(1,2,:,:));
 		call def_variable(bid, 'zbar_n_bot' , (/nod2D/)			, 'node bottom depth', 'm', zbar_n_bot);
 		call def_variable(bid, 'zbar_e_bot' , (/elem2d/)		, 'elem bottom depth', 'm', zbar_e_bot);
 		call def_variable(bid, 'bottom_node_thickness' , (/nod2D/)			, 'node bottom thickness', 'm', bottom_node_thickness);

src/io_restart.F90

@@ -76,10 +76,14 @@ subroutine ini_ocean_io(year, dynamics, tracers, partit, mesh)
 #endif
   call oce_files%def_elem_var('u', 'zonal velocity',        'm/s', dynamics%uv(1,:,:), mesh, partit)
   call oce_files%def_elem_var('v', 'meridional velocity',   'm/s', dynamics%uv(2,:,:), mesh, partit)
-  call oce_files%def_elem_var('urhs_AB', 'Adams–Bashforth for u', 'm/s', dynamics%uv_rhsAB(1,:,:), mesh, partit)
-  call oce_files%def_elem_var('vrhs_AB', 'Adams–Bashforth for v', 'm/s', dynamics%uv_rhsAB(2,:,:), mesh, partit)
+  call oce_files%def_elem_var('urhs_AB', 'Adams-Bashforth for u (n-1 for AB2 and n-2 for AB3)', 'm/s', dynamics%uv_rhsAB(1,1,:,:), mesh, partit)
+  call oce_files%def_elem_var('vrhs_AB', 'Adams-Bashforth for v (n-1 for AB2 and n-2 for AB3)', 'm/s', dynamics%uv_rhsAB(1,2,:,:), mesh, partit)
+  if (dynamics%AB_order==3) then
+       call oce_files%def_elem_var_optional('urhs_AB3', 'Adams-Bashforth for u (n-1) for AB3', 'm/s', dynamics%uv_rhsAB(2,1,:,:), mesh, partit)
+       call oce_files%def_elem_var_optional('vrhs_AB3', 'Adams-Bashforth for v (n-1) for AB3', 'm/s', dynamics%uv_rhsAB(2,2,:,:), mesh, partit)
+  end if
 
-  !___Save restart variables for TKE and IDEMIX_________________________________
+!___Save restart variables for TKE and IDEMIX_________________________________
 !   if (trim(mix_scheme)=='cvmix_TKE' .or. trim(mix_scheme)=='cvmix_TKE+IDEMIX') then
   if (mix_scheme_nmb==5 .or. mix_scheme_nmb==56) then
         call oce_files%def_node_var_optional('tke', 'Turbulent Kinetic Energy', 'm2/s2', tke(:,:), mesh, partit)
@@ -109,8 +113,11 @@ subroutine ini_ocean_io(year, dynamics, tracers, partit, mesh)
          units='none'
      END SELECT
      call oce_files%def_node_var(trim(trname), trim(longname), trim(units), tracers%data(j)%values(:,:), mesh, partit)
-     longname=trim(longname)//', Adams–Bashforth'
-     call oce_files%def_node_var(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:), mesh, partit)
+     longname=trim(longname)//', Adams-Bashforth'
+     call oce_files%def_node_var(trim(trname)//'_AB', trim(longname), trim(units), tracers%data(j)%valuesAB(:,:),    mesh, partit)
+     call oce_files%def_node_var_optional(trim(trname)//'_M1', trim(longname), trim(units), tracers%data(j)%valuesold(1,:,:), mesh, partit)
+     if (tracers%data(j)%AB_order==3) &
+     call oce_files%def_node_var_optional(trim(trname)//'_M2', trim(longname), trim(units), tracers%data(j)%valuesold(2,:,:), mesh, partit)
   end do
   call oce_files%def_node_var('w', 'vertical velocity', 'm/s', dynamics%w, mesh, partit)
   call oce_files%def_node_var('w_expl', 'vertical velocity', 'm/s', dynamics%w_e, mesh, partit)

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(20)
+    type(restart_file_type), public :: files(22)
     integer, public :: nfiles = 0 ! todo: allow dynamically allocated size without messing with shallow copied pointers
   contains
     generic, public :: def_node_var => def_node_var_2d, def_node_var_3d

src/oce_ale_tracer.F90

@@ -148,7 +148,7 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh)
     type(t_partit), intent(inout), target    :: partit
     type(t_mesh)  , intent(in)   , target    :: mesh
     !___________________________________________________________________________
-    integer                                  :: tr_num, node, elem, nzmax, nzmin
+    integer                                  :: i, tr_num, node, elem, nzmax, nzmin
     !___________________________________________________________________________
     ! pointer on necessary derived types
     real(kind=WP), dimension(:,:,:), pointer :: UV, fer_UV
@@ -227,14 +227,6 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh)
         do node=1, myDim_nod2d
            tracers%work%del_ttf(:, node)=tracers%work%del_ttf(:, node)+tracers%work%del_ttf_advhoriz(:, node)+tracers%work%del_ttf_advvert(:, node)
         end do
-!$OMP END PARALLEL DO
-           !___________________________________________________________________________
-           ! AB is not needed after the advection step. Initialize it with the current tracer before it is modified.
-           ! call init_tracers_AB at the beginning of this loop will compute AB for the next time step then.
-!$OMP PARALLEL DO
-        do node=1, myDim_nod2d+eDim_nod2D
-           tracers%data(tr_num)%valuesAB(:, node)=tracers%data(tr_num)%values(:, node) !DS: check that this is the right place!
-        end do
 !$OMP END PARALLEL DO
         ! diffuse tracers
         if (flag_debug .and. mype==0)  print *, achar(27)//'[37m'//'         --> call diff_tracers_ale'//achar(27)//'[0m'