From 260077daa4014709e26bcd5e378c03c0cf9ac4bd Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Thu, 15 Jan 2026 16:41:02 +0000 Subject: [PATCH 1/9] Add support for time-varying masked data --- Utilities/mean_nemo/mean_nemo.f90 | 384 ++++++++++++++++++++---------- 1 file changed, 260 insertions(+), 124 deletions(-) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index baaa757..1950b61 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -37,7 +37,7 @@ PROGRAM mean_nemo INTEGER :: nargs, ifile , iargc, no_fill INTEGER :: ncid, outid, iostat, idim, istop, itime INTEGER :: natts, attid, xtype, varid, icellthick_type - ! ntimes is total number of time points to average over + ! ntimes is total number of time points to average over (for unmasked data) ! ntimes_local is number of time points in each file INTEGER :: jv_loop, jv, jv_thickness, ndims, nvars, dimlen, dimids(4), ntimes, ntimes_local INTEGER :: dimid, unlimitedDimId, unlimitedDimId_local, varunlimitedDimId @@ -46,15 +46,13 @@ PROGRAM mean_nemo INTEGER, ALLOCATABLE :: indimlens(:), start(:) ! Scalars - INTEGER(i1) :: ntimes_i1 - INTEGER(i2) :: ntimes_i2 - INTEGER(i4) :: ntimes_i4, inputdata_fill_value_i4 - REAL(sp) :: ntimes_sp, inputdata_fill_value_sp - REAL(dp) :: ntimes_dp, inputdata_fill_value_dp + INTEGER(i4) :: inputdata_fill_value_i4 + REAL(sp) :: inputdata_fill_value_sp + REAL(dp) :: inputdata_fill_value_dp - ! Logical data masks - LOGICAL, ALLOCATABLE, DIMENSION(:,:,: ) :: l_mask_3d - LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: l_mask_4d + ! Time counters (for masked data) + INTEGER(i2), ALLOCATABLE, DIMENSION(:,:,:) :: ntimes_3d + INTEGER(i2), ALLOCATABLE, DIMENSION(:,:,:,:) :: ntimes_4d !Int 1 versions of the local data arrays INTEGER(i1), ALLOCATABLE, SAVE, DIMENSION(:) :: inputdata_1d_i1 @@ -127,7 +125,7 @@ PROGRAM mean_nemo - !End of definitions + !End of definitions !-------------------------------------------------------------------------------- !1. Read in the arguments (input and output filenames) @@ -152,7 +150,7 @@ PROGRAM mean_nemo !--------------------------------------------------------------------------- !2. Read in the global dimensions from the first input file and set up the output file - + iostat = nf90_open( TRIM(filenames(1)), nf90_share, ncid ) IF( iostat /= nf90_noerr ) THEN WRITE(6,*) TRIM(nf90_strerror(iostat)) @@ -198,7 +196,7 @@ PROGRAM mean_nemo zone iostat = nf90_put_att( outid, nf90_global, "TimeStamp", timestamp) IF (l_verbose) WRITE(6,*)'Writing new TimeStamp attribute' - + !2.2.3 Copy the variable definitions and attributes into the output file. DO jv = 1, nvars iostat = nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts) @@ -222,7 +220,7 @@ PROGRAM mean_nemo iostat = nf90_enddef( outid ) inncids(1) = ncid IF (l_verbose) WRITE(6,*)'Finished defining output file.' - + !--------------------------------------------------------------------------- !3. Read in data from each file for each variable @@ -252,10 +250,10 @@ PROGRAM mean_nemo EXIT ENDIF ENDDO - + !Loop over all variables in first input file DO jv_loop = 0, nvars - + !Need to make sure that we mean up any cell thickness variable first so we can subsequently !use the time-mean cell thickness in the meaning of thickness-weighted variables. IF( jv_loop == 0 ) THEN @@ -272,9 +270,6 @@ PROGRAM mean_nemo ENDIF ENDIF - !Initialise ntimes (number of records to be averaged) - ntimes = 0 - !3.2 Inquire variable to find out name and how many dimensions it has ncid = inncids(1) @@ -329,10 +324,13 @@ PROGRAM mean_nemo IF( l_thckwgt ) WRITE(6,*)'Applying thickness-weighting.' ENDIF - ! Allocate global variables ahead of looping over input files + ! Allocate and initialise arrays used in the calculation of the average IF( ndims == 1 ) THEN + ! Denominator- # of records, always a scalar (no support for masked averages) + ntimes = 0 + ! Numerator- sum over time of the data SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(meandata_1d_i1(outdimlens(dimids(1)))) @@ -355,7 +353,10 @@ PROGRAM mean_nemo END SELECT ELSEIF( ndims == 2 ) THEN + ! Denominator- # of records, always a scalar (no support for masked averages) + ntimes = 0 + ! Numerator- sum over time of the data SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(meandata_2d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)))) @@ -378,7 +379,16 @@ PROGRAM mean_nemo END SELECT ELSEIF( ndims == 3 ) THEN + ! Denominator- # of records, either a scalar (normal average) or array (masked average) + IF( l_ismasked ) THEN + ALLOCATE( ntimes_3d( outdimlens(dimids(1)), outdimlens(dimids(2)), & + & outdimlens(dimids(3)) ) ) + ntimes_3d(:,:,:) = 0 + ELSE + ntimes = 0 + ENDIF + ! Numerator- sum over time of the data SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(meandata_3d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)), & @@ -404,11 +414,19 @@ PROGRAM mean_nemo WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT - IF( l_ismasked ) ALLOCATE(l_mask_3d(outdimlens(dimids(1)), outdimlens(dimids(2)), & - & outdimlens(dimids(3)))) ELSEIF( ndims == 4 ) THEN + IF( l_ismasked ) THEN + ! Initialise ntimes (number of records to be averaged) + ALLOCATE( ntimes_4d( outdimlens(dimids(1)), outdimlens(dimids(2)), & + & outdimlens(dimids(3)), outdimlens(dimids(4)) ) ) + ntimes_4d(:,:,:,:) = 0 + ELSE + ntimes = 0 + ENDIF + + ! Numerator- sum over time of the data (potentially thickness-weighted) SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(meandata_4d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)), & @@ -446,8 +464,6 @@ PROGRAM mean_nemo WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT - IF( l_ismasked ) ALLOCATE(l_mask_4d(outdimlens(dimids(1)), outdimlens(dimids(2)), & - & outdimlens(dimids(3)), outdimlens(dimids(4)))) ELSE WRITE(6,*)'E R R O R: ' WRITE(6,*)'The netcdf variable has more than 4 dimensions which is not taken into account' @@ -478,15 +494,16 @@ PROGRAM mean_nemo indimlens(idim)=dimlen ENDIF END DO - ntimes = ntimes + ntimes_local - - DO itime = 1, ntimes_local !Loop through records in file + + ! Loop through records in file, accumulate the numerator & denominator of the average + DO itime = 1, ntimes_local !start is the offset variable used in call to nf90_get_var start(varunlimitedDimId)=itime IF( ndims == 1 ) THEN - + ntimes = ntimes + 1 + SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(inputdata_1d_i1(outdimlens(dimids(1)))) @@ -524,6 +541,7 @@ PROGRAM mean_nemo END SELECT ELSEIF( ndims == 2 ) THEN + ntimes = ntimes + 1 SELECT CASE( xtype ) CASE( NF90_BYTE ) @@ -571,6 +589,8 @@ PROGRAM mean_nemo iostat = nf90_get_var( ncid, jv, inputdata_3d_i1, start, indimlens ) meandata_3d_i1(:,:,:)=meandata_3d_i1(:,:,:)+inputdata_3d_i1(:,:,:) DEALLOCATE(inputdata_3d_i1) + + ntimes = ntimes + 1 CASE( NF90_SHORT ) ALLOCATE(inputdata_3d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) @@ -578,38 +598,61 @@ PROGRAM mean_nemo iostat = nf90_get_var( ncid, jv, inputdata_3d_i2, start, indimlens ) meandata_3d_i2(:,:,:)=meandata_3d_i2(:,:,:)+inputdata_3d_i2(:,:,:) DEALLOCATE(inputdata_3d_i2) + + ntimes = ntimes + 1 CASE( NF90_INT ) ALLOCATE(inputdata_3d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) inputdata_3d_i4(:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_3d_i4, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_3d(:,:,:) = (inputdata_3d_i4(:,:,:) == inputdata_fill_value_i4) + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + WHERE( inputdata_3d_i4(:,:,:) /= inputdata_fill_value_i4 ) + meandata_3d_i4(:,:,:) = meandata_3d_i4(:,:,:) + inputdata_3d_i4(:,:,:) + ntimes_3d(:,:,:) = ntimes_3d(:,:,:) + 1 + ENDWHERE + ELSE + meandata_3d_i4(:,:,:) = meandata_3d_i4(:,:,:) + inputdata_3d_i4(:,:,:) + ntimes = ntimes + 1 ENDIF - meandata_3d_i4(:,:,:)=meandata_3d_i4(:,:,:)+inputdata_3d_i4(:,:,:) + DEALLOCATE(inputdata_3d_i4) CASE( NF90_FLOAT ) ALLOCATE(inputdata_3d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) inputdata_3d_sp(:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_3d_sp, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_3d(:,:,:) = (inputdata_3d_sp(:,:,:) == inputdata_fill_value_sp) + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + WHERE( inputdata_3d_sp(:,:,:) /= inputdata_fill_value_sp ) + meandata_3d_sp(:,:,:) = meandata_3d_sp(:,:,:) + inputdata_3d_sp(:,:,:) + ntimes_3d(:,:,:) = ntimes_3d(:,:,:) + 1 + ENDWHERE + ELSE + meandata_3d_sp(:,:,:) = meandata_3d_sp(:,:,:) + inputdata_3d_sp(:,:,:) + ntimes = ntimes + 1 ENDIF - meandata_3d_sp(:,:,:)=meandata_3d_sp(:,:,:)+inputdata_3d_sp(:,:,:) + DEALLOCATE(inputdata_3d_sp) CASE( NF90_DOUBLE ) ALLOCATE(inputdata_3d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) inputdata_3d_dp(:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_3d_dp, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_3d(:,:,:) = (inputdata_3d_dp(:,:,:) == inputdata_fill_value_dp) + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + WHERE( inputdata_3d_dp(:,:,:) /= inputdata_fill_value_dp ) + meandata_3d_dp(:,:,:) = meandata_3d_dp(:,:,:) + inputdata_3d_dp(:,:,:) + ntimes_3d(:,:,:) = ntimes_3d(:,:,:) + 1 + ENDWHERE + ELSE + meandata_3d_dp(:,:,:) = meandata_3d_dp(:,:,:) + inputdata_3d_dp(:,:,:) + ntimes = ntimes + 1 ENDIF - meandata_3d_dp(:,:,:)=meandata_3d_dp(:,:,:)+inputdata_3d_dp(:,:,:) + DEALLOCATE(inputdata_3d_dp) CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype @@ -626,6 +669,8 @@ PROGRAM mean_nemo iostat = nf90_get_var( ncid, jv, inputdata_4d_i1, start, indimlens ) meandata_4d_i1(:,:,:,:)=meandata_4d_i1(:,:,:,:)+inputdata_4d_i1(:,:,:,:) DEALLOCATE(inputdata_4d_i1) + + ntimes = ntimes + 1 CASE( NF90_SHORT ) ALLOCATE(inputdata_4d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) @@ -633,62 +678,101 @@ PROGRAM mean_nemo iostat = nf90_get_var( ncid, jv, inputdata_4d_i2, start, indimlens ) meandata_4d_i2(:,:,:,:)=meandata_4d_i2(:,:,:,:)+inputdata_4d_i2(:,:,:,:) DEALLOCATE(inputdata_4d_i2) + + ntimes = ntimes + 1 CASE( NF90_INT ) ALLOCATE(inputdata_4d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) inputdata_4d_i4(:,:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_4d_i4, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_4d(:,:,:,:) = (inputdata_4d_i4(:,:,:,:) == inputdata_fill_value_i4) + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + WHERE( inputdata_4d_i4(:,:,:,:) /= inputdata_fill_value_i4 ) + meandata_4d_i4(:,:,:,:) = meandata_4d_i4(:,:,:,:) + inputdata_4d_i4(:,:,:,:) + ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + ENDWHERE + ELSE + meandata_4d_i4(:,:,:,:) = meandata_4d_i4(:,:,:,:) + inputdata_4d_i4(:,:,:,:) + ntimes = ntimes + 1 ENDIF - meandata_4d_i4(:,:,:,:)=meandata_4d_i4(:,:,:,:)+inputdata_4d_i4(:,:,:,:) + DEALLOCATE(inputdata_4d_i4) CASE( NF90_FLOAT ) ALLOCATE(inputdata_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) inputdata_4d_sp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_4d_sp, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_4d(:,:,:,:) = (inputdata_4d_sp(:,:,:,:) == inputdata_fill_value_sp) - ENDIF - ! Thickness-weighting + + ! Thickness-weighting- get cell thickness IF( l_thckwgt ) THEN ALLOCATE(cellthick_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) cellthick_4d_sp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv_thickness, cellthick_4d_sp, start, indimlens ) - ! Avoid floating point errors - IF( l_ismasked ) THEN ; WHERE( l_mask_4d ) cellthick_4d_sp(:,:,:,:) = 1. ; ENDIF - meandata_4d_sp(:,:,:,:)=meandata_4d_sp(:,:,:,:)+inputdata_4d_sp(:,:,:,:)*cellthick_4d_sp(:,:,:,:) - DEALLOCATE(cellthick_4d_sp) + ENDIF + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + IF( l_thckwgt ) THEN + WHERE( inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp ) + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) * cellthick_4d_sp(:,:,:,:) + ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + ENDWHERE + ELSE + WHERE( inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp ) + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) + ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + ENDWHERE + ENDIF ELSE - meandata_4d_sp(:,:,:,:)=meandata_4d_sp(:,:,:,:)+inputdata_4d_sp(:,:,:,:) + IF( l_thckwgt ) THEN + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) * cellthick_4d_sp(:,:,:,:) + ELSE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) + ENDIF + ntimes = ntimes + 1 ENDIF + + IF( l_thckwgt ) DEALLOCATE(cellthick_4d_sp) DEALLOCATE(inputdata_4d_sp) CASE( NF90_DOUBLE ) ALLOCATE(inputdata_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) inputdata_4d_dp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_4d_dp, start, indimlens ) - ! If the variable is masked, get the mask for the first time and file (assume constant in time) - IF( l_ismasked .AND. ifile == 1 .AND. itime == 1 ) THEN - l_mask_4d(:,:,:,:) = (inputdata_4d_dp(:,:,:,:) == inputdata_fill_value_dp) - ENDIF - ! Thickness-weighting + + ! Thickness-weighting- get cell thickness IF( l_thckwgt ) THEN ALLOCATE(cellthick_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) cellthick_4d_dp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv_thickness, cellthick_4d_dp, start, indimlens ) - ! Avoid floating point errors - IF( l_ismasked ) THEN ; WHERE( l_mask_4d ) cellthick_4d_dp(:,:,:,:) = 1. ; ENDIF - meandata_4d_dp(:,:,:,:)=meandata_4d_dp(:,:,:,:)+inputdata_4d_dp(:,:,:,:)*cellthick_4d_dp(:,:,:,:) - DEALLOCATE(cellthick_4d_dp) + ENDIF + + ! Do not include masked data in the average + IF( l_ismasked ) THEN + IF( l_thckwgt ) THEN + WHERE( inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp ) + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) * cellthick_4d_dp(:,:,:,:) + ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + ENDWHERE + ELSE + WHERE( inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp ) + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) + ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + ENDWHERE + ENDIF ELSE - meandata_4d_dp(:,:,:,:)=meandata_4d_dp(:,:,:,:)+inputdata_4d_dp(:,:,:,:) + IF( l_thckwgt ) THEN + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) * cellthick_4d_dp(:,:,:,:) + ELSE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) + ENDIF + ntimes = ntimes + 1 ENDIF + + IF( l_thckwgt ) DEALLOCATE(cellthick_4d_dp) DEALLOCATE(inputdata_4d_dp) CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype @@ -716,26 +800,19 @@ PROGRAM mean_nemo END DO !loop over files !Divide by number of records to get mean - ! cast ntimes to appropriate types to allow for the divisions - ntimes_i1 = INT(ntimes) - ntimes_i2 = INT(ntimes) - ntimes_i4 = INT(ntimes) - ntimes_sp = REAL(ntimes) - ntimes_dp = REAL(ntimes) - IF( ndims == 1 ) THEN SELECT CASE( xtype ) CASE( NF90_BYTE ) - meandata_1d_i1(:)=meandata_1d_i1(:)/(ntimes_i1) + meandata_1d_i1(:)=meandata_1d_i1(:)/ntimes CASE( NF90_SHORT ) - meandata_1d_i2(:)=meandata_1d_i2(:)/(ntimes_i2) + meandata_1d_i2(:)=meandata_1d_i2(:)/ntimes CASE( NF90_INT ) - meandata_1d_i4(:)=meandata_1d_i4(:)/(ntimes_i4) + meandata_1d_i4(:)=meandata_1d_i4(:)/ntimes CASE( NF90_FLOAT ) - meandata_1d_sp(:)=meandata_1d_sp(:)/(ntimes_sp) + meandata_1d_sp(:)=meandata_1d_sp(:)/ntimes CASE( NF90_DOUBLE ) - meandata_1d_dp(:)=meandata_1d_dp(:)/(ntimes_dp) + meandata_1d_dp(:)=meandata_1d_dp(:)/ntimes CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 @@ -745,15 +822,15 @@ PROGRAM mean_nemo SELECT CASE( xtype ) CASE( NF90_BYTE ) - meandata_2d_i1(:,:)=meandata_2d_i1(:,:)/(ntimes_i1) + meandata_2d_i1(:,:)=meandata_2d_i1(:,:)/ntimes CASE( NF90_SHORT ) - meandata_2d_i2(:,:)=meandata_2d_i2(:,:)/(ntimes_i2) + meandata_2d_i2(:,:)=meandata_2d_i2(:,:)/ntimes CASE( NF90_INT ) - meandata_2d_i4(:,:)=meandata_2d_i4(:,:)/(ntimes_i4) + meandata_2d_i4(:,:)=meandata_2d_i4(:,:)/ntimes CASE( NF90_FLOAT ) - meandata_2d_sp(:,:)=meandata_2d_sp(:,:)/(ntimes_sp) + meandata_2d_sp(:,:)=meandata_2d_sp(:,:)/ntimes CASE( NF90_DOUBLE ) - meandata_2d_dp(:,:)=meandata_2d_dp(:,:)/(ntimes_dp) + meandata_2d_dp(:,:)=meandata_2d_dp(:,:)/ntimes CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 @@ -763,81 +840,142 @@ PROGRAM mean_nemo SELECT CASE( xtype ) CASE( NF90_BYTE ) - meandata_3d_i1(:,:,:)=meandata_3d_i1(:,:,:)/(ntimes_i1) + meandata_3d_i1(:,:,:)=meandata_3d_i1(:,:,:)/ntimes CASE( NF90_SHORT ) - meandata_3d_i2(:,:,:)=meandata_3d_i2(:,:,:)/(ntimes_i2) + meandata_3d_i2(:,:,:)=meandata_3d_i2(:,:,:)/ntimes CASE( NF90_INT ) - meandata_3d_i4(:,:,:)=meandata_3d_i4(:,:,:)/(ntimes_i4) - ! Restore masked points IF( l_ismasked ) THEN - WHERE( l_mask_3d ) meandata_3d_i4(:,:,:) = inputdata_fill_value_i4 + WHERE( ntimes_3d(:,:,:) == 0 ) + meandata_3d_i4(:,:,:) = inputdata_fill_value_i4 + ELSEWHERE + meandata_3d_i4(:,:,:) = meandata_3d_i4(:,:,:) / ntimes_3d(:,:,:) + ENDWHERE + ELSE + meandata_3d_i4(:,:,:) = meandata_3d_i4(:,:,:) / ntimes ENDIF CASE( NF90_FLOAT ) - meandata_3d_sp(:,:,:)=meandata_3d_sp(:,:,:)/(ntimes_sp) - ! Restore masked points IF( l_ismasked ) THEN - WHERE( l_mask_3d ) meandata_3d_sp(:,:,:) = inputdata_fill_value_sp + WHERE( ntimes_3d(:,:,:) == 0 ) + meandata_3d_sp(:,:,:) = inputdata_fill_value_sp + ELSEWHERE + meandata_3d_sp(:,:,:) = meandata_3d_sp(:,:,:) / ntimes_3d(:,:,:) + ENDWHERE + ELSE + meandata_3d_sp(:,:,:) = meandata_3d_sp(:,:,:) / ntimes ENDIF CASE( NF90_DOUBLE ) - meandata_3d_dp(:,:,:)=meandata_3d_dp(:,:,:)/(ntimes_dp) - ! Restore masked points IF( l_ismasked ) THEN - WHERE( l_mask_3d ) meandata_3d_dp(:,:,:) = inputdata_fill_value_dp + WHERE( ntimes_3d(:,:,:) == 0 ) + meandata_3d_dp(:,:,:) = inputdata_fill_value_dp + ELSEWHERE + meandata_3d_dp(:,:,:) = meandata_3d_dp(:,:,:) / ntimes_3d(:,:,:) + ENDWHERE + ELSE + meandata_3d_dp(:,:,:) = meandata_3d_dp(:,:,:) / ntimes ENDIF CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT + IF( l_ismasked ) DEALLOCATE(ntimes_3d) ELSEIF( ndims == 4 ) THEN SELECT CASE( xtype ) CASE( NF90_BYTE ) - meandata_4d_i1(:,:,:,:)=meandata_4d_i1(:,:,:,:)/(ntimes_i1) + meandata_4d_i1(:,:,:,:)=meandata_4d_i1(:,:,:,:)/ntimes CASE( NF90_SHORT ) - meandata_4d_i2(:,:,:,:)=meandata_4d_i2(:,:,:,:)/(ntimes_i2) + meandata_4d_i2(:,:,:,:)=meandata_4d_i2(:,:,:,:)/ntimes CASE( NF90_INT ) - meandata_4d_i4(:,:,:,:)=meandata_4d_i4(:,:,:,:)/(ntimes_i4) - ! Restore masked points IF( l_ismasked ) THEN - WHERE( l_mask_4d ) meandata_4d_i4(:,:,:,:) = inputdata_fill_value_i4 + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_i4(:,:,:,:) = inputdata_fill_value_i4 + ELSEWHERE + meandata_4d_i4(:,:,:,:) = meandata_4d_i4(:,:,:,:) / ntimes_4d(:,:,:,:) + ENDWHERE + ELSE + meandata_4d_i4(:,:,:,:) = meandata_4d_i4(:,:,:,:) / ntimes ENDIF CASE( NF90_FLOAT ) - IF(l_thckwgt) THEN - IF( icellthick_type == NF90_DOUBLE ) THEN - meandata_4d_sp(:,:,:,:)=meandata_4d_sp(:,:,:,:)/( REAL(meancellthick_4d_dp, sp) * ntimes_sp ) + IF( l_ismasked ) THEN + IF( l_thckwgt ) THEN + IF( icellthick_type == NF90_DOUBLE ) THEN + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp + ELSEWHERE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (REAL(meancellthick_4d_dp, sp) * ntimes_4d(:,:,:,:)) + ENDWHERE + ELSE + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp + ELSEWHERE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (meancellthick_4d_sp(:,:,:,:) * ntimes_4d(:,:,:,:)) + ENDWHERE + ENDIF ELSE - meandata_4d_sp(:,:,:,:)=meandata_4d_sp(:,:,:,:)/(meancellthick_4d_sp(:,:,:,:) * ntimes_sp) + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp + ELSEWHERE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / ntimes_4d(:,:,:,:) + ENDWHERE ENDIF ELSE - meandata_4d_sp(:,:,:,:)=meandata_4d_sp(:,:,:,:)/(ntimes_sp) - ENDIF - ! Restore masked points - IF( l_ismasked ) THEN - WHERE( l_mask_4d ) meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp + IF( l_thckwgt ) THEN + IF( icellthick_type == NF90_DOUBLE ) THEN + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (REAL(meancellthick_4d_dp, sp) * ntimes) + ELSE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (meancellthick_4d_sp(:,:,:,:) * ntimes) + ENDIF + ELSE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / ntimes + ENDIF ENDIF + ! If this is the cell thickness, save to normalise thickness-weighted time-means - IF( jv == jv_thickness ) meancellthick_4d_sp = meandata_4d_sp + IF( jv == jv_thickness ) meancellthick_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) CASE( NF90_DOUBLE ) - IF(l_thckwgt) THEN - IF( icellthick_type == NF90_FLOAT ) THEN - meandata_4d_dp(:,:,:,:)=meandata_4d_dp(:,:,:,:)/( REAL(meancellthick_4d_sp, dp) * ntimes_dp ) + IF( l_ismasked ) THEN + IF( l_thckwgt ) THEN + IF( icellthick_type == NF90_FLOAT ) THEN + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp + ELSEWHERE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (REAL(meancellthick_4d_sp, dp) * ntimes_4d(:,:,:,:)) + ENDWHERE + ELSE + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp + ELSEWHERE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (meancellthick_4d_dp(:,:,:,:) * ntimes_4d(:,:,:,:)) + ENDWHERE + ENDIF ELSE - meandata_4d_dp(:,:,:,:)=meandata_4d_dp(:,:,:,:)/(meancellthick_4d_dp(:,:,:,:) * ntimes_dp) + WHERE( ntimes_4d(:,:,:,:) == 0 ) + meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp + ELSEWHERE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / ntimes_4d(:,:,:,:) + ENDWHERE ENDIF ELSE - meandata_4d_dp(:,:,:,:)=meandata_4d_dp(:,:,:,:)/(ntimes_dp) - ENDIF - ! Restore masked points - IF( l_ismasked ) THEN - WHERE( l_mask_4d ) meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp + IF( l_thckwgt ) THEN + IF( icellthick_type == NF90_FLOAT ) THEN + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (REAL(meancellthick_4d_sp, dp) * ntimes) + ELSE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (meancellthick_4d_dp(:,:,:,:) * ntimes) + ENDIF + ELSE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / ntimes + ENDIF ENDIF + ! If this is the cell thickness, save to normalise thickness-weighted time-means - IF( jv == jv_thickness ) meancellthick_4d_dp = meandata_4d_dp + IF( jv == jv_thickness ) meancellthick_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT + IF( l_ismasked ) DEALLOCATE(ntimes_4d) + ENDIF ELSE @@ -957,7 +1095,7 @@ PROGRAM mean_nemo DEALLOCATE(meandata_1d_dp) END SELECT - ELSEIF( ndims == 2 ) THEN + ELSEIF( ndims == 2 ) THEN SELECT CASE( xtype ) CASE( NF90_BYTE ) @@ -979,7 +1117,7 @@ PROGRAM mean_nemo WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT - + ELSEIF( ndims == 3 ) THEN SELECT CASE( xtype ) @@ -1001,8 +1139,7 @@ PROGRAM mean_nemo CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 - END SELECT - IF( l_ismasked ) DEALLOCATE( l_mask_3d ) + END SELECT ELSEIF( ndims == 4 ) THEN @@ -1025,8 +1162,7 @@ PROGRAM mean_nemo CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 - END SELECT - IF( l_ismasked ) DEALLOCATE( l_mask_4d ) + END SELECT ENDIF @@ -1034,7 +1170,7 @@ PROGRAM mean_nemo WRITE(6,*) 'E R R O R writing variable '//TRIM(varname) STOP 17 ENDIF - + END DO !loop over variables IF( allocated(meancellthick_4d_sp) ) DEALLOCATE(meancellthick_4d_sp) @@ -1062,5 +1198,5 @@ PROGRAM mean_nemo WRITE(6,*)'E R R O R closing output file' STOP 19 ENDIF - + END PROGRAM mean_nemo From f307cfb692a5de6c798a8b92233e4d7e43fcf2a4 Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Mon, 19 Jan 2026 19:40:20 +0000 Subject: [PATCH 2/9] Add support for cell thickness masking not matching that of thickness-weighted data --- Utilities/mean_nemo/mean_nemo.f90 | 292 +++++++++++++++++------------- 1 file changed, 169 insertions(+), 123 deletions(-) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index 1950b61..f32f45a 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -32,23 +32,26 @@ PROGRAM mean_nemo CHARACTER(LEN=nf90_max_name), ALLOCATABLE :: filenames(:), indimnames(:) CHARACTER(LEN=256) :: standard_name,cell_methods - LOGICAL :: l_thckwgt, l_doavg, l_ismasked + LOGICAL :: l_thckwgt, l_doavg, l_ismasked, l_inputdata_ismasked, l_cellthick_ismasked INTEGER :: nargs, ifile , iargc, no_fill INTEGER :: ncid, outid, iostat, idim, istop, itime - INTEGER :: natts, attid, xtype, varid, icellthick_type + INTEGER :: natts, attid, xtype, varid ! ntimes is total number of time points to average over (for unmasked data) ! ntimes_local is number of time points in each file - INTEGER :: jv_loop, jv, jv_thickness, ndims, nvars, dimlen, dimids(4), ntimes, ntimes_local + INTEGER :: jv, jv_thickness, ndims, nvars, dimlen, dimids(4), ntimes, ntimes_local INTEGER :: dimid, unlimitedDimId, unlimitedDimId_local, varunlimitedDimId INTEGER :: chunksize = 32000000 INTEGER, ALLOCATABLE :: outdimids(:), outdimlens(:), inncids(:) INTEGER, ALLOCATABLE :: indimlens(:), start(:) ! Scalars - INTEGER(i4) :: inputdata_fill_value_i4 - REAL(sp) :: inputdata_fill_value_sp - REAL(dp) :: inputdata_fill_value_dp + INTEGER(i4) :: inputdata_fill_value_i4, outputdata_fill_value_i4 + REAL(sp) :: inputdata_fill_value_sp, outputdata_fill_value_sp, cellthick_fill_value_sp + REAL(dp) :: inputdata_fill_value_dp, outputdata_fill_value_dp, cellthick_fill_value_dp + + ! Logical data masks (for 4d data only, which may be thickness-weighted and/or masked) + LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: l_mask_4d ! Time counters (for masked data) INTEGER(i2), ALLOCATABLE, DIMENSION(:,:,:) :: ntimes_3d @@ -197,6 +200,38 @@ PROGRAM mean_nemo iostat = nf90_put_att( outid, nf90_global, "TimeStamp", timestamp) IF (l_verbose) WRITE(6,*)'Writing new TimeStamp attribute' + jv_thickness = -1 + l_cellthick_ismasked = .FALSE. + + ! Find out if there is a cell thickness variable in this set of files in case we need to do thickness weighting + DO jv = 1, nvars + iostat = nf90_get_att(ncid, jv, "standard_name", standard_name) + + IF( (iostat == nf90_noerr) .AND. (TRIM(standard_name) == "cell_thickness") ) THEN + jv_thickness = jv + iostat = nf90_inquire_variable( ncid, jv, xtype=xtype ) + + ! Save cell thickness fill value if masked + SELECT CASE( xtype ) + CASE( NF90_FLOAT ) + iostat = nf90_get_att(ncid, jv, "_FillValue", cellthick_fill_value_sp ) + IF( iostat == nf90_noerr ) cellthick_fill_value_dp = REAL(cellthick_fill_value_sp, dp) + CASE( NF90_DOUBLE ) + iostat = nf90_get_att(ncid, jv, "_FillValue", cellthick_fill_value_dp ) + IF( iostat == nf90_noerr ) cellthick_fill_value_sp = REAL(cellthick_fill_value_dp, sp) + CASE DEFAULT + WRITE(6,*) "ERROR : Cell thickness variable must be single or double precision float" + STOP 13 + END SELECT + + ! Is the cell thickness data masked? + l_cellthick_ismasked = (iostat == nf90_noerr) + + ! Exit loop if the cell thickness has been found + EXIT + ENDIF + END DO + !2.2.3 Copy the variable definitions and attributes into the output file. DO jv = 1, nvars iostat = nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts) @@ -204,6 +239,7 @@ PROGRAM mean_nemo DO idim = 1, ndims outdimids(idim) = dimids(idim) END DO + iostat = nf90_def_var( outid, varname, xtype, outdimids, varid ) DEALLOCATE(outdimids) IF (l_verbose) WRITE(6,*)'Defining variable '//TRIM(varname)//'...' @@ -213,6 +249,23 @@ PROGRAM mean_nemo iostat = nf90_copy_att( ncid, varid, attname, outid, varid ) END DO ENDIF + + ! For thickness weighting, output data will use the fill value of the cell thickness if the input data is not masked + IF( jv /= jv_thickness ) THEN + iostat = nf90_get_att(ncid, jv, "cell_methods", cell_methods) + l_thckwgt = ( (iostat == nf90_noerr) .AND. (TRIM(cell_methods) == "time: mean (thickness weighted)") ) + iostat = nf90_inquire_attribute(ncid, jv, "_FillValue") + l_inputdata_ismasked = (iostat == nf90_noerr .AND. ndims >= 3) + + IF( l_thckwgt .AND. l_cellthick_ismasked .AND. (.NOT. l_inputdata_ismasked) ) THEN + SELECT CASE( xtype ) + CASE( NF90_FLOAT ) + iostat = nf90_put_att( outid, varid, '_FillValue', cellthick_fill_value_sp ) + CASE( NF90_DOUBLE ) + iostat = nf90_put_att( outid, varid, '_FillValue', cellthick_fill_value_dp ) + END SELECT + ENDIF + ENDIF END DO !2.3 End definitions in output file and copy 1st file ncid to the inncids array @@ -239,36 +292,8 @@ PROGRAM mean_nemo END DO IF (l_verbose) WRITE(6,*)'All input files open.' - !Find out if there is a cell thickness variable in this set of files - !in case we need to do thickness weighting. - jv_thickness = -1 - DO jv = 1, nvars - iostat = nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts) - iostat = nf90_get_att(ncid, jv, "standard_name", standard_name) - IF( iostat == nf90_noerr .AND. TRIM(standard_name) == "cell_thickness" ) THEN - jv_thickness = jv - EXIT - ENDIF - ENDDO - !Loop over all variables in first input file - DO jv_loop = 0, nvars - - !Need to make sure that we mean up any cell thickness variable first so we can subsequently - !use the time-mean cell thickness in the meaning of thickness-weighted variables. - IF( jv_loop == 0 ) THEN - IF( jv_thickness /= -1 ) THEN - jv = jv_thickness - ELSE - CYCLE - ENDIF - ELSE - IF( jv_loop == jv_thickness ) THEN - CYCLE - ELSE - jv = jv_loop - ENDIF - ENDIF + DO jv = 1, nvars !3.2 Inquire variable to find out name and how many dimensions it has @@ -285,20 +310,42 @@ PROGRAM mean_nemo ! Does averaging need to account for thickness-weighting and masked data? l_thckwgt = .FALSE. + l_inputdata_ismasked = .FALSE. l_ismasked = .FALSE. IF( l_doavg ) THEN + ! Thickness-weighting- if an unsupported data type/shape has this attribute, raise an error below l_thckwgt = ( iostat == nf90_noerr .AND. TRIM(cell_methods) == "time: mean (thickness weighted)" ) - SELECT CASE( xtype ) - CASE( NF90_INT ) - iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_i4 ) - CASE( NF90_FLOAT ) - iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_sp ) - CASE( NF90_DOUBLE ) - iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_dp ) - CASE DEFAULT - iostat = nf90_noerr + 1 - END SELECT - l_ismasked = (iostat == nf90_noerr .AND. ndims >= 3) + + ! Masked data- get the fill value of the input data and/or cell thickness, and set the fill value of the + ! output file (that of the input data is prioritised). Unsupported data types/shapes proceed silently using + ! the unmasked averaging algorithm. + IF( ndims >= 3 ) THEN + SELECT CASE( xtype ) + CASE( NF90_INT ) + iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_i4 ) + IF( iostat == nf90_noerr ) outputdata_fill_value_i4 = inputdata_fill_value_i4 + CASE( NF90_FLOAT ) + iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_sp ) + IF( iostat == nf90_noerr ) THEN + outputdata_fill_value_sp = inputdata_fill_value_sp + ELSE IF( l_thckwgt .AND. l_cellthick_ismasked ) THEN + outputdata_fill_value_sp = cellthick_fill_value_sp + ENDIF + CASE( NF90_DOUBLE ) + iostat = nf90_get_att(ncid, jv, "_FillValue", inputdata_fill_value_dp ) + IF( iostat == nf90_noerr ) THEN + outputdata_fill_value_dp = inputdata_fill_value_dp + ELSE IF( l_thckwgt .AND. l_cellthick_ismasked ) THEN + outputdata_fill_value_dp = cellthick_fill_value_dp + ENDIF + CASE DEFAULT + iostat = nf90_noerr + 1 + END SELECT + ! Is the input data masked? + l_inputdata_ismasked = iostat == nf90_noerr + ! We need to perform masked averaging if the input data and/or cell thickness data is masked + l_ismasked = l_inputdata_ismasked .OR. (l_thckwgt .AND. l_cellthick_ismasked) + ENDIF ENDIF IF( l_thckwgt .AND. jv_thickness == -1 ) THEN @@ -417,8 +464,20 @@ PROGRAM mean_nemo ELSEIF( ndims == 4 ) THEN - IF( l_ismasked ) THEN - ! Initialise ntimes (number of records to be averaged) + ! Denominator- either the sum over time of the cell thickness (thickness-weighted average), + ! or # of records which is either a scalar (normal average) or array (masked average) + IF( l_thckwgt ) THEN + SELECT CASE( xtype ) + CASE( NF90_FLOAT ) + ALLOCATE(meancellthick_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & + & outdimlens(dimids(3)),outdimlens(dimids(4)))) + meancellthick_4d_sp(:,:,:,:)=0.0 + CASE( NF90_DOUBLE ) + ALLOCATE(meancellthick_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & + & outdimlens(dimids(3)),outdimlens(dimids(4)))) + meancellthick_4d_dp(:,:,:,:)=0.0 + END SELECT + ELSE IF( l_ismasked ) THEN ALLOCATE( ntimes_4d( outdimlens(dimids(1)), outdimlens(dimids(2)), & & outdimlens(dimids(3)), outdimlens(dimids(4)) ) ) ntimes_4d(:,:,:,:) = 0 @@ -444,22 +503,10 @@ PROGRAM mean_nemo ALLOCATE(meandata_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) meandata_4d_sp(:,:,:,:)=0.0 - IF( jv == jv_thickness ) THEN - ALLOCATE(meancellthick_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & - & outdimlens(dimids(3)),outdimlens(dimids(4)))) - meancellthick_4d_sp(:,:,:,:)=0.0 - icellthick_type = NF90_FLOAT - ENDIF CASE( NF90_DOUBLE ) ALLOCATE(meandata_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) meandata_4d_dp(:,:,:,:)=0.0 - IF( jv == jv_thickness ) THEN - ALLOCATE(meancellthick_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & - & outdimlens(dimids(3)),outdimlens(dimids(4)))) - meancellthick_4d_dp(:,:,:,:)=0.0 - icellthick_type = NF90_DOUBLE - ENDIF CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 @@ -469,7 +516,7 @@ PROGRAM mean_nemo WRITE(6,*)'The netcdf variable has more than 4 dimensions which is not taken into account' STOP 15 ENDIF - + istop = 0 ! If this variable is a function of the unlimited dimension then @@ -714,24 +761,39 @@ PROGRAM mean_nemo ! Do not include masked data in the average IF( l_ismasked ) THEN + ! Use the union of the input data and cell thickness masks + ALLOCATE(l_mask_4d(outdimlens(dimids(1)),outdimlens(dimids(2)), & + & outdimlens(dimids(3)),outdimlens(dimids(4)))) + IF( l_inputdata_ismasked .AND. (l_thckwgt .AND. l_cellthick_ismasked) ) THEN + l_mask_4d(:,:,:,:) = (inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp) .AND. & + & (cellthick_4d_sp(:,:,:,:) /= cellthick_fill_value_sp) + ELSE IF( l_inputdata_ismasked ) THEN + l_mask_4d(:,:,:,:) = (inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp) + ELSE + l_mask_4d(:,:,:,:) = (cellthick_4d_sp(:,:,:,:) /= cellthick_fill_value_sp) + ENDIF + IF( l_thckwgt ) THEN - WHERE( inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp ) + WHERE( l_mask_4d ) meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) * cellthick_4d_sp(:,:,:,:) - ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + meancellthick_4d_sp(:,:,:,:) = meancellthick_4d_sp(:,:,:,:) + cellthick_4d_sp(:,:,:,:) ENDWHERE ELSE - WHERE( inputdata_4d_sp(:,:,:,:) /= inputdata_fill_value_sp ) + WHERE( l_mask_4d ) meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 ENDWHERE ENDIF + + DEALLOCATE(l_mask_4d) ELSE IF( l_thckwgt ) THEN meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) * cellthick_4d_sp(:,:,:,:) + meancellthick_4d_sp(:,:,:,:) = meancellthick_4d_sp(:,:,:,:) + cellthick_4d_sp(:,:,:,:) ELSE meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + inputdata_4d_sp(:,:,:,:) + ntimes = ntimes + 1 ENDIF - ntimes = ntimes + 1 ENDIF IF( l_thckwgt ) DEALLOCATE(cellthick_4d_sp) @@ -752,24 +814,39 @@ PROGRAM mean_nemo ! Do not include masked data in the average IF( l_ismasked ) THEN + ! Use the union of the input data and cell thickness masks + ALLOCATE(l_mask_4d(outdimlens(dimids(1)),outdimlens(dimids(2)), & + & outdimlens(dimids(3)),outdimlens(dimids(4)))) + IF( l_inputdata_ismasked .AND. (l_thckwgt .AND. l_cellthick_ismasked) ) THEN + l_mask_4d(:,:,:,:) = (inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp) .AND. & + & (cellthick_4d_dp(:,:,:,:) /= cellthick_fill_value_dp) + ELSE IF( l_inputdata_ismasked ) THEN + l_mask_4d(:,:,:,:) = (inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp) + ELSE + l_mask_4d(:,:,:,:) = (cellthick_4d_dp(:,:,:,:) /= cellthick_fill_value_dp) + ENDIF + IF( l_thckwgt ) THEN - WHERE( inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp ) + WHERE( l_mask_4d ) meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) * cellthick_4d_dp(:,:,:,:) - ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 + meancellthick_4d_dp(:,:,:,:) = meancellthick_4d_dp(:,:,:,:) + cellthick_4d_dp(:,:,:,:) ENDWHERE ELSE - WHERE( inputdata_4d_dp(:,:,:,:) /= inputdata_fill_value_dp ) + WHERE( l_mask_4d ) meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) ntimes_4d(:,:,:,:) = ntimes_4d(:,:,:,:) + 1 ENDWHERE ENDIF + + DEALLOCATE(l_mask_4d) ELSE IF( l_thckwgt ) THEN meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) * cellthick_4d_dp(:,:,:,:) + meancellthick_4d_dp(:,:,:,:) = meancellthick_4d_dp(:,:,:,:) + cellthick_4d_dp(:,:,:,:) ELSE meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + inputdata_4d_dp(:,:,:,:) + ntimes = ntimes + 1 ENDIF - ntimes = ntimes + 1 ENDIF IF( l_thckwgt ) DEALLOCATE(cellthick_4d_dp) @@ -799,7 +876,8 @@ PROGRAM mean_nemo END DO !loop over files - !Divide by number of records to get mean + ! Divide numerator (sum of data, possibly thickness-weighted) by denominator (sum of cell thickness + ! or number of time records) to get the mean, setting masked results to a fill value IF( ndims == 1 ) THEN SELECT CASE( xtype ) @@ -846,7 +924,7 @@ PROGRAM mean_nemo CASE( NF90_INT ) IF( l_ismasked ) THEN WHERE( ntimes_3d(:,:,:) == 0 ) - meandata_3d_i4(:,:,:) = inputdata_fill_value_i4 + meandata_3d_i4(:,:,:) = outputdata_fill_value_i4 ELSEWHERE meandata_3d_i4(:,:,:) = meandata_3d_i4(:,:,:) / ntimes_3d(:,:,:) ENDWHERE @@ -856,7 +934,7 @@ PROGRAM mean_nemo CASE( NF90_FLOAT ) IF( l_ismasked ) THEN WHERE( ntimes_3d(:,:,:) == 0 ) - meandata_3d_sp(:,:,:) = inputdata_fill_value_sp + meandata_3d_sp(:,:,:) = outputdata_fill_value_sp ELSEWHERE meandata_3d_sp(:,:,:) = meandata_3d_sp(:,:,:) / ntimes_3d(:,:,:) ENDWHERE @@ -866,7 +944,7 @@ PROGRAM mean_nemo CASE( NF90_DOUBLE ) IF( l_ismasked ) THEN WHERE( ntimes_3d(:,:,:) == 0 ) - meandata_3d_dp(:,:,:) = inputdata_fill_value_dp + meandata_3d_dp(:,:,:) = outputdata_fill_value_dp ELSEWHERE meandata_3d_dp(:,:,:) = meandata_3d_dp(:,:,:) / ntimes_3d(:,:,:) ENDWHERE @@ -889,7 +967,7 @@ PROGRAM mean_nemo CASE( NF90_INT ) IF( l_ismasked ) THEN WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_i4(:,:,:,:) = inputdata_fill_value_i4 + meandata_4d_i4(:,:,:,:) = outputdata_fill_value_i4 ELSEWHERE meandata_4d_i4(:,:,:,:) = meandata_4d_i4(:,:,:,:) / ntimes_4d(:,:,:,:) ENDWHERE @@ -899,83 +977,54 @@ PROGRAM mean_nemo CASE( NF90_FLOAT ) IF( l_ismasked ) THEN IF( l_thckwgt ) THEN - IF( icellthick_type == NF90_DOUBLE ) THEN - WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp - ELSEWHERE - meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (REAL(meancellthick_4d_dp, sp) * ntimes_4d(:,:,:,:)) - ENDWHERE - ELSE - WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp - ELSEWHERE - meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (meancellthick_4d_sp(:,:,:,:) * ntimes_4d(:,:,:,:)) - ENDWHERE - ENDIF + WHERE( meancellthick_4d_sp(:,:,:,:) == 0._sp ) + meandata_4d_sp(:,:,:,:) = outputdata_fill_value_sp + ELSEWHERE + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / meancellthick_4d_sp(:,:,:,:) + ENDWHERE ELSE WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_sp(:,:,:,:) = inputdata_fill_value_sp + meandata_4d_sp(:,:,:,:) = outputdata_fill_value_sp ELSEWHERE meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / ntimes_4d(:,:,:,:) ENDWHERE ENDIF ELSE IF( l_thckwgt ) THEN - IF( icellthick_type == NF90_DOUBLE ) THEN - meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (REAL(meancellthick_4d_dp, sp) * ntimes) - ELSE - meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / (meancellthick_4d_sp(:,:,:,:) * ntimes) - ENDIF + meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / meancellthick_4d_sp(:,:,:,:) ELSE meandata_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) / ntimes ENDIF ENDIF - - ! If this is the cell thickness, save to normalise thickness-weighted time-means - IF( jv == jv_thickness ) meancellthick_4d_sp(:,:,:,:) = meandata_4d_sp(:,:,:,:) + IF( l_thckwgt ) DEALLOCATE(meancellthick_4d_sp) CASE( NF90_DOUBLE ) IF( l_ismasked ) THEN IF( l_thckwgt ) THEN - IF( icellthick_type == NF90_FLOAT ) THEN - WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp - ELSEWHERE - meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (REAL(meancellthick_4d_sp, dp) * ntimes_4d(:,:,:,:)) - ENDWHERE - ELSE - WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp - ELSEWHERE - meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (meancellthick_4d_dp(:,:,:,:) * ntimes_4d(:,:,:,:)) - ENDWHERE - ENDIF + WHERE( meancellthick_4d_dp(:,:,:,:) == 0._dp ) + meandata_4d_dp(:,:,:,:) = outputdata_fill_value_dp + ELSEWHERE + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / meancellthick_4d_dp(:,:,:,:) + ENDWHERE ELSE WHERE( ntimes_4d(:,:,:,:) == 0 ) - meandata_4d_dp(:,:,:,:) = inputdata_fill_value_dp + meandata_4d_dp(:,:,:,:) = outputdata_fill_value_dp ELSEWHERE meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / ntimes_4d(:,:,:,:) ENDWHERE ENDIF ELSE IF( l_thckwgt ) THEN - IF( icellthick_type == NF90_FLOAT ) THEN - meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (REAL(meancellthick_4d_sp, dp) * ntimes) - ELSE - meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / (meancellthick_4d_dp(:,:,:,:) * ntimes) - ENDIF + meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / meancellthick_4d_dp(:,:,:,:) ELSE meandata_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) / ntimes ENDIF ENDIF - - ! If this is the cell thickness, save to normalise thickness-weighted time-means - IF( jv == jv_thickness ) meancellthick_4d_dp(:,:,:,:) = meandata_4d_dp(:,:,:,:) + IF( l_thckwgt ) DEALLOCATE(meancellthick_4d_dp) CASE DEFAULT WRITE(6,*)'Unknown nf90 type: ', xtype STOP 14 END SELECT - IF( l_ismasked ) DEALLOCATE(ntimes_4d) - + IF( l_ismasked .AND. .NOT. l_thckwgt ) DEALLOCATE(ntimes_4d) ENDIF ELSE @@ -1173,9 +1222,6 @@ PROGRAM mean_nemo END DO !loop over variables - IF( allocated(meancellthick_4d_sp) ) DEALLOCATE(meancellthick_4d_sp) - IF( allocated(meancellthick_4d_dp) ) DEALLOCATE(meancellthick_4d_dp) - !4.1 Close all input files IF (l_verbose) WRITE(6,*)'Closing input files...' From 13dd91ed7b80e47a01fbd5f7b7125973767cbffa Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Mon, 19 Jan 2026 19:47:19 +0000 Subject: [PATCH 3/9] Ensure all stop errors resulting directly from NetCDF Fortran API call errors are consistently reported --- Utilities/mean_nemo/mean_nemo.f90 | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index f32f45a..f9d0a53 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -156,7 +156,8 @@ PROGRAM mean_nemo iostat = nf90_open( TRIM(filenames(1)), nf90_share, ncid ) IF( iostat /= nf90_noerr ) THEN - WRITE(6,*) TRIM(nf90_strerror(iostat)) + WRITE(6,*)'E R R O R opening input file '//TRIM(filenames(1))//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 11 ENDIF iostat = nf90_inquire( ncid, ndims, nvars, natts ) @@ -283,8 +284,8 @@ PROGRAM mean_nemo DO ifile = 2, nargs-1 iostat = nf90_open( TRIM(filenames(ifile)), nf90_share, ncid, chunksize=chunksize) IF( iostat /= nf90_noerr ) THEN - WRITE(6,*) TRIM(nf90_strerror(iostat)) - WRITE(6,*)'E R R O R opening input file '//TRIM(filenames(ifile)) + WRITE(6,*)'E R R O R opening input file '//TRIM(filenames(ifile))//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 12 ELSE inncids(ifile) = ncid @@ -863,8 +864,8 @@ PROGRAM mean_nemo ENDIF !End of if statement over number of dimensions IF( iostat /= nf90_noerr ) THEN - WRITE(6,*) TRIM(nf90_strerror(iostat)) - WRITE(6,*) 'E R R O R reading variable '//TRIM(varname)//' from file '//TRIM(filenames(ifile)) + WRITE(6,*) 'E R R O R reading variable '//TRIM(varname)//' from file '//TRIM(filenames(ifile))//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) istop = 1 ENDIF @@ -1110,8 +1111,8 @@ PROGRAM mean_nemo ENDIF !End of ndims if statements IF( iostat /= nf90_noerr ) THEN - WRITE(6,*) TRIM(nf90_strerror(iostat)) - WRITE(6,*) 'E R R O R reading variable '//TRIM(varname)//' from file '//TRIM(filenames(ifile)) + WRITE(6,*) 'E R R O R reading variable '//TRIM(varname)//' from file '//TRIM(filenames(ifile))//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 16 ENDIF @@ -1215,8 +1216,9 @@ PROGRAM mean_nemo ENDIF - IF( iostat /= 0 ) THEN - WRITE(6,*) 'E R R O R writing variable '//TRIM(varname) + IF( iostat /= nf90_noerr ) THEN + WRITE(6,*) 'E R R O R writing variable '//TRIM(varname)//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 17 ENDIF @@ -1230,7 +1232,8 @@ PROGRAM mean_nemo iostat = nf90_close( ncid ) IF( iostat /= nf90_noerr ) THEN WRITE(6,*) TRIM(nf90_strerror(iostat)) - WRITE(6,*)'E R R O R closing input file '//TRIM(filenames(ifile)) + WRITE(6,*)'E R R O R closing input file '//TRIM(filenames(ifile))//':' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 18 ENDIF END DO @@ -1241,7 +1244,8 @@ PROGRAM mean_nemo iostat = nf90_close( outid ) IF( iostat /= nf90_noerr ) THEN WRITE(6,*) TRIM(nf90_strerror(iostat)) - WRITE(6,*)'E R R O R closing output file' + WRITE(6,*)'E R R O R closing output file:' + WRITE(6,*) ' '//TRIM(nf90_strerror(iostat)) STOP 19 ENDIF From bb899590c6c3b2b58b61adb3343f6e440bce3a8d Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Thu, 29 Jan 2026 15:27:36 +0000 Subject: [PATCH 4/9] Update compilation scripts --- Utilities/mean_nemo/compiler.xc40 | 42 ++++++++++++++++++++++--- Utilities/mean_nemo/compiler.xc40_login | 34 +++++++++++++++++--- 2 files changed, 67 insertions(+), 9 deletions(-) diff --git a/Utilities/mean_nemo/compiler.xc40 b/Utilities/mean_nemo/compiler.xc40 index 6a67dc8..ed80d48 100755 --- a/Utilities/mean_nemo/compiler.xc40 +++ b/Utilities/mean_nemo/compiler.xc40 @@ -1,7 +1,41 @@ #!/bin/bash +#PBS -N compile_mean_nemo +#PBS -o compile_mean_nemo.o +#PBS -e compile_mean_nemo.e +#PBS -q shared +#PBS -l ncpus=1 +#PBS -l walltime=00:01:00 +#PBS -l mem=500mb -rm -f mean_nemo.exe +set -eu -module unload cray-netcdf -module load cray-netcdf -ftn -O0 mean_nemo.f90 -o mean_nemo.exe +NAME=mean_nemo + +MODE=prod +#MODE=dev +#MODE=debug + +case ${MODE} in + prod) + opts="-O2 -Ovector1 -hfp0 -hflex_mp=intolerant" + ;; + dev) + opts="-O0" + ;; + debug) + opts="-O0 -Ovector0 -hflex_mp=intolerant -e CID -Ktrap=fp -g" + ;; + *) + echo "Compilation mode \"${MODE}\" not supported" + exit 1 + ;; +esac + +cd ${PBS_O_WORKDIR} + +[ -f ${NAME}.exe ] && rm -f ${NAME}.exe + +module unload cray-hdf5 cray-netcdf 2>/dev/null +module load cray-hdf5 cray-netcdf + +ftn ${opts} mean_nemo.f90 -o ${NAME}.exe diff --git a/Utilities/mean_nemo/compiler.xc40_login b/Utilities/mean_nemo/compiler.xc40_login index 343759f..2e211f0 100755 --- a/Utilities/mean_nemo/compiler.xc40_login +++ b/Utilities/mean_nemo/compiler.xc40_login @@ -1,8 +1,32 @@ #!/bin/bash -rm -f mean_nemo.exe +set -eu -module swap craype-haswell craype-ivybridge -module unload cray-netcdf -module load cray-netcdf -ftn -O0 mean_nemo.f90 -o mean_nemo_login.exe +NAME=mean_nemo_login + +MODE=prod +#MODE=dev +#MODE=debug + +case ${MODE} in + prod) + opts="-O2 -Ovector1 -hfp0 -hflex_mp=intolerant" + ;; + dev) + opts="-O0" + ;; + debug) + opts="-O0 -Ovector0 -hflex_mp=intolerant -e CID -Ktrap=fp -g" + ;; + *) + echo "Mode \"${MODE}\" not supported" + exit 1 + ;; +esac + +[ -f ${NAME}.exe ] && rm -f ${NAME}.exe + +module unload cray-hdf5 cray-netcdf 2>/dev/null +module load cray-hdf5 cray-netcdf + +ftn ${opts} mean_nemo.f90 -o ${NAME}.exe From a3896c146f0dcc71fef89d0ed5edb36a118f13c3 Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Thu, 29 Jan 2026 15:27:48 +0000 Subject: [PATCH 5/9] Remove unnecessary zeroing of input data arrays --- Utilities/mean_nemo/mean_nemo.f90 | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index f9d0a53..108ec7e 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -555,31 +555,26 @@ PROGRAM mean_nemo SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(inputdata_1d_i1(outdimlens(dimids(1)))) - inputdata_1d_i1(:)=0 iostat = nf90_get_var( ncid, jv, inputdata_1d_i1, start, indimlens) meandata_1d_i1(:)=meandata_1d_i1(:)+inputdata_1d_i1(:) DEALLOCATE(inputdata_1d_i1) CASE( NF90_SHORT ) ALLOCATE(inputdata_1d_i2(outdimlens(dimids(1)))) - inputdata_1d_i2(:)=0 iostat = nf90_get_var( ncid, jv, inputdata_1d_i2, start, indimlens) meandata_1d_i2(:)=meandata_1d_i2(:)+inputdata_1d_i2(:) DEALLOCATE(inputdata_1d_i2) CASE( NF90_INT ) ALLOCATE(inputdata_1d_i4(outdimlens(dimids(1)))) - inputdata_1d_i4(:)=0 iostat = nf90_get_var( ncid, jv, inputdata_1d_i4, start, indimlens) meandata_1d_i4(:)=meandata_1d_i4(:)+inputdata_1d_i4(:) DEALLOCATE(inputdata_1d_i4) CASE( NF90_FLOAT ) ALLOCATE(inputdata_1d_sp(outdimlens(dimids(1)))) - inputdata_1d_sp(:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_1d_sp, start, indimlens) meandata_1d_sp(:)=meandata_1d_sp(:)+inputdata_1d_sp(:) DEALLOCATE(inputdata_1d_sp) CASE( NF90_DOUBLE ) ALLOCATE(inputdata_1d_dp(outdimlens(dimids(1)))) - inputdata_1d_dp(:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_1d_dp, start, indimlens) meandata_1d_dp(:)=meandata_1d_dp(:)+inputdata_1d_dp(:) DEALLOCATE(inputdata_1d_dp) @@ -594,31 +589,26 @@ PROGRAM mean_nemo SELECT CASE( xtype ) CASE( NF90_BYTE ) ALLOCATE(inputdata_2d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)))) - inputdata_2d_i1(:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_2d_i1, start, indimlens) meandata_2d_i1(:,:)=meandata_2d_i1(:,:)+inputdata_2d_i1(:,:) DEALLOCATE(inputdata_2d_i1) CASE( NF90_SHORT ) ALLOCATE(inputdata_2d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)))) - inputdata_2d_i2(:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_2d_i2, start, indimlens ) meandata_2d_i2(:,:)=meandata_2d_i2(:,:)+inputdata_2d_i2(:,:) DEALLOCATE(inputdata_2d_i2) CASE( NF90_INT ) ALLOCATE(inputdata_2d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)))) - inputdata_2d_i4(:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_2d_i4, start, indimlens ) meandata_2d_i4(:,:)=meandata_2d_i4(:,:)+inputdata_2d_i4(:,:) DEALLOCATE(inputdata_2d_i4) CASE( NF90_FLOAT ) ALLOCATE(inputdata_2d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)))) - inputdata_2d_sp(:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_2d_sp, start, indimlens ) meandata_2d_sp(:,:)=meandata_2d_sp(:,:)+inputdata_2d_sp(:,:) DEALLOCATE(inputdata_2d_sp) CASE( NF90_DOUBLE ) ALLOCATE(inputdata_2d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)))) - inputdata_2d_dp(:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_2d_dp, start, indimlens ) meandata_2d_dp(:,:)=meandata_2d_dp(:,:)+inputdata_2d_dp(:,:) DEALLOCATE(inputdata_2d_dp) @@ -633,7 +623,6 @@ PROGRAM mean_nemo CASE( NF90_BYTE ) ALLOCATE(inputdata_3d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) - inputdata_3d_i1(:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_3d_i1, start, indimlens ) meandata_3d_i1(:,:,:)=meandata_3d_i1(:,:,:)+inputdata_3d_i1(:,:,:) DEALLOCATE(inputdata_3d_i1) @@ -642,7 +631,6 @@ PROGRAM mean_nemo CASE( NF90_SHORT ) ALLOCATE(inputdata_3d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) - inputdata_3d_i2(:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_3d_i2, start, indimlens ) meandata_3d_i2(:,:,:)=meandata_3d_i2(:,:,:)+inputdata_3d_i2(:,:,:) DEALLOCATE(inputdata_3d_i2) @@ -651,7 +639,6 @@ PROGRAM mean_nemo CASE( NF90_INT ) ALLOCATE(inputdata_3d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) - inputdata_3d_i4(:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_3d_i4, start, indimlens ) ! Do not include masked data in the average @@ -669,7 +656,6 @@ PROGRAM mean_nemo CASE( NF90_FLOAT ) ALLOCATE(inputdata_3d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) - inputdata_3d_sp(:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_3d_sp, start, indimlens ) ! Do not include masked data in the average @@ -687,7 +673,6 @@ PROGRAM mean_nemo CASE( NF90_DOUBLE ) ALLOCATE(inputdata_3d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)))) - inputdata_3d_dp(:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_3d_dp, start, indimlens ) ! Do not include masked data in the average @@ -713,7 +698,6 @@ PROGRAM mean_nemo CASE( NF90_BYTE ) ALLOCATE(inputdata_4d_i1(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - inputdata_4d_i1(:,:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_4d_i1, start, indimlens ) meandata_4d_i1(:,:,:,:)=meandata_4d_i1(:,:,:,:)+inputdata_4d_i1(:,:,:,:) DEALLOCATE(inputdata_4d_i1) @@ -722,7 +706,6 @@ PROGRAM mean_nemo CASE( NF90_SHORT ) ALLOCATE(inputdata_4d_i2(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - inputdata_4d_i2(:,:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_4d_i2, start, indimlens ) meandata_4d_i2(:,:,:,:)=meandata_4d_i2(:,:,:,:)+inputdata_4d_i2(:,:,:,:) DEALLOCATE(inputdata_4d_i2) @@ -731,7 +714,6 @@ PROGRAM mean_nemo CASE( NF90_INT ) ALLOCATE(inputdata_4d_i4(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - inputdata_4d_i4(:,:,:,:)=0 iostat = nf90_get_var( ncid, jv, inputdata_4d_i4, start, indimlens ) ! Do not include masked data in the average @@ -749,14 +731,12 @@ PROGRAM mean_nemo CASE( NF90_FLOAT ) ALLOCATE(inputdata_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - inputdata_4d_sp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_4d_sp, start, indimlens ) ! Thickness-weighting- get cell thickness IF( l_thckwgt ) THEN ALLOCATE(cellthick_4d_sp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - cellthick_4d_sp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv_thickness, cellthick_4d_sp, start, indimlens ) ENDIF @@ -802,14 +782,12 @@ PROGRAM mean_nemo CASE( NF90_DOUBLE ) ALLOCATE(inputdata_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - inputdata_4d_dp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv, inputdata_4d_dp, start, indimlens ) ! Thickness-weighting- get cell thickness IF( l_thckwgt ) THEN ALLOCATE(cellthick_4d_dp(outdimlens(dimids(1)),outdimlens(dimids(2)), & & outdimlens(dimids(3)),outdimlens(dimids(4)))) - cellthick_4d_dp(:,:,:,:)=0.0 iostat = nf90_get_var( ncid, jv_thickness, cellthick_4d_dp, start, indimlens ) ENDIF From 32d866398b127242ef1e5acd9559db0b2e1631bf Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Fri, 30 Jan 2026 11:50:36 +0000 Subject: [PATCH 6/9] Add basic timing calipers --- Utilities/mean_nemo/mean_nemo.f90 | 71 ++++++++++++++++++++++++++++++- 1 file changed, 69 insertions(+), 2 deletions(-) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index 108ec7e..f5932a3 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -26,8 +26,9 @@ PROGRAM mean_nemo INTEGER,PARAMETER :: sp=SELECTED_REAL_KIND(6,37) INTEGER,PARAMETER :: dp=SELECTED_REAL_KIND(12,307) - LOGICAL, PARAMETER :: l_verbose = .true. - + LOGICAL, PARAMETER :: l_verbose = .true. + LOGICAL, PARAMETER :: l_timing = .false. + CHARACTER(LEN=nf90_max_name) :: outfile, attname, dimname, varname, time, date, zone, timestamp CHARACTER(LEN=nf90_max_name), ALLOCATABLE :: filenames(:), indimnames(:) CHARACTER(LEN=256) :: standard_name,cell_methods @@ -126,7 +127,12 @@ PROGRAM mean_nemo REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: meandata_4d_dp REAL(dp), ALLOCATABLE, DIMENSION(:,:,:,:) :: meancellthick_4d_dp + ! Timing-related scalars + INTEGER(i8) :: t_section, t_total, t_variable + REAL(dp) :: secondclock + ! Initialise timing + CALL timing_init(secondclock) !End of definitions @@ -154,6 +160,9 @@ PROGRAM mean_nemo !--------------------------------------------------------------------------- !2. Read in the global dimensions from the first input file and set up the output file + CALL timing_start(t_total) ! Total time taken + CALL timing_start(t_section) + iostat = nf90_open( TRIM(filenames(1)), nf90_share, ncid ) IF( iostat /= nf90_noerr ) THEN WRITE(6,*)'E R R O R opening input file '//TRIM(filenames(1))//':' @@ -274,6 +283,7 @@ PROGRAM mean_nemo iostat = nf90_enddef( outid ) inncids(1) = ncid IF (l_verbose) WRITE(6,*)'Finished defining output file.' + CALL timing_stop(t_section, 'define output file') ; CALL timing_start(t_section) !--------------------------------------------------------------------------- !3. Read in data from each file for each variable @@ -292,9 +302,11 @@ PROGRAM mean_nemo ENDIF END DO IF (l_verbose) WRITE(6,*)'All input files open.' + CALL timing_stop(t_section, 'open input files') ; CALL timing_start(t_section) !Loop over all variables in first input file DO jv = 1, nvars + CALL timing_start(t_variable) ! Total time taken for variable !3.2 Inquire variable to find out name and how many dimensions it has @@ -518,11 +530,14 @@ PROGRAM mean_nemo STOP 15 ENDIF + CALL timing_stop(t_section, 'query input files and allocate arrays') + istop = 0 ! If this variable is a function of the unlimited dimension then ! Average over unlimited dimension IF( l_doavg ) THEN + CALL timing_start(t_section) DO ifile = 1, nargs-1 !Loop through input files @@ -855,6 +870,8 @@ PROGRAM mean_nemo END DO !loop over files + CALL timing_stop(t_section, 'variable '//TRIM(varname)//'- calculate sums for mean') ; CALL timing_start(t_section) + ! Divide numerator (sum of data, possibly thickness-weighted) by denominator (sum of cell thickness ! or number of time records) to get the mean, setting masked results to a fill value IF( ndims == 1 ) THEN @@ -1006,10 +1023,13 @@ PROGRAM mean_nemo IF( l_ismasked .AND. .NOT. l_thckwgt ) DEALLOCATE(ntimes_4d) ENDIF + CALL timing_stop(t_section, 'variable '//TRIM(varname)//'- calculate mean') + ELSE ! Else if the variable does not contain the unlimited dimension just read ! in from first file to be copied to outfile as it should be the same in all ! files (e.g. coordinates) + CALL timing_start(t_section) ncid = inncids(1) iostat = nf90_inquire_variable( ncid, jv, varname, xtype, ndims, dimids, natts) @@ -1094,10 +1114,13 @@ PROGRAM mean_nemo STOP 16 ENDIF + CALL timing_stop(t_section, 'variable '//TRIM(varname)//'- copy data without time coordinate') + ENDIF !End of check for unlimited dimension !--------------------------------------------------------------------------- !4. Write data to output file and close files + CALL timing_start(t_section) IF (l_verbose) WRITE(6,*)'Writing variable '//TRIM(varname)//'...' @@ -1200,9 +1223,13 @@ PROGRAM mean_nemo STOP 17 ENDIF + CALL timing_stop(t_section, 'variable '//TRIM(varname)//'- write to file') + CALL timing_stop(t_variable, 'variable '//TRIM(varname)//'- total time') ! Total time taken for variable + END DO !loop over variables !4.1 Close all input files + CALL timing_start(t_section) IF (l_verbose) WRITE(6,*)'Closing input files...' DO ifile = 1, nargs-1 @@ -1227,4 +1254,44 @@ PROGRAM mean_nemo STOP 19 ENDIF + CALL timing_stop(t_section, 'close files') + CALL timing_stop(t_total, 'TOTAL') ! Total time taken + + CONTAINS + + SUBROUTINE timing_init( secondclock ) + ! Timing initialisation- get the processor clock count rate + REAL(dp), INTENT(out) :: secondclock + INTEGER(i8) :: count_rate + + IF( .NOT. l_timing ) RETURN + + CALL SYSTEM_CLOCK(COUNT_RATE=count_rate) + secondclock = 1._dp / REAL(count_rate, dp) + END SUBROUTINE + + SUBROUTINE timing_start( count ) + ! Start a timing section + INTEGER(i8), INTENT(out) :: count + + IF( .NOT. l_timing ) RETURN + + CALL SYSTEM_CLOCK(COUNT=count) + END SUBROUTINE + + SUBROUTINE timing_stop( count, sect ) + ! Stop a timing section and report the time + INTEGER(i8), INTENT(in) :: count ! Counter returned by timing_start + CHARACTER(len=*), INTENT(in) :: sect + INTEGER(i8) :: count2 + REAL(dp) :: elapsed + CHARACTER(len=128) :: clfmt = "(a,f0.6,'s')" + + IF( .NOT. l_timing ) RETURN + + CALL SYSTEM_CLOCK(COUNT=count2) + elapsed = REAL(count2-count, dp) * secondclock + WRITE(6,clfmt) 'TIMING ('//TRIM(sect)//'): ', elapsed + END SUBROUTINE + END PROGRAM mean_nemo From f8164a3044f81254e9ce405c96bc76705742f972 Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Fri, 30 Jan 2026 13:34:58 +0000 Subject: [PATCH 7/9] Add script for generating a simple dataset for testing mean_nemo --- Utilities/mean_nemo/simple_test_data.py | 378 ++++++++++++++++++++++++ 1 file changed, 378 insertions(+) create mode 100755 Utilities/mean_nemo/simple_test_data.py diff --git a/Utilities/mean_nemo/simple_test_data.py b/Utilities/mean_nemo/simple_test_data.py new file mode 100755 index 0000000..2fe0c4a --- /dev/null +++ b/Utilities/mean_nemo/simple_test_data.py @@ -0,0 +1,378 @@ +#! /usr/bin/env python3 + +import os +import xarray as xr +import numpy as np + +""" +Script to generate a simple dataset for testing mean_nemo + +Two types of files will be generated- input data (intended as input to mean_nemo) and expected results +(the results expected from running mean_nemo with the input data). The following `nccmp` command is recommended +for comparing the latter file with the output from mean_nemo: + + `nccmp -F -f -c 1 -d -m ` + +The expected results are obtained by using xarray/numpy averaging methods. + +Both types of file will contain several variables that test different components of the mean_nemo code, named: + + {data_type}_{data_shape}_{weighting}_{time_coordinate}_{masking} + +where: + + * `data_type`- the data precision, using their numpy alias names ("byte", "short", "intc", "single", "double") + * `data_shape`- the shape of the data ("1d", "2d", "3d", "4d") + * `weighting`- the weighting method used in the average: + * "unweighted"- no weighting is used + * "weighted"- the average is weighted by an unmasked cell thickness + * "weighted-masked"- the average is weighted by a masked cell thickness + * `time_coordinate`- whether the data has a time coordinate ("time") or not ("notime") + * `masking`- whether the data is masked ("masked") or not ("unmasked") +""" + +# User input +# ============================================================================================================================ + +# Paths for files generated by this script +dir_out = '' # Output directory +file_out = f'{dir_out}/input.nc' # Input data (unmasked cell thickness) +file_out_cellthk_masked = f'{dir_out}/input_cellthk_masked.nc' # Input data (masked cell thickness) +file_expected = f'{dir_out}/expected.nc' # Expected results (unmasked cell thickness) +file_expected_cellthk_masked = f'{dir_out}/expected_cellthk_masked.nc' # Expected results (masked cell thickness) + +# Properties of the test data +data_dims = ('t', 'b', 'z', 'y', 'x') # Order of dimensions in the file +data_shape = {'t': 6, 'x': 10, 'y': 10, 'z': 10, 'b': 5} # Shape of the FULL test data (sliced by data_isel entries + # to generate the different data shapes) +data_time_split = [slice(0, 2), slice(2, 4), slice(4, 6)] # How data should be split in time between files +data_masked_isel = [ # Points to be masked in the FULL test data + {'x': [2, 4], 'y': [2, 4, 6]}, + {'z': [7, 8, 9]}, + {'t': [0], 'x': [0], 'y': [0], 'z': [0]}, + {'t': [2], 'x': [1], 'y': [0], 'z': [0]}, + {'t': [4], 'x': [1], 'y': [1], 'z': [0]}, + ] +cellthk_masked_isel = [ # Points to be masked in the FULL cell thickness data + {'x': [2], 'y': [4]}, + {'t': [0], 'x': [0], 'y': [1, 3, 5]}, + {'t': [3], 'x': [1], 'y': [1, 3, 5]}, + {'t': [5], 'x': [2], 'y': [1, 3, 5]}, + ] +data_isel = { # Data shapes to be generated and the coordinates + '1d': {'b': 0, 'z': 0, 'y': 0, 'x': 0}, # used to generate them with respect to the FULL data + '2d': {'b': 0, 'z': 0, 'y': 0}, + '3d': {'b': 0, 'z': 0}, + '4d': {'b': 0}, + # '5d': None, + } +data_fillvals = { # Data types to be generated and the fill values + 'byte': -1, # to be used for masked data + 'short': -1, + 'intc': -1, + 'single': 1e20, + 'double': 1e20, + } + +# ============================================================================================================================ + +# Numpy data type objects +data_dtypes = {i: getattr(np, i) for i in data_fillvals} + +# Dimensions +dict_da_dims = {k: xr.DataArray(range(v), dims=k) for k, v in data_shape.items()} + +# Full test data (double precision)- linear increase over each dimension +data_full = np.mgrid[*[slice(None, data_shape[i]) for i in data_dims]].sum(axis=0).astype('double') +da_full = xr.DataArray(data_full, + coords=[dict_da_dims[i] for i in data_dims]) + +# Mask for full test data (True where valid) +da_isvalid_full = xr.DataArray(True, coords=[dict_da_dims[i] for i in data_dims]) +for isel in data_masked_isel: + da_isvalid_full[isel] = False + +# Full cell thickness data (single precision)- linear increase with time +da_cellthk_full = xr.DataArray(0, + coords=[dict_da_dims[i] for i in ('t', 'z', 'y', 'x')], + attrs={'standard_name': 'cell_thickness'}, + name='cellthk').astype('single') +da_cellthk_full += da_cellthk_full.t + 1 + +# Mask for full cell thickness data (True where valid) +da_cellthk_isvalid_full = xr.DataArray(True, coords=[dict_da_dims[i] for i in ('t', 'z', 'y', 'x')]) +for isel in cellthk_masked_isel: + da_cellthk_isvalid_full[isel] = False + +# DataArray lists to be combined into Datasets later +da_list = [] +da_list_cellthk_masked = [] +da_expected_list = [] +da_expected_list_cellthk_masked = [] + + +def input_data(da: xr.DataArray, da_isvalid: xr.DataArray | None, + dtype: np.dtype, slices: dict | None, fillval: np.typing.ArrayLike, + drop_time_coord: bool = False, is_weighted: bool = False, masked_in_file: bool = False + ) -> (xr.DataArray, xr.DataArray | None): + """ + Generate input data. + + Parameters + ---------- + da : xr.DataArray + Full input data to use as the basis + da_isvalid : xr.DataArray | None + Mask to apply to `da` (`False` where masked) + dtype : np.dtype + Numpy dtype to convert `da` to + slices : dict | None + Indexing slices to apply to `da` + fillval : np.typing.ArrayLike + Value to set when applying `da_isvalid` mask + drop_time_coord : bool + Whether the time coordinate should be dropped from `da` + is_weighted : bool + Whether a `cell_methods` attribute will be applied to `da` to indicate it is weighted by the cell thickness + masked_in_file : bool + Whether a `_FillValue` attribute will be applied to `da` + + Returns + ------- + da : xr.DataArray + Input data + xr.DataArray | None + `masked_in_file == True:` + Input data mask (`True` where valid) corresponding to `da` + `masked_in_file == False`: + None + """ + # Slice full data and mask + if slices is not None: + da = da.isel(**slices, drop=True) + da_isvalid = da_isvalid.isel(**slices, drop=True) + + # Convert to given data type + da = da.astype(dtype) + + # Set masked data to given fill value (np.nan not representable for non-floats) + da = da.where(da_isvalid, fillval) + + # No time coordinate- use first time and drop coordinate + if drop_time_coord: + da = da.isel(t=0, drop=True) + + # Thickness-weighting- assign cell_methods recognised by mean_nemo + if is_weighted: + da = da.assign_attrs(cell_methods='time: mean (thickness weighted)') + + # Set _FillValue if data is to be read in as masked by mean_nemo + da.encoding = {'_FillValue': fillval if masked_in_file else None} + + return da, (da_isvalid if masked_in_file else None) + + +def cellthk_data(da: xr.DataArray, slices: dict | None = None, da_isvalid: xr.DataArray | None = None + ) -> xr.DataArray: + """ + Generate cell thickness data. + + Parameters + ---------- + da : xr.DataArray + Full cell thickness data to use as the basis + slices : dict | None + Indexing slices to apply to `da` + da_isvalid : xr.DataArray | None + Mask to apply to `da` (`False` where masked) + + Returns + ------- + da: xr.DataArray + Cell thickness data + """ + # If data will be masked, set masked points to np.nan and set a _FillValue for these points to be written to file + if da_isvalid is not None: + da = da.where(da_isvalid) + da.encoding = {'_FillValue': np.single(2e20)} + else: + da.encoding = {'_FillValue': None} + + # Slice full data + if slices is not None: + da = da.isel(**slices, drop=True, missing_dims='ignore') + + return da + + +def calc_mean(da: xr.DataArray, da_isvalid: xr.DataArray | None = None, da_cellthk: xr.DataArray | None = None + ) -> xr.DataArray: + """ + Calculate the time-average (i.e. the expected result) using xarray/numpy routines. + + Parameters + ---------- + da : xr.DataArray + Data to be averaged + da_isvalid : xr.DataArray | None + Mask to apply to `da` for the average (`False` where masked) + da_cellthk : xr.DataArray | None + Weighting to apply to `da` for the average (`np.nan` where masked) + + Returns + ------- + da: xr.DataArray + Time average of `da` + """ + dtype = da.dtype + fillval = da.encoding['_FillValue'] + attrs = da.attrs + + # Time coordinate to be added to average data + t_mean = da.t.mean().astype(np.intc).expand_dims('t') + + # Thickness-weighting + if da_cellthk is not None: + + # When the input data has no fill value, use the cell thickness fill value instead + if fillval is None: + fillval = da_cellthk.encoding['_FillValue'] + + # Set masked cell thickness data to 0 (exclude from weighted average- np.nan not allowed) + da_cellthk = da_cellthk.fillna(0) + + # Masked averaging + if da_isvalid is not None: + + # We can use the normal masking-aware aggregators for floats, which support np.nan + if dtype in (np.single, np.double): + if da_cellthk is not None: + da = da.where(da_isvalid.values).weighted(da_cellthk).mean('t') + else: + da = da.where(da_isvalid.values).mean('t') + + # Otherwise, we have to pass a `where` argument since np.nan isn't correctly represented for non-floats + else: + da = da.mean('t', where=da_isvalid.values, skipna=False) + + # The mean results in zeroes where there is no data over time, so we must restore the fill values + da = da.where(da_isvalid.any('t'), fillval) + + # Normal averaging + else: + if da_cellthk is not None: + + # The weighted mean will contain np.nan where there is no data, but this isn't supported for non-floats. + # Stop so we can see what needs to be done. + if dtype not in (np.single, np.double): + raise TypeError(f'Weighted unmasked average not supported for non-floats') + + da = da.weighted(da_cellthk).mean('t') + else: + da = da.mean('t') + + # Recast the result back to the original precision + da = da.astype(dtype) + + # Restore _FillValue and other attributes (averaging removes them) + da.encoding = {'_FillValue': fillval} + da = da.assign_attrs(**attrs) + + # Assign a time coordinate (dropped by averaging)- use an unweighted average since this is what mean_nemo does + da = da.expand_dims('t').assign_coords(t=t_mean) + + return da + + +# Generate data for each of the: +# Data types +for dtype_name in data_dtypes: + # Data shapes + for shape_name in data_isel: + # Data to be used in thickness-weighted (where cell thickness is masked or unmasked) or normal averages + for thk_name in ('unweighted', 'weighted', 'weighted-masked'): + # Data with/without a time coordinate + for time_name in ('time', 'notime'): + # Data that will be read in by mean_nemo as masked or unmasked + for mask_name in ('masked', 'unmasked'): + + dtype = data_dtypes[dtype_name] + fillval = dtype(data_fillvals[dtype_name]) + slices = data_isel[shape_name] + + # Skip thickness-weighted cases not supported by mean_nemo + if (thk_name.startswith('weighted') and + not (time_name == 'time' and shape_name in ('4d',) and dtype in (np.single, np.double))): + continue + # Skip scalars, which aren't supported by mean_nemo + if (time_name == 'notime') and (shape_name == '1d'): + continue + # Skip masked cases not supported by mean_nemo + if (mask_name == 'masked') and not (shape_name in ('3d', '4d') and dtype in (np.intc, np.single, np.double)): + continue + + # Get input data and its mask + da_input, da_isvalid = input_data(da_full, da_isvalid_full, dtype, slices, fillval, + drop_time_coord=(time_name == 'notime'), + is_weighted=(thk_name.startswith('weighted')), + masked_in_file=(mask_name == 'masked'), + ) + + # Variable name + da_input.name = f'{dtype_name}_{shape_name}_{thk_name}_{time_name}_{mask_name}' + + # Get the cell thickness if required for thickness-weighting + if thk_name.startswith('weighted'): + da_cellthk = cellthk_data(da_cellthk_full, slices=slices, + da_isvalid=(da_cellthk_isvalid_full if thk_name == 'weighted-masked' else None), + ) + else: + da_cellthk = None + + # Calculate the expected result + if time_name == 'time': + da_expected = calc_mean(da_input, da_isvalid=da_isvalid, da_cellthk=da_cellthk) + else: + da_expected = da_input + + # We can only have one cell thickness variable per file (mean_nemo limitation)- + # put the cell thickness and the variables weighted by it in separate files + if thk_name == 'weighted-masked': + da_list_cellthk_masked.append(da_input) + da_expected_list_cellthk_masked.append(da_expected) + else: + da_list.append(da_input) + da_expected_list.append(da_expected) + +# Add unmasked cell thickness to data +da_cellthk = cellthk_data(da_cellthk_full) +da_expected = calc_mean(da_cellthk) +da_list.append(da_cellthk) +da_expected_list.append(da_expected) + +# Add masked cell thickness to data +da_cellthk = cellthk_data(da_cellthk_full, da_isvalid=da_cellthk_isvalid_full) +da_expected = calc_mean(da_cellthk) +da_list_cellthk_masked.append(da_cellthk) +da_expected_list_cellthk_masked.append(da_expected) + +# Combine DataArrays to Datasets +ds_out = xr.merge(da_list) +ds_out_cellthk_masked = xr.merge(da_list_cellthk_masked) +ds_expected = xr.merge(da_expected_list) +ds_expected_cellthk_masked = xr.merge(da_expected_list_cellthk_masked) + +# Create output directory if needed +if not os.path.isdir(dir_out): + os.makedirs(dir_out) + +# Write input data to files- potentially split over multiple files (file_0.nc, file_1.nc, etc) +for i_file, time_slice in enumerate(data_time_split): + fname = f'{file_out.replace('.nc', '')}_{i_file}.nc' + ds_out.isel(t=time_slice).to_netcdf(fname, mode='w', format='NETCDF4_CLASSIC', unlimited_dims=['t']) + fname = f'{file_out_cellthk_masked.replace('.nc', '')}_{i_file}.nc' + ds_out_cellthk_masked.isel(t=time_slice).to_netcdf(fname, mode='w', format='NETCDF4_CLASSIC', unlimited_dims=['t']) + +# Write expected results to files- currently NC3 with 64-bit offset +ds_expected.to_netcdf(file_expected, mode='w', format='NETCDF3_64BIT', unlimited_dims=['t']) +ds_expected_cellthk_masked.to_netcdf(file_expected_cellthk_masked, mode='w', format='NETCDF3_64BIT', unlimited_dims=['t']) From 527bac03b8e2d96652d81215a568a97e68520d44 Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Wed, 4 Feb 2026 18:30:40 +0000 Subject: [PATCH 8/9] Update CONTRIBUTORS.md --- CONTRIBUTORS.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 2278e5c..794d23d 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -1,6 +1,7 @@ # Contributors -| GitHub user | Real Name | Affiliation | Date | -| ----------- | --------- | ----------- | ---- | -| james-bruten-mo | James Bruten | Met Office | 2025-12-09 | -| Pierre-siddall | Pierre Siddall | Met Office | 2025-01-19 | +| GitHub user | Real Name | Affiliation | Date | +|-----------------|----------------|-------------|------------| +| james-bruten-mo | James Bruten | Met Office | 2025-12-09 | +| Pierre-siddall | Pierre Siddall | Met Office | 2025-01-19 | +| dcalve | Daley Calvert | Met Office | 2026-02-04 | From f29bd2383749c0cfd792d4ebd67e63df2b358905 Mon Sep 17 00:00:00 2001 From: Daley Calvert Date: Fri, 6 Feb 2026 15:26:24 +0000 Subject: [PATCH 9/9] Add Momentum Copyright templates --- Utilities/mean_nemo/mean_nemo.f90 | 6 ++++++ Utilities/mean_nemo/simple_test_data.py | 6 ++++++ 2 files changed, 12 insertions(+) diff --git a/Utilities/mean_nemo/mean_nemo.f90 b/Utilities/mean_nemo/mean_nemo.f90 index f5932a3..c2013fe 100644 --- a/Utilities/mean_nemo/mean_nemo.f90 +++ b/Utilities/mean_nemo/mean_nemo.f90 @@ -1,3 +1,9 @@ +! ----------------------------------------------------------------------------- +! (C) Crown copyright Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +! ----------------------------------------------------------------------------- + PROGRAM mean_nemo !----------------------------------------------------------------------------- diff --git a/Utilities/mean_nemo/simple_test_data.py b/Utilities/mean_nemo/simple_test_data.py index 2fe0c4a..aa7fd6d 100755 --- a/Utilities/mean_nemo/simple_test_data.py +++ b/Utilities/mean_nemo/simple_test_data.py @@ -1,5 +1,11 @@ #! /usr/bin/env python3 +# ----------------------------------------------------------------------------- +# (C) Crown copyright Met Office. All rights reserved. +# The file LICENCE, distributed with this code, contains details of the terms +# under which the code may be used. +# ----------------------------------------------------------------------------- + import os import xarray as xr import numpy as np