Draft: Refactoring addsshsubcycl

config/namelist.dyn

@@ -22,5 +22,15 @@ wsplit_maxcfl= 1.0     ! maximum allowed CFL criteria in vertical (0.5 < w_max_c
                        ! in older FESOM it used to be w_exp_max=1.e-3
 ldiag_KE=.false.       ! activates energy diagnostics
 AB_order=2
+
+use_ssh_se_subcycl = .false.
+se_BTsteps     = 50   
+se_BTtheta     = 0.14   ! default: 0.14, 
+se_bottdrag    = .true.
+se_bdrag_si    = .true. ! bottomdrag semi-implicite/explicite
+se_visc        = .true.
+se_visc_gamma0 = 10  
+se_visc_gamma1 = 19500  !19500 (core2@32spd), 2750 (core2@72spd)
+se_visc_gamma2 = 0
 /
 

src/MOD_DYN.F90

@@ -42,7 +42,7 @@ MODULE MOD_DYN
 ! set main structure for dynamicss, contains viscosity options and parameters +
 ! option for momentum advection
 TYPE T_DYN
-!___________________________________________________________________________
+    !___________________________________________________________________________
     ! instant zonal merdional velocity & Adams-Bashfort rhs
     real(kind=WP), allocatable, dimension(:,:,:)   :: uv, uv_rhs, fer_uv
     real(kind=WP), allocatable, dimension(:,:,:,:) :: uv_rhsAB
@@ -55,6 +55,24 @@ MODULE MOD_DYN
 
     ! sea surface height arrays
     real(kind=WP), allocatable, dimension(:)    :: eta_n, d_eta, ssh_rhs, ssh_rhs_old
+
+    !___arrays for split explicite ssh computation______________________________
+    ! se_uvh...transport velocity, 
+    real(kind=WP), allocatable, dimension(:,:,:):: se_uvh 
+    !se_uv_rhs...vertical integral of transport velocity rhs, se_uvBT_4AB...
+    ! barotropic transport velocities (vertically integrated), contains actual 
+    ! timestep (1:2) and previous timestep (3:4) for adams-bashfort interpolation
+    real(kind=WP), allocatable, dimension(:,:)  :: se_uvBT_rhs, se_uvBT_4AB 
+
+    ! se_uvBT...barotropic trnasport velocities from barotropic time stepping
+    ! se_uvBT_theta...velocities for dissipative time stepping of thickness equation
+    ! UBTmean_mean... Mean BT velocity to trim 3D velocity in tracers
+    real(kind=WP), allocatable, dimension(:,:)  :: se_uvBT, se_uvBT_theta, se_uvBT_mean, se_uvBT_12 
+
+    ! array that are needed for viscosity and bottomdrag stabilization of 
+    ! split-expl subcycling method
+    real(kind=WP), allocatable, dimension(:,:)  :: se_uvBT_stab_hvisc 
+    real(kind=WP), allocatable, dimension(:)    :: se_uvBT_stab_bdrag
 
     ! LA: 2023-05-17 iceberg arrays
     real(kind=WP), allocatable, dimension(:)    :: eta_n_ib ! kh 18.03.21 additional array for asynchronous iceberg computations
@@ -99,7 +117,24 @@ MODULE MOD_DYN
     logical                                     :: use_wsplit    = .false.
     ! 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
-    real(kind=WP)                               :: wsplit_maxcfl= 1.0
+    real(kind=WP)                               :: wsplit_maxcfl = 1.0
+
+    ! switch between ssh computation, by solver or split explicite subcycling
+    ! use_ssh_se_subcycl = .false. --> solver
+    ! use_ssh_se_subcycl = .true.  --> split explicite subcycling
+    logical                                     :: use_ssh_se_subcycl = .false.
+
+    ! barotropic subcycling time-steps and dissipation parameter
+    integer                                     :: se_BTsteps     = 50
+    real(kind=WP)                               :: se_BTtheta     = 0.14_WP
+    logical                                     :: se_bottdrag    = .true.
+    logical                                     :: se_bdrag_si    = .true.
+    logical                                     :: se_visc        = .true.
+    real(kind=WP)                               :: se_visc_gamma0 = 10
+    real(kind=WP)                               :: se_visc_gamma1 = 2750
+    real(kind=WP)                               :: se_visc_gamma2 = 0
+
+    !___________________________________________________________________________
     ! energy diagnostic part: will be computed inside the model ("hard integration"):
     logical                                      :: ldiag_ke       = .true.
     ! different contributions to velocity change. will be computed inside the code.
@@ -229,7 +264,8 @@ subroutine WRITE_T_DYN(dynamics, unit, iostat, iomsg)
     write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_freeslip
     write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_wsplit
     write(unit, iostat=iostat, iomsg=iomsg) dynamics%wsplit_maxcfl
-
+    write(unit, iostat=iostat, iomsg=iomsg) dynamics%use_ssh_se_subcycl
+
     !___________________________________________________________________________
     call dynamics%solverinfo%WRITE_T_SOLVERINFO(unit)
 
@@ -249,6 +285,17 @@ subroutine WRITE_T_DYN(dynamics, unit, iostat, iomsg)
         call write_bin_array(dynamics%fer_w , unit, iostat, iomsg)
         call write_bin_array(dynamics%fer_uv, unit, iostat, iomsg)
     end if
+    if (dynamics%use_ssh_se_subcycl) then
+        call write_bin_array(dynamics%se_uvh       , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_rhs , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_4AB , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT      , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_theta, unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_mean , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_12 , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_stab_hvisc , unit, iostat, iomsg)
+        call write_bin_array(dynamics%se_uvBT_stab_bdrag , unit, iostat, iomsg)
+    end if 
 
 
 end subroutine WRITE_T_DYN
@@ -275,6 +322,7 @@ subroutine READ_T_DYN(dynamics, unit, iostat, iomsg)
     read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_freeslip
     read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_wsplit
     read(unit, iostat=iostat, iomsg=iomsg) dynamics%wsplit_maxcfl
+    read(unit, iostat=iostat, iomsg=iomsg) dynamics%use_ssh_se_subcycl
 
     !___________________________________________________________________________
     call dynamics%solverinfo%READ_T_SOLVERINFO(unit)
@@ -295,6 +343,17 @@ subroutine READ_T_DYN(dynamics, unit, iostat, iomsg)
         call read_bin_array(dynamics%fer_w     , unit, iostat, iomsg)
         call read_bin_array(dynamics%fer_uv    , unit, iostat, iomsg)
     end if
+    if (dynamics%use_ssh_se_subcycl) then
+        call read_bin_array(dynamics%se_uvh       , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_rhs , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_4AB , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT      , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_theta, unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_mean , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_12 , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_stab_hvisc , unit, iostat, iomsg)
+        call read_bin_array(dynamics%se_uvBT_stab_bdrag , unit, iostat, iomsg)
+    end if 
 
 end subroutine READ_T_DYN
 

src/associate_mesh_ass.h

@@ -57,14 +57,14 @@ Z_3d_n(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D)              => mesh%Z_3d_n(:,:)
 #if defined(__async_icebergs)
 Z_3d_n_ib(1:mesh%nl-1, 1:myDim_nod2D+eDim_nod2D)           => mesh%Z_3d_n_ib(:,:)
 #endif
-helem(1:mesh%nl-1, 1:myDim_elem2D)                         => mesh%helem(:,:)
-bottom_elem_thickness(1:myDim_elem2D)                      => mesh%bottom_elem_thickness(:)
+helem(1:mesh%nl-1, 1:myDim_elem2D+eDim_elem2D)             => mesh%helem(:,:)
+bottom_elem_thickness(1:myDim_elem2D+eDim_elem2D)          => mesh%bottom_elem_thickness(:)
 bottom_node_thickness(1:myDim_nod2D+eDim_nod2D)            => mesh%bottom_node_thickness(:)
 dhe(1:myDim_elem2D)                                        => mesh%dhe(:)
 hbar(1:myDim_nod2D+eDim_nod2D)                             => mesh%hbar(:)
 hbar_old(1:myDim_nod2D+eDim_nod2D)                         => mesh%hbar_old(:)
-!zbar_n(1:mesh%nl)                                          => mesh%zbar_n
-!Z_n(1:mesh%nl-1)                                           => mesh%Z_n
+!zbar_n(1:mesh%nl)                                         => mesh%zbar_n
+!Z_n(1:mesh%nl-1)                                          => mesh%Z_n
 zbar_n_bot(1:myDim_nod2D+eDim_nod2D)                       => mesh%zbar_n_bot(:)
 zbar_e_bot(1:myDim_elem2D+eDim_elem2D)                     => mesh%zbar_e_bot(:)
 zbar_n_srf(1:myDim_nod2D+eDim_nod2D)                       => mesh%zbar_n_srf(:)

src/gen_modules_diag.F90

@@ -80,9 +80,10 @@ module diagnostics
   logical                                       :: ldiag_extflds    =.false.
   logical                                       :: ldiag_ice        =.false.
   logical                                       :: ldiag_trflx      =.false.
-
-  namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, &
-                       ldiag_salt3D, ldiag_forc, ldiag_vorticity, ldiag_extflds, ldiag_trflx
+
+  namelist /diag_list/ ldiag_solver, lcurt_stress_surf, ldiag_curl_vel3, ldiag_Ri, & 
+                       ldiag_TurbFlux, ldiag_dMOC, ldiag_DVD, ldiag_salt3D, ldiag_forc, &
+                       ldiag_vorticity, ldiag_extflds, ldiag_trflx, ldiag_ice
 
   contains
 
@@ -181,69 +182,78 @@ subroutine diag_curl_vel3(mode, dynamics, partit, mesh)
 #include "associate_mesh_ass.h" 
     UV => dynamics%uv(:,:,:)
 
-    !=====================
+    !___________________________________________________________________________
     if (firstcall) then  !allocate the stuff at the first call
         allocate(curl_vel3(nl-1, myDim_nod2D+eDim_nod2D))
         firstcall=.false.
         if (mode==0) return
     end if
 
+    !___________________________________________________________________________
     curl_vel3=0.
-
-    DO ed=1,myDim_edge2D
+    do ed=1,myDim_edge2D
         enodes=edges(:,ed)
         el=edge_tri(:,ed)
+        !_______________________________________________________________________
         nl1=nlevels(el(1))-1
         nu1=ulevels(el(1))
         deltaX1=edge_cross_dxdy(1,ed)
         deltaY1=edge_cross_dxdy(2,ed)
-        nl2=0
-        nu2=0
+
+        !_______________________________________________________________________
         if (el(2)>0) then
             deltaX2=edge_cross_dxdy(3,ed)
             deltaY2=edge_cross_dxdy(4,ed)
             nl2=nlevels(el(2))-1
             nu2=ulevels(el(2))
-        end if   
-
-        nl12 = min(nl1,nl2)
-        nu12 = min(nu1,nu2)
-        DO nz=nu1,nu12-1
-            c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))
-            curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
-            curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
-        END DO
-        if (nu2>0) then
-            DO nz=nu2,nu12-1
+            nl12 = min(nl1,nl2)
+            nu12 = max(nu1,nu2)
+            !___________________________________________________________________
+            do nz=nu1,nu12-1
+                c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))
+                curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
+                curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
+            end do
+            do nz=nu2,nu12-1
                 c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2))
                 curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
                 curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
-            END DO
+            end do
+            !___________________________________________________________________
+            do nz=nu12,nl12
+                c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))- &
+                deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2))
+                curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
+                curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
+            end do
+            !___________________________________________________________________
+            do nz=nl12+1,nl1
+                c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))
+                curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
+                curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
+            end do
+            do nz=nl12+1,nl2
+                c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2))
+                curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
+                curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
+            end do
+        !_______________________________________________________________________    
+        else
+            do nz=nu1,nl1
+                c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))
+                curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
+                curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
+            end do
         end if
-        DO nz=nu12,nl12
-            c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))- &
-            deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2))
-            curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
-            curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
-        END DO
-        DO nz=nl12+1,nl1
-            c1=deltaX1*UV(1,nz,el(1))+deltaY1*UV(2,nz,el(1))
-            curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
-            curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
-        END DO
-        DO nz=nl12+1,nl2
-            c1= -deltaX2*UV(1,nz,el(2))-deltaY2*UV(2,nz,el(2))
-            curl_vel3(nz,enodes(1))=curl_vel3(nz,enodes(1))+c1
-            curl_vel3(nz,enodes(2))=curl_vel3(nz,enodes(2))-c1
-        END DO
-    END DO
-
-    DO n=1, myDim_nod2D
-        !!PS DO nz=1, nlevels_nod2D(n)-1
-        DO nz=ulevels_nod2D(n), nlevels_nod2D(n)-1
+    end do
+
+    !___________________________________________________________________________
+    do n=1, myDim_nod2D
+        do nz=ulevels_nod2D(n), nlevels_nod2D(n)-1
             curl_vel3(nz,n)=curl_vel3(nz,n)/areasvol(nz,n)
-        END DO
-    END DO
+        end do
+    end do
+
 end subroutine diag_curl_vel3
 ! ==============================================================
 ! 

src/io_blowup.F90

@@ -102,12 +102,27 @@ subroutine ini_blowup_io(year, ice, dynamics, tracers, partit, mesh)
 		!___Define the netCDF variables for 2D fields_______________________________
 		!___SSH_____________________________________________________________________
 		call def_variable(bid, 'eta_n'		, (/nod2D/)			, 'sea surface elevation', 'm', dynamics%eta_n);
-		call def_variable(bid, 'd_eta'		, (/nod2D/)			, 'change in ssh from solver', 'm', dynamics%d_eta);
 		!___ALE related fields______________________________________________________
 		call def_variable(bid, 'hbar'		, (/nod2D/)			, 'ALE surface elevation hbar_n+0.5', 'm', hbar);
 !!PS 		call def_variable(bid, 'hbar_old'	, (/nod2D/)			, 'ALE surface elevation hbar_n-0.5', 'm', hbar_old);
-		call def_variable(bid, 'ssh_rhs'	, (/nod2D/)			, 'RHS for the elevation', '?', dynamics%ssh_rhs);
-		call def_variable(bid, 'ssh_rhs_old', (/nod2D/)			, 'RHS for the elevation', '?', dynamics%ssh_rhs_old);
+        if (.not. dynamics%use_ssh_se_subcycl) then
+            call def_variable(bid, 'd_eta'		, (/nod2D/)			, 'change in ssh from solver', 'm', dynamics%d_eta);
+            call def_variable(bid, 'ssh_rhs'	, (/nod2D/)			, 'RHS for the elevation', '?', dynamics%ssh_rhs);
+            call def_variable(bid, 'ssh_rhs_old', (/nod2D/)			, 'RHS for the elevation', '?', dynamics%ssh_rhs_old);
+
+        else
+            call def_variable(bid, 'ubt_rhs'  , (/elem2D/), 'zonal RHS barotr. transp. equation' , '?'  , dynamics%se_uvBT_rhs(  1,:));
+            call def_variable(bid, 'vbt_rhs'  , (/elem2D/), 'merid. RHS barotr. transp. equation', '?'  , dynamics%se_uvBT_rhs(  2,:));
+            call def_variable(bid, 'ubt'	  , (/elem2D/), 'zonal barotr. transp.'              , '?'  , dynamics%se_uvBT(      1,:));
+            call def_variable(bid, 'vbt'	  , (/elem2D/), 'merid. barotr. transp.'             , '?'  , dynamics%se_uvBT(      2,:));
+            call def_variable(bid, 'ubt_theta', (/elem2D/), 'zonal barotr. theta term.'          , '?'  , dynamics%se_uvBT_theta(1,:));
+            call def_variable(bid, 'vbt_theta', (/elem2D/), 'merid. barotr. theta term'          , '?'  , dynamics%se_uvBT_theta(2,:));
+            call def_variable(bid, 'ubt_mean' , (/elem2D/), 'zonal barotr. mean term.'           , '?'  , dynamics%se_uvBT_mean( 1,:));
+            call def_variable(bid, 'vbt_mean' , (/elem2D/), 'merid. barotr. mean term'           , '?'  , dynamics%se_uvBT_mean( 2,:));
+            call def_variable(bid, 'uh'       , (/nl-1, elem2D/), 'zonal velocity'               , 'm/s', dynamics%se_uvh(1,:,:));
+            call def_variable(bid, 'vh'       , (/nl-1, elem2D/), 'meridional velocity'          , 'm/s', dynamics%se_uvh(2,:,:));
+		end if 
+
 		!___Define the netCDF variables for 3D fields_______________________________
 		call def_variable(bid, 'hnode'		, (/nl-1,  nod2D/)	, 'ALE stuff', '?', hnode);
 		call def_variable(bid, 'helem'		, (/nl-1, elem2D/)	, 'Element layer thickness', 'm/s', helem(:,:));