Skip to content

Add memory-light online streamfunction mode (lbas-indexed) #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hdrake wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Add memory-light online streamfunction mode (lbas-indexed) #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6837d34
Generalize kill-zone/tracer limits; size streamfunction zones adaptively
hdrake Jun 22, 2026
ba720ee
Fix _rerun.csv file handling so rerun can re-read its labels
hdrake Jun 23, 2026
daf0eec
Add memory-light online streamfunction mode (lbas-indexed)
hdrake Jun 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions projects/AVISO/kill_zones.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ SUBROUTINE kill_zones
! Exit domain defined by the boxes [ienw, iene]x[jens, jenn]
! in the namelist
CASE(1)
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1
Expand Down Expand Up @@ -66,7 +66,7 @@ SUBROUTINE kill_zones
END DO

! Next the geographical killing zone
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1 + numexit
Expand Down
4 changes: 2 additions & 2 deletions projects/IFS/kill_zones.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ SUBROUTINE kill_zones
! Exit domain defined by the boxes [ienw, iene]x[jens, jenn]
! in the namelist
CASE(1)
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1
Expand Down Expand Up @@ -66,7 +66,7 @@ SUBROUTINE kill_zones
END DO

! Next the geographical killing zone
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1 + numexit
Expand Down
4 changes: 2 additions & 2 deletions projects/NEMO/kill_zones.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ SUBROUTINE kill_zones
! Exit domain defined by the boxes [ienw, iene]x[jens, jenn]
! in the namelist
CASE(1)
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1
Expand Down Expand Up @@ -66,7 +66,7 @@ SUBROUTINE kill_zones
END DO

! Next the geographical killing zone
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1 + numexit
Expand Down
4 changes: 2 additions & 2 deletions projects/ROMS/kill_zones.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ SUBROUTINE kill_zones
! Exit domain defined by the boxes [ienw, iene]x[jens, jenn]
! in the namelist
CASE(1)
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1
Expand Down Expand Up @@ -66,7 +66,7 @@ SUBROUTINE kill_zones
END DO

! Next the geographical killing zone
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1 + numexit
Expand Down
4 changes: 2 additions & 2 deletions projects/Theoretical/kill_zones.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ SUBROUTINE kill_zones
! Exit domain defined by the boxes [ienw, iene]x[jens, jenn]
! in the namelist
CASE(1)
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1
Expand Down Expand Up @@ -66,7 +66,7 @@ SUBROUTINE kill_zones
END DO

! Next the geographical killing zone
DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)
IF(ienw(nexit) <= x1 .AND. x1 <= iene(nexit) .AND. &
jens(nexit) <= y1 .AND. y1 <= jenn(nexit) ) THEN
nend = nexit +1 + numexit
Expand Down
2 changes: 1 addition & 1 deletion src/mod_diffusion.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ SUBROUTINE kill_zones_diffusion(xb, xa, yb, ya)
REAL(DP), INTENT(INOUT) :: xb, xa, yb, ya
INTEGER :: nexit

DO nexit = 1, 10
DO nexit = 1, SIZE(ienw)

! Latitude band
IF (jens(nexit) == jenn(nexit)) THEN
Expand Down
25 changes: 23 additions & 2 deletions src/mod_init.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,10 +69,13 @@ SUBROUTINE init_namelist()
namelist /INIT_TRACERS_SEEDING/ tracer0min, tracer0max
namelist /INIT_KILLZONES/ timax, l_nosurface, exitType, ienw, iene, jens, jenn, &
tracerchoice, tracere, maxormin
namelist /INIT_POSTPROCESS/ l_psi, l_offline, dirpsi, xyflux, &
namelist /INIT_POSTPROCESS/ l_psi, l_offline, l_psi_rerun, dirpsi, xyflux, &
l_divergence, divconst
namelist /INIT_ACTIVE/ l_diffusion, ah, av

! Loop index for scanning geographic killing-zone slots
INTEGER :: izone

! Read namelist
OPEN (8,file='namelist.in', &
& status='OLD', delim='APOSTROPHE')
Expand All @@ -91,6 +94,24 @@ SUBROUTINE init_namelist()
READ (8,nml=INIT_ACTIVE)
CLOSE(8)

! Default the output-file prefix here, before it is first used.
! read_rerun runs before open_outfiles, so applying the default only
! inside open_outfiles left read_rerun looking for the wrong filename.
IF (TRIM(outDataFile) == '') outDataFile = 'TRACMASS'

! Detect the highest occupied geographic killing-zone slot from the
! namelist, before reverse()/zeroindx mutate the arrays below. Boxes
! live in up to 10 slots [ienw,iene]x[jens,jenn]; kill_zones sets
! nend = slot+1, so the slot INDEX (not the count) sets the largest
! geographic lbas. Undefined slots are all-zero (initialised in
! mod_vars). Subdomain wall zones and the final maxlbas are added in
! init_subdomain, once exitType and all killing zones are finalised.
ngeozones = 0
DO izone = 1, SIZE(ienw)
IF (ienw(izone) /= 0 .OR. iene(izone) /= 0 .OR. &
jens(izone) /= 0 .OR. jenn(izone) /= 0) ngeozones = izone
END DO

! If input data is on a A grid
#ifdef A_grid
jmt = jmt-1
Expand Down Expand Up @@ -174,7 +195,7 @@ SUBROUTINE reverse()
INTEGER :: ii

IF (griddir(2) == -1) THEN
DO ii = 1, 10
DO ii = 1, SIZE(jenn)
IF (isec == 2) THEN
jenn(ii) = jmt - jenn(ii) ! Meridional reverse
jens(ii) = jmt - jens(ii) ! Meridional reverse
Expand Down
45 changes: 31 additions & 14 deletions src/mod_stream.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ SUBROUTINE compute_stream()
IF (l_tracers) CALL open_outstream('yr')
IF (l_tracers) CALL open_outstream('rr')

DO ilvar1 = 1, 21
DO ilvar1 = 1, maxlbas

psi_xy(:,:) = 0.; psi_xz(:,:) = 0.; psi_yz(:,:) = 0.
IF (l_tracers) THEN
Expand All @@ -59,8 +59,10 @@ SUBROUTINE compute_stream()
psi_rr(:,:) = 0.
END IF

! For online calculation - adding fluxes
IF (l_offline .EQV. .FALSE.) THEN
! Legacy per-trajectory online calculation: sum each trajectory's
! fluxes into slot 0 by killing zone. Offline and the lbas-indexed
! rerun mode already have fluxes binned by zone, so they skip this.
IF ((l_offline .EQV. .FALSE.) .AND. (l_psi_rerun .EQV. .FALSE.)) THEN

! Cleaning of fluxes
fluxes_xy(:,:,0) = 0.d0; fluxes_xz(:,:,0) = 0.d0; fluxes_yz(:,:,0) = 0.d0
Expand Down Expand Up @@ -94,8 +96,10 @@ SUBROUTINE compute_stream()
END DO intrajLoop
END IF

! Pre-binned fluxes (offline / rerun) read zone ilvar1 directly;
! legacy online reads the summed slot 0.
ilvar3 = 0
IF (l_offline) ilvar3 = ilvar1
IF (l_offline .OR. l_psi_rerun) ilvar3 = ilvar1

IF (dirpsi(ilvar1) == 1) THEN
DO ilvar2 = 2, MAX(imtdom,jmtdom,km,resolution)
Expand Down Expand Up @@ -170,8 +174,11 @@ SUBROUTINE init_stream()

INTEGER :: index

index = 21
IF (l_offline .EQV. .FALSE.) THEN
! Legacy per-trajectory online accumulation needs one slot per
! trajectory (lbas unknown during the run); offline and the
! lbas-indexed rerun mode size to the killing-zone count maxlbas.
index = maxlbas
IF ((l_offline .EQV. .FALSE.) .AND. (l_psi_rerun .EQV. .FALSE.)) THEN
index = ntracmax
END IF

Expand Down Expand Up @@ -208,9 +215,19 @@ SUBROUTINE update_fluxes(indx1, indx2, dir, psicase, indt1, indt2)
INTEGER, INTENT(IN) :: indx1, indx2, dir
INTEGER, INTENT(IN), OPTIONAL :: indt1, indt2
CHARACTER(LEN=2), INTENT(IN) :: psicase
INTEGER :: index1, index2, indm1, indm2
INTEGER :: index1, index2, indm1, indm2, ipsi
REAL(DP) :: slope

! Select the flux-array slot. Legacy online accumulation stores one
! slot per trajectory (lbas not yet known); the lbas-indexed rerun
! mode stores directly into the trajectory's killing-zone slot.
IF (l_psi_rerun) THEN
ipsi = trajectories(ntrac)%lbas
IF (ipsi <= 0) RETURN
ELSE
ipsi = ntrac
END IF

!
index1 = indx1
index2 = indx2
Expand All @@ -230,15 +247,15 @@ SUBROUTINE update_fluxes(indx1, indx2, dir, psicase, indt1, indt2)
END IF

! Geographical streamfunctions
IF (psicase=='xy') fluxes_xy(index1,index2,ntrac) = fluxes_xy(index1,index2,ntrac) + dir*subvol
IF (psicase=='xz') fluxes_xz(index1,index2,ntrac) = fluxes_xz(index1,index2,ntrac) + dir*subvol
IF (psicase=='yz') fluxes_yz(index1,index2,ntrac) = fluxes_yz(index1,index2,ntrac) + dir*subvol
IF (psicase=='xy') fluxes_xy(index1,index2,ipsi) = fluxes_xy(index1,index2,ipsi) + dir*subvol
IF (psicase=='xz') fluxes_xz(index1,index2,ipsi) = fluxes_xz(index1,index2,ipsi) + dir*subvol
IF (psicase=='yz') fluxes_yz(index1,index2,ipsi) = fluxes_yz(index1,index2,ipsi) + dir*subvol

! Geographical + tracer
IF (psicase=='xr' .AND. PRESENT(indt1)) THEN
fluxes_xr(index1,index2,ntrac,indt1) = fluxes_xr(index1,index2,ntrac,indt1) + dir*subvol
fluxes_xr(index1,index2,ipsi,indt1) = fluxes_xr(index1,index2,ipsi,indt1) + dir*subvol
ELSE IF (psicase=='yr' .AND. PRESENT(indt1)) THEN
fluxes_yr(index1,index2,ntrac,indt1) = fluxes_yr(index1,index2,ntrac,indt1) + dir*subvol
fluxes_yr(index1,index2,ipsi,indt1) = fluxes_yr(index1,index2,ipsi,indt1) + dir*subvol
END IF

! Tracer-tracer equation
Expand All @@ -249,15 +266,15 @@ SUBROUTINE update_fluxes(indx1, indx2, dir, psicase, indt1, indt2)
slope = (FLOAT(indm1)-FLOAT(index1))*(FLOAT(indt2)-FLOAT(indt1))/(FLOAT(index2)-FLOAT(index1))
indm2 = NINT( slope + indt1)

fluxes_rr(indm1,indm2,ntrac) = fluxes_rr(indm1,indm2,ntrac) + subvol
fluxes_rr(indm1,indm2,ipsi) = fluxes_rr(indm1,indm2,ipsi) + subvol
END DO

DO indm1 = index2, index1-1

slope = (FLOAT(indm1)-FLOAT(index1))*(FLOAT(indt2)-FLOAT(indt1))/(FLOAT(index2)-FLOAT(index1))
indm2 = NINT( slope + indt1)

fluxes_rr(indm1,indm2,ntrac) = fluxes_rr(indm1,indm2,ntrac) - subvol
fluxes_rr(indm1,indm2,ipsi) = fluxes_rr(indm1,indm2,ipsi) - subvol
END DO

END IF
Expand Down
79 changes: 63 additions & 16 deletions src/mod_subdomain.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ MODULE mod_subdomain

USE mod_domain
USE mod_grid
USE mod_tracervars
USE mod_postprocessvars

IMPLICIT NONE

Expand All @@ -27,6 +29,15 @@ SUBROUTINE init_subdomain()
!
! --------------------------------------------------

INTEGER :: ntracerkill ! number of tracer-based killing zones
INTEGER :: subgeomax ! highest subdomain wall killing-zone slot
INTEGER :: ngeo ! number of geographic killing-zone slots

subgeomax = 0
! Subdomain wall killing zones occupy the LAST 4 geographic slots,
! leaving the lower slots for user-defined namelist zones.
ngeo = SIZE(ienw)

! Make sure killing zones are on
IF (exitType==2 .AND. l_subdom) THEN
exitType = 3 ! Include both thermodynamic and geographical killing zones
Expand Down Expand Up @@ -55,22 +66,25 @@ SUBROUTINE init_subdomain()

! These killing zones will not be activated for hemispheric cap subdomains
IF ((iperio == 1 .AND. imt == imtdom .AND. jmaxdom == 1) .EQV. .FALSE.) THEN
! 7 - south wall
ienw(7) = -1; iene(7) = imtdom + 1; jens(7) = jmindom + 1; jenn(7) = jmindom + 1
! south wall (4th-from-last slot)
ienw(ngeo-3) = -1; iene(ngeo-3) = imtdom + 1; jens(ngeo-3) = jmindom + 1; jenn(ngeo-3) = jmindom + 1
subgeomax = ngeo-3
END IF

IF ((iperio == 1 .AND. imt == imtdom .AND. jmaxdom == jmtdom) .EQV. .FALSE.) THEN
! 8 - north wall
ienw(8) = -1; iene(8) = imtdom + 1; jens(8) = jmaxdom - 1; jenn(8) = jmaxdom - 1
! north wall (3rd-from-last slot)
ienw(ngeo-2) = -1; iene(ngeo-2) = imtdom + 1; jens(ngeo-2) = jmaxdom - 1; jenn(ngeo-2) = jmaxdom - 1
subgeomax = ngeo-2
END IF

! These killing zones will not be activated if iperio = 1
! and imindom = 1 and imaxdom = imt
IF ((iperio == 1 .AND. imt == imtdom) .EQV. .FALSE.) THEN
! 9 - east wall
ienw(9) = imindom + 1; iene(9) = imindom + 1; jens(9) = -1; jenn(9) = jmtdom + 1
! 10 - west wall
ienw(10) = imaxdom - 1; iene(10) = imaxdom - 1; jens(10) = - 1; jenn(10) = jmtdom + 1
! east wall (2nd-from-last slot)
ienw(ngeo-1) = imindom + 1; iene(ngeo-1) = imindom + 1; jens(ngeo-1) = -1; jenn(ngeo-1) = jmtdom + 1
! west wall (last slot)
ienw(ngeo) = imaxdom - 1; iene(ngeo) = imaxdom - 1; jens(ngeo) = - 1; jenn(ngeo) = jmtdom + 1
subgeomax = ngeo
END IF

! Redefine the killing zones in the new reference system
Expand All @@ -89,14 +103,15 @@ SUBROUTINE init_subdomain()
jmt = jmaxdom - jmindom + 1

! The last 4 kill zones are reserved to the Subdomain
! 7 - south wall
ienw(7) = -1; iene(7) = imt + 1; jens(7) = jmindom + 1; jenn(7) = jmindom + 1
! 8 - north wall
ienw(8) = -1; iene(8) = imt + 1; jens(8) = jmaxdom - 1; jenn(8) = jmaxdom - 1
! 9 - east wall
ienw(9) = imindom + 1; iene(9) = imindom + 1; jens(9) = -1; jenn(9) = jmt + 1
! 10 - west wall
ienw(10) = imaxdom - 1; iene(10) = imaxdom - 1; jens(10) = - 1; jenn(10) = jmt + 1
! south wall (4th-from-last slot)
ienw(ngeo-3) = -1; iene(ngeo-3) = imt + 1; jens(ngeo-3) = jmindom + 1; jenn(ngeo-3) = jmindom + 1
! north wall (3rd-from-last slot)
ienw(ngeo-2) = -1; iene(ngeo-2) = imt + 1; jens(ngeo-2) = jmaxdom - 1; jenn(ngeo-2) = jmaxdom - 1
! east wall (2nd-from-last slot)
ienw(ngeo-1) = imindom + 1; iene(ngeo-1) = imindom + 1; jens(ngeo-1) = -1; jenn(ngeo-1) = jmt + 1
! west wall (last slot)
ienw(ngeo) = imaxdom - 1; iene(ngeo) = imaxdom - 1; jens(ngeo) = - 1; jenn(ngeo) = jmt + 1
subgeomax = ngeo

! Redefine the killing zones in the new reference system
ienw = ienw - imindom + 1; iene = iene - imindom + 1;
Expand All @@ -117,6 +132,38 @@ SUBROUTINE init_subdomain()

END IF

! Finalise maxlbas now that exitType (possibly promoted above) and all
! killing zones are known. Size the per-zone streamfunction/summary
! arrays to the actual run. nend (lbas) encoding in the project
! kill_zones.F90 routines:
! nend = 0 -> time limit
! nend = 1 -> reaching the surface
! nend = nexit+1 -> geographic killing zone nexit
! nend = nexit+1+numexit -> tracer killing zone (exitType 3)
! ngeozones (highest geographic slot) was seeded from the namelist in
! init_namelist; bump it for the subdomain wall zones added above.
! Tracer zones are counted via the 999 sentinel default.
ngeozones = MAX(ngeozones, subgeomax)
ntracerkill = COUNT(tracerchoice /= 999)

SELECT CASE (exitType)
CASE (1) ! geographic killing zones only
maxlbas = 1 + ngeozones
CASE (2) ! tracer-based killing zones only
maxlbas = 1 + ntracerkill
CASE (3) ! tracer + geographic killing zones
maxlbas = 1 + ntracerkill + ngeozones
CASE DEFAULT ! exitType 4 (hard coded) or unset: full safety
maxlbas = MAXZONES
END SELECT

IF (maxlbas > MAXZONES) THEN
PRINT*, 'ERROR: number of killing zones (maxlbas =', maxlbas, &
') exceeds MAXZONES =', MAXZONES
PRINT*, 'Increase MAXZONES in mod_precdef (src/mod_vars.F90).'
STOP
END IF

END SUBROUTINE init_subdomain

SUBROUTINE update_subindex(ji,jj)
Expand Down
4 changes: 2 additions & 2 deletions src/mod_tracers.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ MODULE mod_tracers

IMPLICIT NONE

INTEGER, DIMENSION(10) :: numtracerarray = 0
INTEGER, DIMENSION(MAXTRACERS) :: numtracerarray = 0

PRIVATE :: tracers_default

Expand All @@ -39,7 +39,7 @@ SUBROUTINE init_tracer()

! Calculate the number of tracers
WHERE (tracername==' ') numtracerarray = 1
numtracers = 10 - SUM(numtracerarray)
numtracers = SIZE(numtracerarray) - SUM(numtracerarray)

! Allocate the tracer array
ALLOCATE(tracers(numtracers), tracervalue(numtracers))
Expand Down
Loading