From 6837d34731d04915f722fa02ed1e23be84a62a5a Mon Sep 17 00:00:00 2001 From: hdrake Date: Mon, 22 Jun 2026 19:04:19 -0400 Subject: [PATCH] Generalize kill-zone/tracer limits; size streamfunction zones adaptively The streamfunction flux arrays and their per-zone loops were hard-coded to a fixed 21 kill zones, and the kill-zone (10) and tracer (10) capacities were scattered as bare literals across many files (including per-project kill_zones.F90). This generalizes those limits and sizes the streamfunction allocation to the zones a run actually uses. - Add named capacity constants in mod_precdef: MAXGEOZONES, MAXTRACERS, and MAXZONES = 1 + MAXTRACERS + MAXGEOZONES, used to dimension the kill-zone and tracer configuration arrays. - Compute maxlbas (the actual streamfunction zone count) from the run: namelist geographic zones + subdomain wall zones + tracer kill zones, finalised in init_subdomain once exitType and all zones are known. The flux arrays (offline) and the compute_stream loop now size to maxlbas instead of a fixed 21. - Derive runtime loop bounds and the subdomain "last 4" wall slots from SIZE() of the arrays, removing the hard-coded slot numbers (7-10) in mod_subdomain and the 1,10 loops in kill_zones.F90 (x5), mod_diffusion, and mod_init. - Initialise ienw/iene/jens/jenn to 0 so undefined slots are well-defined. Fixes a latent out-of-bounds when the configured zone count differs from 21, and an under-allocation when an exitType=2 run uses a subdomain (which promotes exitType to 3 and adds geographic walls). Behaviour is otherwise unchanged: streamfunction output now contains the actual number of zone blocks, verified byte-identical to the previous output for the real zones. Co-Authored-By: Claude Opus 4.8 --- projects/AVISO/kill_zones.F90 | 4 +- projects/IFS/kill_zones.F90 | 4 +- projects/NEMO/kill_zones.F90 | 4 +- projects/ROMS/kill_zones.F90 | 4 +- projects/Theoretical/kill_zones.F90 | 4 +- src/mod_diffusion.F90 | 2 +- src/mod_init.F90 | 18 ++++++- src/mod_stream.F90 | 4 +- src/mod_subdomain.F90 | 79 +++++++++++++++++++++++------ src/mod_tracers.F90 | 4 +- src/mod_vars.F90 | 30 +++++++---- 11 files changed, 115 insertions(+), 42 deletions(-) diff --git a/projects/AVISO/kill_zones.F90 b/projects/AVISO/kill_zones.F90 index 25193ea..ee87a9c 100644 --- a/projects/AVISO/kill_zones.F90 +++ b/projects/AVISO/kill_zones.F90 @@ -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 @@ -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 diff --git a/projects/IFS/kill_zones.F90 b/projects/IFS/kill_zones.F90 index 25193ea..ee87a9c 100644 --- a/projects/IFS/kill_zones.F90 +++ b/projects/IFS/kill_zones.F90 @@ -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 @@ -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 diff --git a/projects/NEMO/kill_zones.F90 b/projects/NEMO/kill_zones.F90 index 25193ea..ee87a9c 100644 --- a/projects/NEMO/kill_zones.F90 +++ b/projects/NEMO/kill_zones.F90 @@ -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 @@ -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 diff --git a/projects/ROMS/kill_zones.F90 b/projects/ROMS/kill_zones.F90 index 25193ea..ee87a9c 100644 --- a/projects/ROMS/kill_zones.F90 +++ b/projects/ROMS/kill_zones.F90 @@ -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 @@ -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 diff --git a/projects/Theoretical/kill_zones.F90 b/projects/Theoretical/kill_zones.F90 index 25193ea..ee87a9c 100644 --- a/projects/Theoretical/kill_zones.F90 +++ b/projects/Theoretical/kill_zones.F90 @@ -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 @@ -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 diff --git a/src/mod_diffusion.F90 b/src/mod_diffusion.F90 index 4cf73de..5816a8b 100644 --- a/src/mod_diffusion.F90 +++ b/src/mod_diffusion.F90 @@ -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 diff --git a/src/mod_init.F90 b/src/mod_init.F90 index dfe76e3..82d122b 100644 --- a/src/mod_init.F90 +++ b/src/mod_init.F90 @@ -73,6 +73,9 @@ SUBROUTINE init_namelist() 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') @@ -91,6 +94,19 @@ SUBROUTINE init_namelist() READ (8,nml=INIT_ACTIVE) CLOSE(8) + ! 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 @@ -174,7 +190,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 diff --git a/src/mod_stream.F90 b/src/mod_stream.F90 index 8db9425..e545795 100644 --- a/src/mod_stream.F90 +++ b/src/mod_stream.F90 @@ -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 @@ -170,7 +170,7 @@ SUBROUTINE init_stream() INTEGER :: index - index = 21 + index = maxlbas IF (l_offline .EQV. .FALSE.) THEN index = ntracmax END IF diff --git a/src/mod_subdomain.F90 b/src/mod_subdomain.F90 index 7501bf0..7530cab 100644 --- a/src/mod_subdomain.F90 +++ b/src/mod_subdomain.F90 @@ -14,6 +14,8 @@ MODULE mod_subdomain USE mod_domain USE mod_grid + USE mod_tracervars + USE mod_postprocessvars IMPLICIT NONE @@ -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 @@ -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 @@ -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; @@ -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) diff --git a/src/mod_tracers.F90 b/src/mod_tracers.F90 index 19f4a13..43fa099 100644 --- a/src/mod_tracers.F90 +++ b/src/mod_tracers.F90 @@ -20,7 +20,7 @@ MODULE mod_tracers IMPLICIT NONE - INTEGER, DIMENSION(10) :: numtracerarray = 0 + INTEGER, DIMENSION(MAXTRACERS) :: numtracerarray = 0 PRIVATE :: tracers_default @@ -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)) diff --git a/src/mod_vars.F90 b/src/mod_vars.F90 index 6fcc26a..5d40b1d 100755 --- a/src/mod_vars.F90 +++ b/src/mod_vars.F90 @@ -31,6 +31,15 @@ MODULE mod_precdef INTEGER, PARAMETER :: PP = selected_real_kind(6 , 37) INTEGER, PARAMETER :: DP = selected_real_kind(15, 307) INTEGER, PARAMETER :: QP = selected_real_kind(33, 4931) + + ! Maximum capacities for killing zones and tracers. These statically + ! dimension the configuration arrays; runtime loops should derive their + ! bounds from SIZE() of the array, and the actual number of streamfunction + ! zones used is the computed variable maxlbas. + INTEGER, PARAMETER :: MAXGEOZONES = 10 ! geographic killing-zone slots (ienw/iene/jens/jenn) + INTEGER, PARAMETER :: MAXTRACERS = 10 ! tracer slots (tracerchoice, etc.) + ! Max killing zones (lbas): 1 (surface) + tracer killing zones + geographic. + INTEGER, PARAMETER :: MAXZONES = 1 + MAXTRACERS + MAXGEOZONES ENDMODULE mod_precdef ! Verbose options @@ -321,7 +330,8 @@ MODULE mod_domain LOGICAL :: l_nosurface = .FALSE. ! Can trajectories reach the surface? INTEGER :: exitType - INTEGER, DIMENSION(10) :: ienw ,iene, jens ,jenn ! Horizontal killing zones + INTEGER, DIMENSION(MAXGEOZONES) :: ienw = 0, iene = 0, jens = 0, jenn = 0 ! Horizontal killing zones + INTEGER :: ngeozones = 0 ! Highest occupied geographic killing-zone slot REAL(PP) :: timax ! Maximum time to run a trajectory @@ -374,15 +384,15 @@ MODULE mod_tracervars INTEGER :: numtracers = 0 ! Tracer choice - INTEGER, DIMENSION(10) :: tracerchoice = 999, maxormin = 1 + INTEGER, DIMENSION(MAXTRACERS) :: tracerchoice = 999, maxormin = 1 ! Tracer characteristics - CHARACTER(len=100), DIMENSION(10) :: tracername = '', tracerunit, & + CHARACTER(len=100), DIMENSION(MAXTRACERS) :: tracername = '', tracerunit, & tracervarname,traceraction - CHARACTER(len=2), DIMENSION(10) :: tracerdimension = '3D' + CHARACTER(len=2), DIMENSION(MAXTRACERS) :: tracerdimension = '3D' - REAL(DP), DIMENSION(10) :: tracermin, tracermax, & + REAL(DP), DIMENSION(MAXTRACERS) :: tracermin, tracermax, & tracer0min=-9999.d0, tracer0max=9999.d0, & tracershift=0.d0, tracerscale=1.d0, & tracere @@ -418,7 +428,7 @@ MODULE mod_psi INTEGER :: xyflux = 1 ! Direction integration - INTEGER , DIMENSION(21) :: dirpsi = 1 + INTEGER , DIMENSION(MAXZONES) :: dirpsi = 1 ! Barotropic streamfunction REAL(DP), ALLOCATABLE, DIMENSION(:,:,:) :: fluxes_xy @@ -457,11 +467,11 @@ MODULE mod_postprocessvars INTEGER, DIMENSION(:),ALLOCATABLE :: nsavewrite INTEGER :: nsave = 0 - INTEGER, DIMENSION(0:21) :: ntrajout = 0 + INTEGER, DIMENSION(0:MAXZONES) :: ntrajout = 0 INTEGER :: ntrajtot = 0 - INTEGER :: maxlbas = 0 + INTEGER :: maxlbas = MAXZONES - REAL(DP), DIMENSION(0:21) :: volout = 0 + REAL(DP), DIMENSION(0:MAXZONES) :: volout = 0 REAL(DP) :: voltot = 0 ! Temporary trajectory information @@ -495,6 +505,6 @@ MODULE mod_divvars LOGICAL :: l_divergence = .FALSE. REAL(DP), DIMENSION(:,:,:,:), ALLOCATABLE :: tracerdiv - REAL(DP), DIMENSION(10) :: divconst = 1.d0 + REAL(DP), DIMENSION(MAXTRACERS) :: divconst = 1.d0 END MODULE mod_divvars