From 8208f290b2ae01656096f20f7fd45fd8d90438c0 Mon Sep 17 00:00:00 2001 From: oceandie Date: Tue, 16 Jan 2024 17:06:45 +0000 Subject: [PATCH 01/15] starting modifying the README to describe the purpose of this branch --- README.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 5f92759..eb63b32 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,7 @@ -# CDFTOOLS +# CDFTOOLS - JMMP fork CDFTOOLS is a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). It is now available on GitHub under the CeCILL license (). - CDFTOOLS is an open source package and contributions from other group are welcomed. The Git workflow policy is still to be defined. - - Actual version of CDFTOOLS is known as version 4. (See changes in paragraph *New features in version 4*, below). + This is the main development branch of the JMMP FORK of the [main repository of CDFTOOLS](https://github.com/meom-group/CDFTOOLS). The actual version of CDFTOOLS is known as version 4. ## Using CDFTOOLS From 74d6252c545fb0b2c23544929930dbebeb21a90d Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 16 Jan 2024 17:14:40 +0000 Subject: [PATCH 02/15] Update README.md --- README.md | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index eb63b32..2801d78 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,6 @@ -# CDFTOOLS - JMMP fork - CDFTOOLS is a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). It is now available on GitHub under the CeCILL license (). - - This is the main development branch of the JMMP FORK of the [main repository of CDFTOOLS](https://github.com/meom-group/CDFTOOLS). The actual version of CDFTOOLS is known as version 4. - +# CDFTOOLS - JMMP development branch + This is the main development branch of the JMMP FORK of [CDFTOOLS](https://github.com/meom-group/CDFTOOLS), a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). The main repository is available on GitHub () under the CeCILL license (). The current supported version of CDFTOOLS is version 4. + ## Using CDFTOOLS #### Cloning the git repository From c7cb01e68830240f251b8949c22bffbdbc753f48 Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 16 Jan 2024 17:19:24 +0000 Subject: [PATCH 03/15] Update README.md --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 2801d78..a65dd5c 100644 --- a/README.md +++ b/README.md @@ -4,9 +4,9 @@ ## Using CDFTOOLS #### Cloning the git repository -To retrieve a copy of the CDFTOOLS source code and create a working directory, run the following on the command line: +To retrieve a local copy of the CDFTOOLS `dev_jmmp` development branch run the following on the command line: -```> git clone https://github.com/meom-group/CDFTOOLS ``` +```> git clone -b dev_jmmp --single-branch https://github.com/JMMP-Group/CDFTOOLS.git CDFTOOLS_jmmp``` #### Compiling CDFTOOLS All the fortran source are in src/ subdirectory. In src/ there is a Makefile for compiling the sources. The compiler/machines related definitions are supposed to be collected in a `make.macro` file. Some examples of `make.macro` are given in the Macrolib directory and can be used as template for a new compiler or new machine. Then the good practice is to make a link From f98bdf21839b5c488b05ca9b400065a5c0b2a0f0 Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 16 Jan 2024 17:29:01 +0000 Subject: [PATCH 04/15] Update README.md --- README.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index a65dd5c..7ebf72c 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,9 @@ # CDFTOOLS - JMMP development branch - This is the main development branch of the JMMP FORK of [CDFTOOLS](https://github.com/meom-group/CDFTOOLS), a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). The main repository is available on GitHub () under the CeCILL license (). The current supported version of CDFTOOLS is version 4. + This is the main development branch of the JMMP FORK of [CDFTOOLS](https://github.com/meom-group/CDFTOOLS), a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). The main repository of CDFTOOLS is available on GitHub () under the CeCILL license (). The current supported version of CDFTOOLS is version 4. -## Using CDFTOOLS +## Using the JMMP development branch of CDFTOOLS + +The following are instructions specific for using the JMMP development branch of CDFTOOLS in the main computational resources of the UK used for numerical ocean modelling (e.g., Met Office HPC, NOC, ARCHER, JASMIN ...). These instructions closely follow the more general instructions provided on the main repository of CDFTOOLS (). #### Cloning the git repository To retrieve a local copy of the CDFTOOLS `dev_jmmp` development branch run the following on the command line: From d974bcb195a23b036df4d5cdb16adce8dad3e801 Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 16 Jan 2024 17:34:02 +0000 Subject: [PATCH 05/15] Update README.md --- README.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 7ebf72c..8e3732c 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,9 @@ # CDFTOOLS - JMMP development branch This is the main development branch of the JMMP FORK of [CDFTOOLS](https://github.com/meom-group/CDFTOOLS), a diagnostic package written in fortran 90 for the analysis of NEMO model output, initialized in the frame of the DRAKKAR project (). The main repository of CDFTOOLS is available on GitHub () under the CeCILL license (). The current supported version of CDFTOOLS is version 4. -## Using the JMMP development branch of CDFTOOLS +## Using CDFTOOLS -The following are instructions specific for using the JMMP development branch of CDFTOOLS in the main computational resources of the UK used for numerical ocean modelling (e.g., Met Office HPC, NOC, ARCHER, JASMIN ...). These instructions closely follow the more general instructions provided on the main repository of CDFTOOLS (). - -#### Cloning the git repository +#### Cloning the git JMMP development branch To retrieve a local copy of the CDFTOOLS `dev_jmmp` development branch run the following on the command line: ```> git clone -b dev_jmmp --single-branch https://github.com/JMMP-Group/CDFTOOLS.git CDFTOOLS_jmmp``` From 0ab0d51602290527ccbb53c25e479712d59935e5 Mon Sep 17 00:00:00 2001 From: DaveStorkey Date: Fri, 26 Jul 2024 10:01:06 +0100 Subject: [PATCH 06/15] Bug fixes (missing variable declarations/definitions) for modcdfnames_CMIP6.h90. --- src/modcdfnames_CMIP6.h90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/modcdfnames_CMIP6.h90 b/src/modcdfnames_CMIP6.h90 index b5350d3..4a42ee8 100644 --- a/src/modcdfnames_CMIP6.h90 +++ b/src/modcdfnames_CMIP6.h90 @@ -32,12 +32,15 @@ ! VVL case [ In CPMIP so far, only one name is used. ] CHARACTER(LEN=256) :: cn_ve3tvvl='thkcello', cn_ve3wvvl='thkcello' !: e3. (3D). CHARACTER(LEN=256) :: cn_ve3uvvl='thkcello', cn_ve3vvvl='thkcello' !: e3. + CHARACTER(LEN=256) :: cn_ve3t0='e3t_0', cn_ve3w0='e3w_0' !: e3. (3D). (at rest) + CHARACTER(LEN=256) :: cn_ve3u0='e3u_0', cn_ve3v0='e3v_0' !: e3. CHARACTER(LEN=256) :: cn_vff='ff' CHARACTER(LEN=256) :: cn_gdept='gdept', cn_gdepw='gdepw' !: 1d dep variable CHARACTER(LEN=256) :: cn_hdept='hdept', cn_hdepw='hdepw' !: 2d dep variable + CHARACTER(LEN=256) :: cn_dept3d='gdept_0' !: initial dept 3D CHARACTER(LEN=256) :: cn_depu3d='depu3d', cn_depw3d='depw3d' !: Local depth U and W in broken line extraction CHARACTER(LEN=256) :: cn_glamt='glamt', cn_gphit='gphit' !: glam gphi From 447a49be6849bae508e7bcaff32efe569c527a97 Mon Sep 17 00:00:00 2001 From: DaveStorkey Date: Tue, 6 Aug 2024 08:12:40 +0100 Subject: [PATCH 07/15] Committing the version of cdfsmooth.f90 with Mike's updated Shapiro filter. Note this still requires some work, but it runs and produces a smoothed field. --- src/cdfsmooth.f90 | 465 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 385 insertions(+), 80 deletions(-) diff --git a/src/cdfsmooth.f90 b/src/cdfsmooth.f90 index 34f7339..0f1a569 100644 --- a/src/cdfsmooth.f90 +++ b/src/cdfsmooth.f90 @@ -18,6 +18,12 @@ PROGRAM cdfsmooth !! 3.0 : 07/2011 : R. Dussin : Add anisotropic box !! : 4.0 : 03/2017 : J.M. Molines !!---------------------------------------------------------------------- + !! *** ROUTINE lisshapiro2d *** Mike Bell 17 Aug 2021 + !! + !! ** Purpose : apply a korder 2D shapiro filter kpass times. The land/sea mask and values at + !! The land/sea mask and values at selected points may be forced iteratively toward their initial values. + !! + !! !!---------------------------------------------------------------------- !! routines : description !! filterinit : initialise weight @@ -28,7 +34,7 @@ PROGRAM cdfsmooth !! initbox : initialize weight for box car average !! lislanczos2d : Lanczos filter !! lishan2d : hanning 2d filter - !! lisshapiro1d : shapiro filter + !! lisshapiro2d : shapiro filter !! lisbox : box car filter !!---------------------------------------------------------------------- USE cdfio @@ -53,6 +59,7 @@ PROGRAM cdfsmooth INTEGER(KIND=4) :: narg, iargc ! browse arguments INTEGER(KIND=4) :: ijarg ! argument index for browsing line INTEGER(KIND=4) :: ncut, nband ! cut period/ length, bandwidth + INTEGER(KIND=4) :: npass ! number of passes of Shapiro filter INTEGER(KIND=4) :: nfilter = jp_lanc ! default value INTEGER(KIND=4) :: nvars, ierr ! number of vars INTEGER(KIND=4) :: ncout ! ncid of output file @@ -66,6 +73,8 @@ PROGRAM cdfsmooth REAL(KIND=4) :: fn, rspval ! cutoff freq/wavelength, spval REAL(KIND=4) :: ranis ! anistropy + LOGICAL :: lap ! choice 1D or 2D (Laplacian) for Shapiro filter + LOGICAL :: lsm ! T => land points not used (for Shapiro filter only) REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdep, gdeptmp ! depth array REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: v2d, w2d ! raw data, filtered result @@ -77,20 +86,22 @@ PROGRAM cdfsmooth CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cv_names ! array of var name CHARACTER(LEN=256) :: cf_in, cf_out ! file names + CHARACTER(LEN=256) :: cf_in_stem, cf_suffix ! file name parts CHARACTER(LEN=256) :: cv_dep, cv_tim ! variable name for depth and time CHARACTER(LEN=256) :: ctyp ! filter type CHARACTER(LEN=256) :: cldum ! dummy character variable CHARACTER(LEN=256) :: clklist ! ciphered k-list of level LOGICAL :: lnc4 = .FALSE. ! flag for netcdf4 output with chinking and deflation + INTEGER(KIND=4) :: ji, jj, jmx, jkx ! dummy loop index !!---------------------------------------------------------------------- CALL ReadCdfNames() narg=iargc() IF ( narg == 0 ) THEN - PRINT *,' usage : cdfsmooth -f IN-file -c ncut [-t FLT-type] [-k LST-level] ...' - PRINT *,' [-anis ratio ] [-nc4 ] ' + PRINT *,' usage : cdfsmooth -f IN-file -c ncut [-t FLT-type] [-npass npass] [-k LST-level] ...' + PRINT *,' [-anis ratio ] [-lap lap] [-lsm lsm] [-nc4 ] ' PRINT *,' ' PRINT *,' PURPOSE :' PRINT *,' Perform a spatial smoothing on the file using a particular filter as' @@ -107,17 +118,20 @@ PROGRAM cdfsmooth PRINT *,' Hanning , H, h' PRINT *,' Shapiro , S, s' PRINT *,' Box , B, b' - PRINT *,' [-anis ratio ] : Specify an anisotropic ratio in case of Box-car filter.' - PRINT *,' With ratio=1, the box is a square 2.ncut x 2.ncut grid points.' - PRINT *,' In general, the box is then a rectangle 2.ncut*ratio x 2.ncut.' + PRINT *,' [-npass npass ] : Number of passes of filter; only used for Shapiro filter' PRINT *,' [-k LST-level ] : levels to be filtered (default = all levels)' PRINT *,' LST-level is a comma-separated list of levels. For example,' PRINT *,' the syntax 1-3,6,9-12 will select 1 2 3 6 9 10 11 12' + PRINT *,' [-anis ratio ] : Specify an anisotropic ratio in case of Box-car filter.' + PRINT *,' With ratio=1, the box is a square 2.ncut x 2.ncut grid points.' + PRINT *,' In general, the box is then a rectangle 2.ncut*ratio x 2.ncut.' + PRINT *,' [-lap lap ] : .TRUE. implies Laplacian; .FALSE. implies 1D lines; only used for Shapiro filter' + PRINT *,' [-lsm lsm ] : .TRUE. implies land pts not used; .FALSE. implies full 2D field (no land/sea mask); only used for Shapiro filter' PRINT *,' [-nc4] : produce netcdf4 output file with chunking and deflation.' PRINT *,' ' PRINT *,' OUTPUT : ' PRINT *,' Output file name is build from input file name with indication' - PRINT *,' of the filter type (1 letter) and of ncut.' + PRINT *,' of the filter type (1 letter) and of ncut .' PRINT *,' netcdf file : IN-file[LHSB]ncut' PRINT *,' variables : same as input variables.' PRINT *,' ' @@ -127,17 +141,23 @@ PROGRAM cdfsmooth ijarg = 1 ilev = 0 ranis = 1 ! anisotropic ratio for Box car filter + lap = .TRUE. + lsm = .TRUE. ctyp = 'L' ncut = 0 ! hence program exit if none specified on command line + npass = 1 ! just one pass DO WHILE (ijarg <= narg ) CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 SELECT CASE (cldum) CASE ( '-f' ) ; CALL getarg ( ijarg, cf_in ) ; ijarg=ijarg+1 CASE ( '-c' ) ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) ncut CASE ( '-t' ) ; CALL getarg ( ijarg, ctyp ) ; ijarg=ijarg+1 + CASE ( '-npass') ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) npass CASE ( '-k' ) ; CALL getarg ( ijarg, clklist ) ; ijarg=ijarg+1 ; CALL GetList (clklist, iklist, ilev ) CASE ('-anis') ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) ranis + CASE ('-lap') ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) lap + CASE ('-lsm') ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) lsm CASE ( '-nc4') ; lnc4 = .TRUE. CASE DEFAULT ; PRINT *,' ERROR :' ,TRIM(cldum),' : unknown option.' ; STOP 99 END SELECT @@ -147,6 +167,13 @@ PROGRAM cdfsmooth ENDIF IF ( chkfile(cf_in) ) STOP 99 ! missing file + IF ( index( cf_in, ".nc" ) > 0 ) THEN + cf_in_stem = cf_in(1:len(trim(cf_in))-3) + cf_suffix = ".nc" + ELSE + cf_in_stem = cf_in + cf_suffix = "" + ENDIF ! remark: for a spatial filter, fn=dx/lambda where dx is spatial step, lamda is cutting wavelength fn = 1./ncut @@ -154,25 +181,25 @@ PROGRAM cdfsmooth ALLOCATE ( dec(0:nband) , de(0:nband) ) - WRITE(cf_out,'(a,a,i3.3)') TRIM(cf_in),'L',ncut ! default name + WRITE(cf_out,'(a,a,i3.3,a)') TRIM(cf_in_stem),'_L',ncut,TRIM(cf_suffix) ! default name SELECT CASE ( ctyp) CASE ( 'Lanczos','L','l') nfilter=jp_lanc - WRITE(cf_out,'(a,a,i3.3)') TRIM(cf_in),'L',ncut + WRITE(cf_out,'(a,a,i3.3,a)') TRIM(cf_in_stem),'_L',ncut,TRIM(cf_suffix) PRINT *,' Working with Lanczos filter' CASE ( 'Hanning','H','h') nfilter=jp_hann ALLOCATE ( dec2d(0:2,0:2) ) - WRITE(cf_out,'(a,a,i3.3)') TRIM(cf_in),'H',ncut + WRITE(cf_out,'(a,a,i3.3,a)') TRIM(cf_in_stem),'_H',ncut,TRIM(cf_suffix) PRINT *,' Working with Hanning filter' CASE ( 'Shapiro','S','s') nfilter=jp_shap - WRITE(cf_out,'(a,a,i3.3)') TRIM(cf_in),'S',ncut + WRITE(cf_out,'(a,a,2i1.1,2l1,a)') TRIM(cf_in_stem),'_S',ncut,npass,lap,lsm,TRIM(cf_suffix) PRINT *,' Working with Shapiro filter' CASE ( 'Box','B','b') nfilter=jp_boxc - WRITE(cf_out,'(a,a,i3.3)') TRIM(cf_in),'B',ncut + WRITE(cf_out,'(a,a,i3.3,a)') TRIM(cf_in_stem),'_B',ncut,TRIM(cf_suffix) IF ( ranis /=1. ) THEN PRINT *, 'Anisotropic box car with ratio Lx = ', ranis, 'x Ly' ELSE @@ -267,7 +294,6 @@ PROGRAM cdfsmooth WHERE ( v2d == rspval ) iw =0 IF ( ncut /= 0 ) CALL filter( nfilter, v2d, iw, w2d) IF ( ncut == 0 ) w2d = v2d - w2d = w2d *iw ! mask filtered data ierr = putvar(ncout, id_varout(jvar), w2d, jk, npiglo, npjglo, ktime=jt) ! END DO @@ -314,7 +340,7 @@ SUBROUTINE filter (kfilter, px, kpx, py) SELECT CASE ( kfilter) CASE ( jp_lanc ) ; CALL lislanczos2d (px, kpx, py, npiglo, npjglo, fn, nband) CASE ( jp_hann ) ; CALL lishan2d (px, kpx, py, ncut, npiglo, npjglo) - CASE ( jp_shap ) ; CALL lisshapiro1d (px, kpx, py, ncut, npiglo, npjglo) + CASE ( jp_shap ) ; CALL lisshapiro2d (px, kpx, py, ncut, npass, lap, lsm, npiglo, npjglo) CASE ( jp_boxc ) ; CALL lisbox (px, kpx, py, npiglo, npjglo, fn, nband, ranis) END SELECT @@ -531,86 +557,296 @@ SUBROUTINE lishan2d(px, kiw, py, korder, kpi, kpj) END SUBROUTINE lishan2d - SUBROUTINE lisshapiro1d(px, kiw, py, korder, kpi, kpj) + SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) !!--------------------------------------------------------------------- - !! *** ROUTINE lisshapiro1d *** + !! *** ROUTINE lisshapiro2d *** !! - !! ** Purpose : compute shapiro filter + !! ** Purpose : apply a korder 2D shapiro filter kpass times. The land/sea mask and + !! values at selected points may be forced to their initial values. !! - !! References : adapted from Mercator code + !! References : Originally adapted from Mercator code. + !! This version rewritten by Mike Bell to implement the algorithm from + !! Francis, P.E., Quart. J. Roy. Met. Soc 101, pp 567-582 (1975) !!---------------------------------------------------------------------- REAL(KIND=4), DIMENSION(:,:), INTENT(in ) :: px ! input data INTEGER(KIND=4), DIMENSION(:,:), INTENT(in ) :: kiw ! validity flags REAL(KIND=4), DIMENSION(:,:), INTENT(out) :: py ! output data INTEGER(KIND=4), INTENT(in ) :: korder ! order of the filter + INTEGER(KIND=4), INTENT(in ) :: kpass ! number of passes of the filter + LOGICAL, INTENT(in ) :: lap ! True => Laplacian; False => lines + LOGICAL, INTENT(in ) :: lsm ! True => land points not used; False => no land/sea mask used (full 2D field) INTEGER(KIND=4), INTENT(in ) :: kpi, kpj ! size of the data - INTEGER(KIND=4) :: jj, ji, jorder ! loop indexes - INTEGER(KIND=4) :: imin, imax, ihalo=0 - REAL(KIND=4), PARAMETER :: rp_aniso_diff_XY = 2.25 ! anisotrope case - REAL(KIND=4) :: zalphax, zalphay, znum - REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: ztmp , zpx , zpy, zkiw - LOGICAL :: ll_cycl = .TRUE. - !!---------------------------------------------------------------------- + INTEGER(KIND=4) :: max_iterations = 10 ! max number of iterations allowed to enforce original land/sea mask and original values at selected points + + INTEGER(KIND=4) :: jn_fix_pts ! number of points where bathymetry is to remain fixed + INTEGER(KIND=4) :: jj, ji, jorder, jpass, jiteration, jpt ! loop indexes + INTEGER(KIND=4) :: ijt ! transposed ji index used for north pole fold + INTEGER(KIND=4) :: ji_min, jj_min ! indices of min field values + INTEGER(KIND=4) :: jcount_shallow, jcount_fixed ! temporary indices for print out + REAL(KIND=4) :: znum + REAL(KIND=4) :: rms_int, znpts + REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: ztmp , zpx , zpx_iteration, zpy, zkiw, zones + INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ji_fix ! i indices of bathymetry pts to hold fixed + INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: jj_fix ! j indices of bathymetry pts to hold fixed + LOGICAL :: l_test_shallow, l_test_fixed ! temporary logicals + +!----------------------------------------------------------------------------------------------- +!! Variables that can be set through the namelist + + LOGICAL :: ll_npol_fold ! north pole fold in grid? + LOGICAL :: ll_cycl ! this filter has only been tested for cyclic grids +!! for limited area models the bathymetry should be smoothed over a wider domain than that of the model; the margin should be at least korder*kpass +!! points and preferably 2*korder*kpass points (the bathymetry in nested models needs to match that of the outer domain in the nesting zone; that +!! might be done by including the nesting zone in the list of fixed points - but simpler solutions might be adequate) + + LOGICAL :: l_single_point_response ! => study response to single non-zero value + LOGICAL :: l_pass_shallow_updates ! => update values where bathy < zmin_val between passes + LOGICAL :: l_pass_fixed_pt_updates ! => update values where bathymetry should remain fixed + + REAL(KIND=4) :: zmin_val, zfactor_shallow, ztol_fixed, ztol_shallow ! see below + + INTEGER (KIND=4) :: ji_single_pt, jj_single_pt ! indices of single point (see l_single_point_response) + INTEGER(KIND=4) :: jst_prt, jend_prt ! + +!! The points that are printed out are set by ji_min_prt, jj_min_prt, ji_max_prt, jj_max_prt - these are local to prt_summary + +!----------------------------------------------------------------------------------------------- + + NAMELIST / nam_shapiro / ll_npol_fold, ll_cycl, l_single_point_response, l_pass_shallow_updates, l_pass_fixed_pt_updates, & + & ji_single_pt, jj_single_pt, jst_prt, jend_prt +!! Namelist default values + ll_npol_fold = .FALSE. + ll_cycl = .FALSE. + l_single_point_response = .FALSE. + l_pass_shallow_updates = .TRUE. + l_pass_fixed_pt_updates = .TRUE. + + zmin_val = -5 ! minimum depth (e.g. 10.0 metres) + ztol_shallow = 1.0 ! tolerance in metres of minimum shallow values (zmin_val - ztol_shallow) + ztol_fixed = 1.0 ! tolerance in metres for fit to bathymetry at point specified to be fixed + zfactor_shallow = 1.5 ! zfactor_shallow needs to be slightly greater than 1.0 ! 1.1 to 1.5 are reasonable values + + ji_single_pt = 622 ; jj_single_pt = 779 + jst_prt = 400 ; jend_prt = 405 + + OPEN(UNIT=20, file = 'namelist_shapiro.txt', form='formatted', status='old' ) + READ(NML=nam_shapiro, UNIT = 20) + WRITE(NML=nam_shapiro, UNIT=6) + CLOSE(20) + + PRINT*, 'korder, kpass, lap, lsm = ', korder, kpass, lap, lsm + +!----------------------------------------------------------------------------------------------- + +! MJB 2021/02/10: This code has been re-written for files in which column 1 is the same as column kpi - 1; col 2 is same as col kpi. +! In other words the halo columns are included in the input file. +! px + + ! we do NOT allocate with an additional ihalo + ALLOCATE( ztmp(kpi,kpj) , zkiw(kpi,kpj) ) + ALLOCATE( zpx (kpi,kpj) , zpy (kpi,kpj) ) + ALLOCATE( zones(kpi,kpj), zpx_iteration(kpi,kpj) ) + + ! print values from top and bottom rows to check they are OK + PRINT *, ' col 1 = ', (px(1,jj), jj = jst_prt, jend_prt) + PRINT *, ' col 2 = ', (px(2,jj), jj = jst_prt, jend_prt) + PRINT *, ' col kpi-1 = ', (px(kpi-1,jj), jj = jst_prt, jend_prt) + PRINT *, ' col kpi = ', (px(kpi,jj), jj = jst_prt, jend_prt) + +! option for testing out response to a single non-zero value + IF ( l_single_point_response ) THEN + zpx(:,:) = 0.0 + zpx(ji_single_pt,jj_single_pt) = 1.0 ! choose a point somewhere in the domain + END IF + +! read the indices of points to remain fixed + IF ( l_pass_fixed_pt_updates ) THEN + OPEN (unit=20, file = 'Notes/list_fixed_points.txt', form='formatted', status='old' ) + READ (20, *) jn_fix_pts + + ALLOCATE( ji_fix(jn_fix_pts), jj_fix(jn_fix_pts) ) + + DO jpt = 1, jn_fix_pts + READ (20, *) ji_fix(jpt), jj_fix(jpt) + ENDDO + CLOSE (20) + END IF + +! write out land/sea mask and initial bathymetry values + + zpx(:,:) = px(: ,:) ! px is only used at interior points and kiw is not used hereafter + zkiw(:,:) = kiw(:,:) ! halos could be inserted in zpx and zkiw here if not present in px and kiw + + IF ( .NOT. lsm ) zkiw(:,:) = 1.0 ! used to test whether filter is stable when all points are included in the filter + + zones(:,:) = 1.0 + PRINT *, ' point 1 zkiw' + CALL prt_summary( zkiw, zones, kpi, kpj) + + IF ( zmin_val > 0.0 ) THEN + DO jj = 2, kpj-1 + DO ji = 2,kpi-1 + IF ( zkiw(ji,jj) > 0.0 .AND. zpx(ji,jj) < zmin_val ) zpx(ji,jj) = zmin_val + ENDDO + ENDDO + ENDIF + + PRINT *, ' point 1 zpx' + CALL prt_summary( zpx, zkiw, kpi, kpj) + +! main calculations start + +! enforce cyclic conditions and north pole fold (southern boundary is assumed to be land) on zpx and zkiw + CALL impose_bcs( zpx, kpi, kpj, ll_cycl, ll_npol_fold) + CALL impose_bcs( zkiw, kpi, kpj, ll_cycl, ll_npol_fold) + + zpy (:,:) = zpx(:,:) ! initialisation of zpy is necessary for row 1 (and other outer rows/columns if non-periodic) + + zpx_iteration(:,:) = zpx(:,:) ! zpx_iteration is only updated outside the jpass and jorder loops + ! zpx is a working array updated within the jpass loop + + jiterationloop: DO jiteration=1,max_iterations + + zpx(:,:) = zpx_iteration(:,:) ! initial values for zpx for this value of jiteration + + jpassloop: DO jpass=1,kpass + + ztmp(:,:) = zpx(:,:) ! initialision of jorder loop + + CALL impose_bcs( ztmp, kpi, kpj, ll_cycl, ll_npol_fold) + + DO jorder=1,korder + + IF ( lap ) THEN + DO jj = 2,kpj-1 + DO ji = 2,kpi-1 + znum = 0.25*(ztmp(ji-1,jj )-ztmp(ji,jj))*zkiw(ji-1,jj ) & + & + 0.25*(ztmp(ji+1,jj )-ztmp(ji,jj))*zkiw(ji+1,jj ) & + & + 0.25*(ztmp(ji ,jj-1)-ztmp(ji,jj))*zkiw(ji ,jj-1) & + & + 0.25*(ztmp(ji ,jj+1)-ztmp(ji,jj))*zkiw(ji ,jj+1) - IF(ll_cycl) ihalo=1 - ! we allocate with an ihalo - ALLOCATE( ztmp(0:kpi+ihalo,kpj) , zkiw(0:kpi+ihalo,kpj) ) - ALLOCATE( zpx (0:kpi+ihalo,kpj) , zpy (0:kpi+ihalo,kpj) ) - - IF(ll_cycl) THEN - zpx(1:kpi,:) = px(: ,:) ; zkiw(1:kpi,:) = kiw(: ,:) - zpx(0 ,:) = px(kpi,:) ; zkiw(0 ,:) = kiw(kpi,:) - zpx(kpi+1,:) = px(1 ,:) ; zkiw(kpi+1,:) = kiw(1 ,:) - ELSE - zpx(: ,:) = px(: ,:) - ENDIF - - zpy (:,:) = zpx(:,:) ! init? - ztmp(:,:) = zpx(:,:) ! init - - zalphax=1./2. - zalphay=1./2. - - ! Dx/Dy=rp_aniso_diff_XY , D_ = vitesse de diffusion - ! 140 passes du fitre, Lx/Ly=1.5, le rp_aniso_diff_XY correspondant est: - - IF ( rp_aniso_diff_XY >= 1. ) zalphay=zalphay/rp_aniso_diff_XY - IF ( rp_aniso_diff_XY < 1. ) zalphax=zalphax*rp_aniso_diff_XY - - DO jorder=1,korder - imin = 2 - ihalo - imax = kpi-1 + ihalo - DO ji = imin,imax - DO jj = 2,kpj-1 - ! We crop on the coast - znum = ztmp(ji,jj) & - & + 0.25*zalphax*(ztmp(ji-1,jj )-ztmp(ji,jj))*zkiw(ji-1,jj ) & - & + 0.25*zalphax*(ztmp(ji+1,jj )-ztmp(ji,jj))*zkiw(ji+1,jj ) & - & + 0.25*zalphay*(ztmp(ji ,jj-1)-ztmp(ji,jj))*zkiw(ji ,jj-1) & - & + 0.25*zalphay*(ztmp(ji ,jj+1)-ztmp(ji,jj))*zkiw(ji ,jj+1) - zpy(ji,jj) = znum*zkiw(ji,jj)+zpx(ji,jj)*(1.-zkiw(ji,jj)) - ENDDO ! end loop ji - ENDDO ! end loop jj - - IF ( ll_cycl ) THEN - zpy(0 ,:) = zpy(kpi,:) - zpy(kpi+1,:) = zpy(1 ,:) - ENDIF - - ! update the tmp array - ztmp(:,:) = zpy(:,:) + zpy(ji,jj) = - 0.5*znum*zkiw(ji,jj) + + ENDDO ! end loop ji + ENDDO ! end loop jj + ELSE + DO jj = 1,kpj + DO ji = 2,kpi-1 + znum = 0.25*(ztmp(ji-1,jj )-ztmp(ji,jj))*zkiw(ji-1,jj ) & + & + 0.25*(ztmp(ji+1,jj )-ztmp(ji,jj))*zkiw(ji+1,jj ) + + zpy(ji,jj) = - znum*zkiw(ji,jj) - ENDDO + ENDDO ! end loop ji + ENDDO ! end loop jj + + ztmp(:,:) = zpy(:,:) + + DO jj = 2,kpj-1 + DO ji = 2,kpi-1 + znum = + 0.25*(ztmp(ji ,jj-1)-ztmp(ji,jj))*zkiw(ji ,jj-1) & + & + 0.25*(ztmp(ji ,jj+1)-ztmp(ji,jj))*zkiw(ji ,jj+1) + + zpy(ji,jj) = - znum*zkiw(ji,jj) + + ENDDO ! end loop ji + ENDDO ! end loop jj + ENDIF + + CALL impose_bcs( zpy, kpi, kpj, ll_cycl, ll_npol_fold) + + PRINT *, 'jorder, point 2 zpy = ', jorder + CALL prt_summary( zpy, zkiw, kpi, kpj) + + ztmp(:,:) = zpy(:,:) ! update ztmp for use with the next value of jorder + + ENDDO ! jorder - ! return this array - IF( ll_cycl ) THEN - py(:,:) = zpy(1:kpi,:) - ELSE - py(:,:) = zpy(: ,:) - ENDIF + zpy(:,:) = zpx(:,:) - zpy(:,:) ! zpy stores k-order filter after jpass iterations + zpx(:,:) = zpy(:,:) ! update zpx for use with the next value of jpass - END SUBROUTINE lisshapiro1d + PRINT *, 'jpass, point 3 zpx = ', jpass + CALL prt_summary( zpx, zkiw, kpi, kpj) + + END DO jpassloop + + IF ( .NOT. ( l_pass_fixed_pt_updates .OR. l_pass_shallow_updates ) ) EXIT jiterationloop + +! Find the number of points where the filtered bathymetry (zpy) is less than zmin_val (to within tolerance ztol_shallow) + + jcount_shallow = 0 + rms_int = 0.0 + PRINT *, ' zpy(ji,jj), ji, jj, jcount where zpy < zmin_val - ztol_shallow' + DO jj = 2, kpj-1 + DO ji = 2,kpi-1 + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < zmin_val - ztol_shallow ) THEN + jcount_shallow = jcount_shallow + 1 + IF ( jcount_shallow < 50 ) PRINT *, zpy(ji,jj), ji, jj, jcount_shallow + ENDIF + rms_int = rms_int + ( zpy(ji,jj) - px(ji,jj) )*( zpy(ji,jj) - px(ji,jj) ) + ENDDO + ENDDO + znpts = (kpj-2)*(kpi-2) + rms_int = SQRT( rms_int / znpts ) + PRINT *, 'jcount_shallow = ', jcount_shallow + PRINT *, 'rms_int = ', rms_int + + IF ( l_pass_fixed_pt_updates ) THEN + jcount_fixed = 0 + PRINT *, ' px(ji,jj), zpy(ji,jj), ji, jj, jcount where ABS( zpy(ji,jj) - px(ji,jj) ) > ztol_fixed ' + DO jpt = 1, jn_fix_pts + ji = ji_fix(jpt) + jj = jj_fix(jpt) + IF ( zkiw(ji,jj) > 0.0 .AND. ABS( zpy(ji,jj) - px(ji,jj) ) > ztol_fixed ) THEN ! zpy is updated value; px is original value + jcount_fixed = jcount_fixed + 1 + IF ( jcount_fixed < 50 ) PRINT *, px(ji,jj), zpy(ji,jj), ji, jj, jcount_fixed + ENDIF + ENDDO + PRINT *, 'jcount_fixed = ', jcount_fixed + ENDIF ! l_pass_fixed_pt_updates + +! If output bathymetry is too small at some sea points or not close enough to the original at the selected points, increment zpx + + l_test_shallow = jcount_shallow == 0 .OR. .NOT. l_pass_shallow_updates + l_test_fixed = jcount_fixed == 0 .OR. .NOT. l_pass_fixed_pt_updates + IF ( l_test_shallow .AND. l_test_fixed ) EXIT jiterationloop + + IF ( l_pass_shallow_updates ) THEN ! increment zpx before next pass at points where zpy < zmin_val + DO jj = 2, kpj-1 + DO ji = 2,kpi-1 + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < zmin_val) THEN + zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + MAX(px(ji,jj), zfactor_shallow*zmin_val) - zpy(ji,jj) + ENDIF + ENDDO + ENDDO + END IF + + IF ( l_pass_fixed_pt_updates ) THEN + DO jpt = 1, jn_fix_pts + ji = ji_fix(jpt) + jj = jj_fix(jpt) + IF ( zkiw(ji,jj) > 0.0 ) zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + px(ji,jj) - zpy(ji,jj) + ENDDO + END IF + + END DO jiterationloop + + IF ( l_pass_shallow_updates .AND. jcount_shallow == 0) THEN + DO jj = 1, kpj + DO ji = 1,kpi + IF ( zkiw(ji,jj) > 0.0 ) zpy(ji,jj) = MAX( zpy(ji,jj), zmin_val ) + ENDDO + ENDDO + ENDIF + + py(:,:) = zpy(: ,:) ! first use of py (halos could easily be removed here) + + DEALLOCATE( ztmp, zkiw ) + DEALLOCATE( zpx, zpy ) + DEALLOCATE( zones ) + + END SUBROUTINE lisshapiro2d SUBROUTINE lisbox(px, kiw, py, kpi, kpj, pfn, knj,panis) !!--------------------------------------------------------------------- @@ -651,4 +887,73 @@ SUBROUTINE lisbox(px, kiw, py, kpi, kpj, pfn, knj,panis) END SUBROUTINE lisbox + SUBROUTINE prt_summary( pa, pkiw, kpi, kpj) + + REAL(KIND=4), DIMENSION(:,:), INTENT(in ) :: pa, pkiw + INTEGER(KIND=4), INTENT(in ) :: kpi,kpj + + INTEGER(KIND=4) :: ji, jj + INTEGER(KIND=4) :: ji_min_prt, jj_min_prt + INTEGER(KIND=4) :: ji_max_prt, jj_max_prt + REAL(KIND=4) :: zmin, zmax + +! User sets these values + ji_min_prt = 620 ; jj_min_prt = 777 + ji_max_prt = 624 ; jj_max_prt = 781 + + + IF ( ji_max_prt > kpi ) ji_max_prt = kpi + IF ( jj_max_prt > kpj ) jj_max_prt = kpj + + DO jj = jj_min_prt, jj_max_prt + PRINT*, jj, (ji, pa(ji,jj), ji = ji_min_prt, ji_max_prt) + END DO + + zmin = 1.E20 ; zmax = -1.E20 + DO jj = 1, kpj + DO ji = 1, kpi + IF ( pkiw(ji,jj) .NE. 0 .AND. pa(ji,jj) < zmin ) THEN + ji_min_prt = ji ; jj_min_prt = jj ; zmin = pa(ji,jj) + END IF + IF ( pkiw(ji,jj) .NE. 0 .AND. pa(ji,jj) > zmax ) THEN + ji_max_prt = ji ; jj_max_prt = jj ; zmax = pa(ji,jj) + END IF + ENDDO + ENDDO + + PRINT*, 'ji, jj, min value = ', ji_min_prt, jj_min_prt, zmin + PRINT*, 'ji, jj, max value = ', ji_max_prt, jj_max_prt, zmax + + RETURN + END SUBROUTINE prt_summary + + SUBROUTINE impose_bcs( pfld, kpi, kpj, ll_cycl, ll_npol_fold) + REAL(KIND=4), DIMENSION(:,:), INTENT(inout) :: pfld + INTEGER, INTENT(in ) :: kpi, kpj + LOGICAL, INTENT(in ) :: ll_cycl, ll_npol_fold + + INTEGER ji, ijt + +! cyclic points in ji + IF ( ll_cycl ) THEN + pfld(1 ,:) = pfld(kpi-1,:) + pfld(kpi,:) = pfld(2 ,:) + ENDIF + +! north pole fold; (southern boundary is assumed to be land) + IF ( ll_npol_fold ) THEN + DO ji = 2, kpi + ijt = kpi-ji+2 + pfld(ji,kpj) = pfld(ijt,kpj-2) + END DO + pfld(1,kpj) = pfld(3,kpj-2) + DO ji = kpi/2+1, kpi + ijt = kpi-ji+2 + pfld(ji,kpj-1) = pfld(ijt,kpj-1) + END DO + END IF + + RETURN + END SUBROUTINE impose_bcs + END PROGRAM cdfsmooth From d17d2ce19b85ec9bef69cbf5067dab56880d18c8 Mon Sep 17 00:00:00 2001 From: DaveStorkey Date: Tue, 6 Aug 2024 09:54:01 +0100 Subject: [PATCH 08/15] Tidy up cdfsmooth.f90: - Add option to set name of namelist file in argument. - Move hardcoded parameters to namelist. - Put in an IF test on l_pass_shallow_updates. --- src/cdfsmooth.f90 | 74 +++++++++++++++++++++++++---------------------- 1 file changed, 39 insertions(+), 35 deletions(-) diff --git a/src/cdfsmooth.f90 b/src/cdfsmooth.f90 index 0f1a569..0bcad9c 100644 --- a/src/cdfsmooth.f90 +++ b/src/cdfsmooth.f90 @@ -85,15 +85,14 @@ PROGRAM cdfsmooth TYPE (variable), DIMENSION(:), ALLOCATABLE :: stypvar ! struture for attribute CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cv_names ! array of var name - CHARACTER(LEN=256) :: cf_in, cf_out ! file names - CHARACTER(LEN=256) :: cf_in_stem, cf_suffix ! file name parts + CHARACTER(LEN=256) :: cf_in, cf_out, cf_nam ! file names + CHARACTER(LEN=256) :: cf_in_stem, cf_suffix ! file name parts CHARACTER(LEN=256) :: cv_dep, cv_tim ! variable name for depth and time CHARACTER(LEN=256) :: ctyp ! filter type CHARACTER(LEN=256) :: cldum ! dummy character variable CHARACTER(LEN=256) :: clklist ! ciphered k-list of level LOGICAL :: lnc4 = .FALSE. ! flag for netcdf4 output with chinking and deflation - INTEGER(KIND=4) :: ji, jj, jmx, jkx ! dummy loop index !!---------------------------------------------------------------------- CALL ReadCdfNames() @@ -114,6 +113,7 @@ PROGRAM cdfsmooth PRINT *,' of iteration of the Shapiro filter.' PRINT *,' ' PRINT *,' OPTIONS :' + PRINT *,' [-n NAM-file] : namelist file. Only required for Shapiro filter.' PRINT *,' [-t FLT-type] : Lanczos , L, l (default)' PRINT *,' Hanning , H, h' PRINT *,' Shapiro , S, s' @@ -138,6 +138,7 @@ PROGRAM cdfsmooth STOP ENDIF ! + cf_nam = "" ijarg = 1 ilev = 0 ranis = 1 ! anisotropic ratio for Box car filter @@ -150,6 +151,7 @@ PROGRAM cdfsmooth CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 SELECT CASE (cldum) CASE ( '-f' ) ; CALL getarg ( ijarg, cf_in ) ; ijarg=ijarg+1 + CASE ( '-n' ) ; CALL getarg ( ijarg, cf_nam ) ; ijarg=ijarg+1 CASE ( '-c' ) ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) ncut CASE ( '-t' ) ; CALL getarg ( ijarg, ctyp ) ; ijarg=ijarg+1 CASE ( '-npass') ; CALL getarg ( ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) npass @@ -604,7 +606,7 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) LOGICAL :: l_pass_shallow_updates ! => update values where bathy < zmin_val between passes LOGICAL :: l_pass_fixed_pt_updates ! => update values where bathymetry should remain fixed - REAL(KIND=4) :: zmin_val, zfactor_shallow, ztol_fixed, ztol_shallow ! see below + REAL(KIND=4) :: rmin_val, rfactor_shallow, rtol_fixed, rtol_shallow ! see below INTEGER (KIND=4) :: ji_single_pt, jj_single_pt ! indices of single point (see l_single_point_response) INTEGER(KIND=4) :: jst_prt, jend_prt ! @@ -614,23 +616,23 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) !----------------------------------------------------------------------------------------------- NAMELIST / nam_shapiro / ll_npol_fold, ll_cycl, l_single_point_response, l_pass_shallow_updates, l_pass_fixed_pt_updates, & - & ji_single_pt, jj_single_pt, jst_prt, jend_prt + & ji_single_pt, jj_single_pt, jst_prt, jend_prt, rmin_val, rtol_shallow, rtol_fixed, rfactor_shallow !! Namelist default values ll_npol_fold = .FALSE. ll_cycl = .FALSE. l_single_point_response = .FALSE. l_pass_shallow_updates = .TRUE. l_pass_fixed_pt_updates = .TRUE. - - zmin_val = -5 ! minimum depth (e.g. 10.0 metres) - ztol_shallow = 1.0 ! tolerance in metres of minimum shallow values (zmin_val - ztol_shallow) - ztol_fixed = 1.0 ! tolerance in metres for fit to bathymetry at point specified to be fixed - zfactor_shallow = 1.5 ! zfactor_shallow needs to be slightly greater than 1.0 ! 1.1 to 1.5 are reasonable values + rmin_val = 5.0 ! minimum depth (e.g. 10.0 metres) + rtol_shallow = 1.0 ! tolerance in metres of minimum shallow values (rmin_val - rtol_shallow) + rtol_fixed = 1.0 ! tolerance in metres for fit to bathymetry at point specified to be fixed + rfactor_shallow = 1.5 ! rfactor_shallow needs to be slightly greater than 1.0 ! 1.1 to 1.5 are reasonable values ji_single_pt = 622 ; jj_single_pt = 779 jst_prt = 400 ; jend_prt = 405 - OPEN(UNIT=20, file = 'namelist_shapiro.txt', form='formatted', status='old' ) + IF( len(trim(cf_nam)) == 0 ) cf_nam='namelist_shapiro.txt' + OPEN(UNIT=20, file = trim(cf_nam), form='formatted', status='old' ) READ(NML=nam_shapiro, UNIT = 20) WRITE(NML=nam_shapiro, UNIT=6) CLOSE(20) @@ -684,10 +686,10 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) PRINT *, ' point 1 zkiw' CALL prt_summary( zkiw, zones, kpi, kpj) - IF ( zmin_val > 0.0 ) THEN + IF ( rmin_val > 0.0 ) THEN DO jj = 2, kpj-1 DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpx(ji,jj) < zmin_val ) zpx(ji,jj) = zmin_val + IF ( zkiw(ji,jj) > 0.0 .AND. zpx(ji,jj) < rmin_val ) zpx(ji,jj) = rmin_val ENDDO ENDDO ENDIF @@ -773,32 +775,34 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) IF ( .NOT. ( l_pass_fixed_pt_updates .OR. l_pass_shallow_updates ) ) EXIT jiterationloop -! Find the number of points where the filtered bathymetry (zpy) is less than zmin_val (to within tolerance ztol_shallow) +! Find the number of points where the filtered bathymetry (zpy) is less than rmin_val (to within tolerance rtol_shallow) - jcount_shallow = 0 - rms_int = 0.0 - PRINT *, ' zpy(ji,jj), ji, jj, jcount where zpy < zmin_val - ztol_shallow' - DO jj = 2, kpj-1 - DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < zmin_val - ztol_shallow ) THEN - jcount_shallow = jcount_shallow + 1 - IF ( jcount_shallow < 50 ) PRINT *, zpy(ji,jj), ji, jj, jcount_shallow - ENDIF - rms_int = rms_int + ( zpy(ji,jj) - px(ji,jj) )*( zpy(ji,jj) - px(ji,jj) ) + IF ( l_pass_shallow_updates ) THEN + jcount_shallow = 0 + rms_int = 0.0 + PRINT *, ' zpy(ji,jj), ji, jj, jcount where zpy < rmin_val - rtol_shallow' + DO jj = 2, kpj-1 + DO ji = 2,kpi-1 + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rmin_val - rtol_shallow ) THEN + jcount_shallow = jcount_shallow + 1 + IF ( jcount_shallow < 50 ) PRINT *, zpy(ji,jj), ji, jj, jcount_shallow + ENDIF + rms_int = rms_int + ( zpy(ji,jj) - px(ji,jj) )*( zpy(ji,jj) - px(ji,jj) ) + ENDDO ENDDO - ENDDO - znpts = (kpj-2)*(kpi-2) - rms_int = SQRT( rms_int / znpts ) - PRINT *, 'jcount_shallow = ', jcount_shallow - PRINT *, 'rms_int = ', rms_int + znpts = (kpj-2)*(kpi-2) + rms_int = SQRT( rms_int / znpts ) + PRINT *, 'jcount_shallow = ', jcount_shallow + PRINT *, 'rms_int = ', rms_int + ENDIF ! l_pass_shallow_updates IF ( l_pass_fixed_pt_updates ) THEN jcount_fixed = 0 - PRINT *, ' px(ji,jj), zpy(ji,jj), ji, jj, jcount where ABS( zpy(ji,jj) - px(ji,jj) ) > ztol_fixed ' + PRINT *, ' px(ji,jj), zpy(ji,jj), ji, jj, jcount where ABS( zpy(ji,jj) - px(ji,jj) ) > rtol_fixed ' DO jpt = 1, jn_fix_pts ji = ji_fix(jpt) jj = jj_fix(jpt) - IF ( zkiw(ji,jj) > 0.0 .AND. ABS( zpy(ji,jj) - px(ji,jj) ) > ztol_fixed ) THEN ! zpy is updated value; px is original value + IF ( zkiw(ji,jj) > 0.0 .AND. ABS( zpy(ji,jj) - px(ji,jj) ) > rtol_fixed ) THEN ! zpy is updated value; px is original value jcount_fixed = jcount_fixed + 1 IF ( jcount_fixed < 50 ) PRINT *, px(ji,jj), zpy(ji,jj), ji, jj, jcount_fixed ENDIF @@ -812,11 +816,11 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) l_test_fixed = jcount_fixed == 0 .OR. .NOT. l_pass_fixed_pt_updates IF ( l_test_shallow .AND. l_test_fixed ) EXIT jiterationloop - IF ( l_pass_shallow_updates ) THEN ! increment zpx before next pass at points where zpy < zmin_val + IF ( l_pass_shallow_updates ) THEN ! increment zpx before next pass at points where zpy < rmin_val DO jj = 2, kpj-1 DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < zmin_val) THEN - zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + MAX(px(ji,jj), zfactor_shallow*zmin_val) - zpy(ji,jj) + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rmin_val) THEN + zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + MAX(px(ji,jj), rfactor_shallow*rmin_val) - zpy(ji,jj) ENDIF ENDDO ENDDO @@ -835,7 +839,7 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) IF ( l_pass_shallow_updates .AND. jcount_shallow == 0) THEN DO jj = 1, kpj DO ji = 1,kpi - IF ( zkiw(ji,jj) > 0.0 ) zpy(ji,jj) = MAX( zpy(ji,jj), zmin_val ) + IF ( zkiw(ji,jj) > 0.0 ) zpy(ji,jj) = MAX( zpy(ji,jj), rmin_val ) ENDDO ENDDO ENDIF From c3a5df136cfd6f3f4e761edaef7b8ee5c6f93781 Mon Sep 17 00:00:00 2001 From: DaveStorkey Date: Tue, 6 Aug 2024 10:22:18 +0100 Subject: [PATCH 09/15] More cdfsmooth.f90 tidying: - change default Shapiro namelist settings to FALSE. - make namelist variables DOCTOR-compliant. - set name of file with list of fixed points in the namelist. --- src/cdfsmooth.f90 | 107 ++++++++++++++++++++++++---------------------- 1 file changed, 56 insertions(+), 51 deletions(-) diff --git a/src/cdfsmooth.f90 b/src/cdfsmooth.f90 index 0bcad9c..c72eeac 100644 --- a/src/cdfsmooth.f90 +++ b/src/cdfsmooth.f90 @@ -596,37 +596,42 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) !----------------------------------------------------------------------------------------------- !! Variables that can be set through the namelist - LOGICAL :: ll_npol_fold ! north pole fold in grid? - LOGICAL :: ll_cycl ! this filter has only been tested for cyclic grids -!! for limited area models the bathymetry should be smoothed over a wider domain than that of the model; the margin should be at least korder*kpass + LOGICAL :: ln_npol_fold ! north pole fold in grid? + LOGICAL :: ln_ew_cycl ! this filter has only been tested for EW cyclic grids +!! MJB: for limited area models the bathymetry should be smoothed over a wider domain than that of the model; the margin should be at least korder*kpass !! points and preferably 2*korder*kpass points (the bathymetry in nested models needs to match that of the outer domain in the nesting zone; that !! might be done by including the nesting zone in the list of fixed points - but simpler solutions might be adequate) - LOGICAL :: l_single_point_response ! => study response to single non-zero value - LOGICAL :: l_pass_shallow_updates ! => update values where bathy < zmin_val between passes - LOGICAL :: l_pass_fixed_pt_updates ! => update values where bathymetry should remain fixed + LOGICAL :: ln_single_point_response ! => study response to single non-zero value + LOGICAL :: ln_pass_shallow_updates ! => update values where bathy < zmin_val between passes + LOGICAL :: ln_pass_fixed_pt_updates ! => update values where bathymetry should remain fixed - REAL(KIND=4) :: rmin_val, rfactor_shallow, rtol_fixed, rtol_shallow ! see below + CHARACTER(LEN=256) :: cn_fixed_points ! file containing list of fixed points + ! for ln_pass_fixed_pt_updates = T - INTEGER (KIND=4) :: ji_single_pt, jj_single_pt ! indices of single point (see l_single_point_response) - INTEGER(KIND=4) :: jst_prt, jend_prt ! + REAL(KIND=4) :: rn_min_val, rn_factor_shallow, rn_tol_fixed, rn_tol_shallow ! see below + + INTEGER (KIND=4) :: ji_single_pt, jj_single_pt ! indices of single point (see ln_single_point_response) + INTEGER(KIND=4) :: jst_prt, jend_prt ! !! The points that are printed out are set by ji_min_prt, jj_min_prt, ji_max_prt, jj_max_prt - these are local to prt_summary !----------------------------------------------------------------------------------------------- - NAMELIST / nam_shapiro / ll_npol_fold, ll_cycl, l_single_point_response, l_pass_shallow_updates, l_pass_fixed_pt_updates, & - & ji_single_pt, jj_single_pt, jst_prt, jend_prt, rmin_val, rtol_shallow, rtol_fixed, rfactor_shallow + NAMELIST / nam_shapiro / ln_npol_fold, ln_ew_cycl, ln_single_point_response, ln_pass_shallow_updates, & + & ln_pass_fixed_pt_updates, cn_fixed_points, ji_single_pt, jj_single_pt, jst_prt, jend_prt, & + & rn_min_val, rn_tol_shallow, rn_tol_fixed, rn_factor_shallow !! Namelist default values - ll_npol_fold = .FALSE. - ll_cycl = .FALSE. - l_single_point_response = .FALSE. - l_pass_shallow_updates = .TRUE. - l_pass_fixed_pt_updates = .TRUE. - rmin_val = 5.0 ! minimum depth (e.g. 10.0 metres) - rtol_shallow = 1.0 ! tolerance in metres of minimum shallow values (rmin_val - rtol_shallow) - rtol_fixed = 1.0 ! tolerance in metres for fit to bathymetry at point specified to be fixed - rfactor_shallow = 1.5 ! rfactor_shallow needs to be slightly greater than 1.0 ! 1.1 to 1.5 are reasonable values + ln_npol_fold = .FALSE. + ln_ew_cycl = .FALSE. + ln_single_point_response = .FALSE. + ln_pass_shallow_updates = .FALSE. + ln_pass_fixed_pt_updates = .FALSE. + cn_fixed_points = "" + rn_min_val = 5.0 ! minimum depth (e.g. 10.0 metres) + rn_tol_shallow = 1.0 ! tolerance in metres of minimum shallow values (rn_min_val - rn_tol_shallow) + rn_tol_fixed = 1.0 ! tolerance in metres for fit to bathymetry at point specified to be fixed + rn_factor_shallow = 1.5 ! rn_factor_shallow needs to be slightly greater than 1.0 ! 1.1 to 1.5 are reasonable values ji_single_pt = 622 ; jj_single_pt = 779 jst_prt = 400 ; jend_prt = 405 @@ -657,14 +662,14 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) PRINT *, ' col kpi = ', (px(kpi,jj), jj = jst_prt, jend_prt) ! option for testing out response to a single non-zero value - IF ( l_single_point_response ) THEN + IF ( ln_single_point_response ) THEN zpx(:,:) = 0.0 zpx(ji_single_pt,jj_single_pt) = 1.0 ! choose a point somewhere in the domain END IF ! read the indices of points to remain fixed - IF ( l_pass_fixed_pt_updates ) THEN - OPEN (unit=20, file = 'Notes/list_fixed_points.txt', form='formatted', status='old' ) + IF ( ln_pass_fixed_pt_updates ) THEN + OPEN (unit=20, file = trim(cn_fixed_points), form='formatted', status='old' ) READ (20, *) jn_fix_pts ALLOCATE( ji_fix(jn_fix_pts), jj_fix(jn_fix_pts) ) @@ -686,10 +691,10 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) PRINT *, ' point 1 zkiw' CALL prt_summary( zkiw, zones, kpi, kpj) - IF ( rmin_val > 0.0 ) THEN + IF ( rn_min_val > 0.0 ) THEN DO jj = 2, kpj-1 DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpx(ji,jj) < rmin_val ) zpx(ji,jj) = rmin_val + IF ( zkiw(ji,jj) > 0.0 .AND. zpx(ji,jj) < rn_min_val ) zpx(ji,jj) = rn_min_val ENDDO ENDDO ENDIF @@ -700,8 +705,8 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) ! main calculations start ! enforce cyclic conditions and north pole fold (southern boundary is assumed to be land) on zpx and zkiw - CALL impose_bcs( zpx, kpi, kpj, ll_cycl, ll_npol_fold) - CALL impose_bcs( zkiw, kpi, kpj, ll_cycl, ll_npol_fold) + CALL impose_bcs( zpx, kpi, kpj, ln_ew_cycl, ln_npol_fold) + CALL impose_bcs( zkiw, kpi, kpj, ln_ew_cycl, ln_npol_fold) zpy (:,:) = zpx(:,:) ! initialisation of zpy is necessary for row 1 (and other outer rows/columns if non-periodic) @@ -716,7 +721,7 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) ztmp(:,:) = zpx(:,:) ! initialision of jorder loop - CALL impose_bcs( ztmp, kpi, kpj, ll_cycl, ll_npol_fold) + CALL impose_bcs( ztmp, kpi, kpj, ln_ew_cycl, ln_npol_fold) DO jorder=1,korder @@ -756,7 +761,7 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) ENDDO ! end loop jj ENDIF - CALL impose_bcs( zpy, kpi, kpj, ll_cycl, ll_npol_fold) + CALL impose_bcs( zpy, kpi, kpj, ln_ew_cycl, ln_npol_fold) PRINT *, 'jorder, point 2 zpy = ', jorder CALL prt_summary( zpy, zkiw, kpi, kpj) @@ -773,17 +778,17 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) END DO jpassloop - IF ( .NOT. ( l_pass_fixed_pt_updates .OR. l_pass_shallow_updates ) ) EXIT jiterationloop + IF ( .NOT. ( ln_pass_fixed_pt_updates .OR. ln_pass_shallow_updates ) ) EXIT jiterationloop -! Find the number of points where the filtered bathymetry (zpy) is less than rmin_val (to within tolerance rtol_shallow) +! Find the number of points where the filtered bathymetry (zpy) is less than rn_min_val (to within tolerance rn_tol_shallow) - IF ( l_pass_shallow_updates ) THEN + IF ( ln_pass_shallow_updates ) THEN jcount_shallow = 0 rms_int = 0.0 - PRINT *, ' zpy(ji,jj), ji, jj, jcount where zpy < rmin_val - rtol_shallow' + PRINT *, ' zpy(ji,jj), ji, jj, jcount where zpy < rn_min_val - rn_tol_shallow' DO jj = 2, kpj-1 DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rmin_val - rtol_shallow ) THEN + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rn_min_val - rn_tol_shallow ) THEN jcount_shallow = jcount_shallow + 1 IF ( jcount_shallow < 50 ) PRINT *, zpy(ji,jj), ji, jj, jcount_shallow ENDIF @@ -794,39 +799,39 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) rms_int = SQRT( rms_int / znpts ) PRINT *, 'jcount_shallow = ', jcount_shallow PRINT *, 'rms_int = ', rms_int - ENDIF ! l_pass_shallow_updates + ENDIF ! ln_pass_shallow_updates - IF ( l_pass_fixed_pt_updates ) THEN + IF ( ln_pass_fixed_pt_updates ) THEN jcount_fixed = 0 - PRINT *, ' px(ji,jj), zpy(ji,jj), ji, jj, jcount where ABS( zpy(ji,jj) - px(ji,jj) ) > rtol_fixed ' + PRINT *, ' px(ji,jj), zpy(ji,jj), ji, jj, jcount where ABS( zpy(ji,jj) - px(ji,jj) ) > rn_tol_fixed ' DO jpt = 1, jn_fix_pts ji = ji_fix(jpt) jj = jj_fix(jpt) - IF ( zkiw(ji,jj) > 0.0 .AND. ABS( zpy(ji,jj) - px(ji,jj) ) > rtol_fixed ) THEN ! zpy is updated value; px is original value + IF ( zkiw(ji,jj) > 0.0 .AND. ABS( zpy(ji,jj) - px(ji,jj) ) > rn_tol_fixed ) THEN ! zpy is updated value; px is original value jcount_fixed = jcount_fixed + 1 IF ( jcount_fixed < 50 ) PRINT *, px(ji,jj), zpy(ji,jj), ji, jj, jcount_fixed ENDIF ENDDO PRINT *, 'jcount_fixed = ', jcount_fixed - ENDIF ! l_pass_fixed_pt_updates + ENDIF ! ln_pass_fixed_pt_updates ! If output bathymetry is too small at some sea points or not close enough to the original at the selected points, increment zpx - l_test_shallow = jcount_shallow == 0 .OR. .NOT. l_pass_shallow_updates - l_test_fixed = jcount_fixed == 0 .OR. .NOT. l_pass_fixed_pt_updates + l_test_shallow = jcount_shallow == 0 .OR. .NOT. ln_pass_shallow_updates + l_test_fixed = jcount_fixed == 0 .OR. .NOT. ln_pass_fixed_pt_updates IF ( l_test_shallow .AND. l_test_fixed ) EXIT jiterationloop - IF ( l_pass_shallow_updates ) THEN ! increment zpx before next pass at points where zpy < rmin_val + IF ( ln_pass_shallow_updates ) THEN ! increment zpx before next pass at points where zpy < rn_min_val DO jj = 2, kpj-1 DO ji = 2,kpi-1 - IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rmin_val) THEN - zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + MAX(px(ji,jj), rfactor_shallow*rmin_val) - zpy(ji,jj) + IF ( zkiw(ji,jj) > 0.0 .AND. zpy(ji,jj) < rn_min_val) THEN + zpx_iteration(ji,jj) = zpx_iteration(ji,jj) + MAX(px(ji,jj), rn_factor_shallow*rn_min_val) - zpy(ji,jj) ENDIF ENDDO ENDDO END IF - IF ( l_pass_fixed_pt_updates ) THEN + IF ( ln_pass_fixed_pt_updates ) THEN DO jpt = 1, jn_fix_pts ji = ji_fix(jpt) jj = jj_fix(jpt) @@ -836,10 +841,10 @@ SUBROUTINE lisshapiro2d(px, kiw, py, korder, kpass, lap, lsm, kpi, kpj) END DO jiterationloop - IF ( l_pass_shallow_updates .AND. jcount_shallow == 0) THEN + IF ( ln_pass_shallow_updates .AND. jcount_shallow == 0) THEN DO jj = 1, kpj DO ji = 1,kpi - IF ( zkiw(ji,jj) > 0.0 ) zpy(ji,jj) = MAX( zpy(ji,jj), rmin_val ) + IF ( zkiw(ji,jj) > 0.0 ) zpy(ji,jj) = MAX( zpy(ji,jj), rn_min_val ) ENDDO ENDDO ENDIF @@ -931,21 +936,21 @@ SUBROUTINE prt_summary( pa, pkiw, kpi, kpj) RETURN END SUBROUTINE prt_summary - SUBROUTINE impose_bcs( pfld, kpi, kpj, ll_cycl, ll_npol_fold) + SUBROUTINE impose_bcs( pfld, kpi, kpj, ln_ew_cycl, ln_npol_fold) REAL(KIND=4), DIMENSION(:,:), INTENT(inout) :: pfld INTEGER, INTENT(in ) :: kpi, kpj - LOGICAL, INTENT(in ) :: ll_cycl, ll_npol_fold + LOGICAL, INTENT(in ) :: ln_ew_cycl, ln_npol_fold INTEGER ji, ijt ! cyclic points in ji - IF ( ll_cycl ) THEN + IF ( ln_ew_cycl ) THEN pfld(1 ,:) = pfld(kpi-1,:) pfld(kpi,:) = pfld(2 ,:) ENDIF ! north pole fold; (southern boundary is assumed to be land) - IF ( ll_npol_fold ) THEN + IF ( ln_npol_fold ) THEN DO ji = 2, kpi ijt = kpi-ji+2 pfld(ji,kpj) = pfld(ijt,kpj-2) From 136b95b6ac4fefad6973a4c15ca953dcae65dfc6 Mon Sep 17 00:00:00 2001 From: DaveStorkey Date: Tue, 6 Aug 2024 13:33:03 +0100 Subject: [PATCH 10/15] cdfsmooth.f90 : remove Mike's header comment for the module and put in a history comment instead. --- src/cdfsmooth.f90 | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/cdfsmooth.f90 b/src/cdfsmooth.f90 index c72eeac..4a4c629 100644 --- a/src/cdfsmooth.f90 +++ b/src/cdfsmooth.f90 @@ -17,15 +17,9 @@ PROGRAM cdfsmooth !! 3.0 : 01/2011 : J.M. Molines : Doctor norm + Lic. !! 3.0 : 07/2011 : R. Dussin : Add anisotropic box !! : 4.0 : 03/2017 : J.M. Molines + !! 2024 : M.J. Bell : Shapiro filter routine rewritten !!---------------------------------------------------------------------- - !! *** ROUTINE lisshapiro2d *** Mike Bell 17 Aug 2021 - !! - !! ** Purpose : apply a korder 2D shapiro filter kpass times. The land/sea mask and values at - !! The land/sea mask and values at selected points may be forced iteratively toward their initial values. - !! - !! - !!---------------------------------------------------------------------- - !! routines : description + !! routines : description !! filterinit : initialise weight !! filter : main routine for filter computation !! initlanc : initialise lanczos weights From f67f5bf2051dfaf336bfa5b217624c3bc61b76fc Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 11 Feb 2025 13:50:40 +0000 Subject: [PATCH 11/15] adding compilation file for gfortran (#5) --- Macrolib/macro.ukmo_gfortran | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 Macrolib/macro.ukmo_gfortran diff --git a/Macrolib/macro.ukmo_gfortran b/Macrolib/macro.ukmo_gfortran new file mode 100644 index 0000000..ab3881f --- /dev/null +++ b/Macrolib/macro.ukmo_gfortran @@ -0,0 +1,15 @@ +# module load netcdf-fortran/4.6.1-gcc-12.2.0-43finqs gcc/13.2.0-gcc-12.2.0-lx4jx7u + +NCDF = -I${NETCDFF_ROOT}/include -L${NETCDFF_ROOT}/lib -lnetcdff + +NC4 = -D key_netcdf4 +#CMIP6 = -D key_CMIP6 +CMIP6 = + +F90=gfortran +#OMP=-fopenmp +OMP= +FFLAGS= -O $(NCDF) $(NC4) -fno-second-underscore -ffree-line-length-256 $(OMP) + +#INSTALL=$(HOME)/local/bin +#INSTALL_MAN=$(HOME)/local/man From f0f4497493e5adcf15ab3da4c95cadf865aa5c1c Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Wed, 12 Feb 2025 11:33:40 +0000 Subject: [PATCH 12/15] bugfix to address issue #6 (#7) --- src/cdfio.F90 | 14 +++++++------- src/cdfsigtrp.f90 | 14 +++++++------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/cdfio.F90 b/src/cdfio.F90 index c421ac2..bb4ae3e 100644 --- a/src/cdfio.F90 +++ b/src/cdfio.F90 @@ -1981,13 +1981,13 @@ FUNCTION getvarxz (cdfile, cdvar, kj, kpi, kpz, kimin, kkmin, ktime) lao=.FALSE. clvar=cdvar - IF ( clvar == cn_ve3v) THEN - SELECT CASE ( cg_zgr_ver ) - CASE ( 'v2.0' ) ; clvar = 'e3v_ps' - CASE ( 'v3.0' ) ; clvar = 'e3v' - CASE ( 'v3.6' ) ; clvar = 'e3v_0' - END SELECT - ENDIF + !IF ( clvar == cn_ve3v) THEN + !SELECT CASE ( cg_zgr_ver ) + !CASE ( 'v2.0' ) ; clvar = 'e3v_ps' + !CASE ( 'v3.0' ) ; clvar = 'e3v' + !CASE ( 'v3.6' ) ; clvar = 'e3v_0' + !END SELECT + !ENDIF CALL ERR_HDL(NF90_OPEN(cdfile,NF90_NOWRITE,incid) ) CALL ERR_HDL(NF90_INQ_VARID ( incid,clvar,id_var)) diff --git a/src/cdfsigtrp.f90 b/src/cdfsigtrp.f90 index 31dd2c5..ed46621 100644 --- a/src/cdfsigtrp.f90 +++ b/src/cdfsigtrp.f90 @@ -285,11 +285,10 @@ PROGRAM cdfsigtrp lchk = lchk .OR. chkfile( cf_vfil ) IF ( lchk ) STOP 99 ! missing file - ! Look for missing value for salinity, U and V - zsps = getspval(cf_sfil, cn_vosaline ) - zspu = getspval(cf_ufil, cn_vozocrtx ) - zspv = getspval(cf_vfil, cn_vomecrty ) - + ! Look for missing value for salinity, U and V + zsps = getspval(cf_sfil, cn_vosaline ) + zspu = getspval(cf_ufil, cn_vozocrtx ) + zspv = getspval(cf_vfil, cn_vomecrty ) IF ( lg_vvl ) THEN cn_fe3u = cf_ufil @@ -327,8 +326,9 @@ PROGRAM cdfsigtrp ! Initialise sections from file ! first call to get nsection and allocate arrays IF ( lbrk ) THEN - npiglo = getdim (cf_brk, cn_x) - nsection = 1 ; iimina=1 ; iimaxa=npiglo ; ijmina=1 ; ijmaxa=1 + !npiglo = getdim (cf_brk, cn_x) + !nsection = 1 ; iimina=1 ; iimaxa=npiglo ; ijmina=1 ; ijmaxa=1 + nsection = 1 ELSE nsection = 0 CALL section_init(cf_section, csection,cvarname,clongname,iimina, iimaxa, ijmina, ijmaxa, nsection) From 41c71128d21a585ae0011672d510a92c89031dd3 Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Mon, 24 Feb 2025 13:06:16 +0000 Subject: [PATCH 13/15] adding template for nam_cdf_names (#11) --- nam_cdf_names | 104 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 nam_cdf_names diff --git a/nam_cdf_names b/nam_cdf_names new file mode 100644 index 0000000..cdcdf91 --- /dev/null +++ b/nam_cdf_names @@ -0,0 +1,104 @@ +! ================================ +! cdftools namelist +! ================================ + +!--------------------------------------- +&namdim ! Dimension names +!--------------------------------------- + cn_x = 'x' !: longitude + cn_y = 'y' !: latitude + cn_z = 'depth' !: depth + cn_t = 'time_counter' !: time +/ + +!--------------------------------------- +&namdimvar ! Dimension variable names +!--------------------------------------- + cn_vlon2d = 'nav_lon' !: longitude + cn_vlat2d = 'nav_lat' !: latitude + + cn_vdeptht = 'deptht' + cn_vdepthu = 'depthu' + cn_vdepthv = 'depthv' + cn_vdepthw = 'depthw' !: depth + + cn_vtimec = 'time_centered' !: time +/ + +!--------------------------------------- +&nammetrics ! Mesh mask variable names +!--------------------------------------- + cn_ve1t = 'e1t' + cn_ve1u = 'e1u' + cn_ve1v = 'e1v' + cn_ve1f = 'e1f' !: e1 + + cn_ve2t = 'e2t' + cn_ve2u = 'e2u' + cn_ve2v = 'e2v' + cn_ve2f = 'e2f' !: e2 + + cn_ve3t1d = 'e3t' !: e3 (1D) + cn_ve3w1d = 'e3w' !: e3 (1D) + cn_ve3t = 'e3t_0' !: e3 (3D) + cn_ve3u = 'e3u_0' !: e3 (3D) + cn_ve3v = 'e3v_0' !: e3 (3D) + cn_ve3w = 'e3w_0' !: e3 (3D) + + cn_vff = 'ff_f' !: f^2 + + cn_glamt = 'glamt' + cn_glamu = 'glamu' + cn_glamv = 'glamv' + cn_glamf = 'glamf' !: longitude + + cn_gphit = 'gphit' + cn_gphiu = 'gphiu' + cn_gphiv = 'gphiv' + cn_gphif = 'gphif' !: latitude + + cn_gdept = 'gdept' + cn_gdepw = 'gdepw' !: 1d dep variable + + cn_hdept = 'hdept' + cn_hdepw = 'hdepw' !: 2d dep variable +/ + +!--------------------------------------- +&namvars ! Variable names +!--------------------------------------- + cn_votemper = 'thetao_con' !: conservative temperature + cn_vosaline = 'so_abs' !: absolute salinity + cn_vozocrtx = 'uo' !: zonal velocity + cn_vomecrty = 'vo' !: meridional velocity + cn_sozotaux = 'tauuo' !: zonal wind stress +/ + +!--------------------------------------- +&nambathy ! Bathymetry file and variable names +!--------------------------------------- + cn_bathymet = 'bathy_metry' +/ + +!--------------------------------------- +&namsqdvar ! Squared variable names +!--------------------------------------- +/ + +!--------------------------------------- +&namcubvar ! Cubed variable names +!--------------------------------------- +/ + +!--------------------------------------- +&nammeshmask ! Mesh mask file names +!--------------------------------------- + cn_fzgr = 'mesh_mask.nc' !: vertical mesh file + cn_fhgr = 'mesh_mask.nc' !: horizontal mesh file + cn_fmsk = 'mesh_mask.nc' !: mesh mask file +/ + +!--------------------------------------- +&nammeshmask_var ! Basin mask names +!--------------------------------------- +/ From fe7e4d71937577a91be8acab6b0719541332c7f8 Mon Sep 17 00:00:00 2001 From: Diego Bruciaferri Date: Tue, 25 Feb 2025 21:32:07 +0000 Subject: [PATCH 14/15] 12 add lonlat option in cdftransport (#13) * passing netcdf libs to the linker * changes to address issue #12 --- Macrolib/macro.ukmo_gfortran | 2 +- src/Makefile | 4 ++-- src/cdftransport.f90 | 34 +++++++++++++++++++++++++--------- 3 files changed, 28 insertions(+), 12 deletions(-) diff --git a/Macrolib/macro.ukmo_gfortran b/Macrolib/macro.ukmo_gfortran index ab3881f..7ac8fb1 100644 --- a/Macrolib/macro.ukmo_gfortran +++ b/Macrolib/macro.ukmo_gfortran @@ -9,7 +9,7 @@ CMIP6 = F90=gfortran #OMP=-fopenmp OMP= -FFLAGS= -O $(NCDF) $(NC4) -fno-second-underscore -ffree-line-length-256 $(OMP) +FFLAGS= -O $(NCDF) $(NC4) -Wl,-rpath=${NETCDFF_ROOT}/lib -fno-second-underscore -ffree-line-length-256 $(OMP) #INSTALL=$(HOME)/local/bin #INSTALL_MAN=$(HOME)/local/man diff --git a/src/Makefile b/src/Makefile index ef2107f..7d50be8 100644 --- a/src/Makefile +++ b/src/Makefile @@ -287,8 +287,8 @@ cdfvtrp: cdfio.o cdfvtrp.f90 cdfpsi: cdfio.o modutils.o cdfpsi.f90 $(F90) cdfpsi.f90 -o $(BINDIR)/cdfpsi cdfio.o modcdfnames.o modutils.o $(FFLAGS) -cdftransport: cdfio.o modutils.o cdftransport.f90 - $(F90) cdftransport.f90 -o $(BINDIR)/cdftransport cdfio.o modcdfnames.o modutils.o $(FFLAGS) +cdftransport: cdfio.o modutils.o cdftools.o cdftransport.f90 + $(F90) cdftransport.f90 -o $(BINDIR)/cdftransport cdfio.o modcdfnames.o modutils.o cdftools.o $(FFLAGS) cdfvFWov: cdfio.o modutils.o cdfvFWov.f90 $(F90) cdfvFWov.f90 -o $(BINDIR)/cdfvFWov cdfio.o modcdfnames.o modutils.o $(FFLAGS) diff --git a/src/cdftransport.f90 b/src/cdftransport.f90 index fefe73e..acfdf88 100644 --- a/src/cdftransport.f90 +++ b/src/cdftransport.f90 @@ -47,6 +47,7 @@ PROGRAM cdftransport !! interm_pt : choose intermediate points on a broken line. !!---------------------------------------------------------------------- USE cdfio + USE cdftools ! for cdf_findij USE modcdfnames USE modutils ! for global attribute !!---------------------------------------------------------------------- @@ -114,11 +115,14 @@ PROGRAM cdftransport INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_varout ! Netcdf output INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipkc, id_varoutc ! Netcdf output - REAL(KIND=4) :: rxi0, ryj0 ! working variables - REAL(KIND=4) :: rxi1, ryj1 ! working variables - REAL(KIND=4) :: ai, bi ! equation of line (y=ai.x +bi) - REAL(KIND=4) :: aj, bj ! equation of line (x=aj.y +bj - REAL(KIND=4) :: rd, rd1, rd2 ! distance between point, between vertical layers + REAL(KIND=4) :: lon_min, lon_max ! longitude-limit of the section + REAL(KIND=4) :: lat_min, lat_max ! latitude-limit of the section + REAL(KIND=4) :: rxi0, ryj0 ! working variables + REAL(KIND=4) :: rxi1, ryj1 ! working variables + REAL(KIND=4) :: ai, bi ! equation of line (y=ai.x +bi) + REAL(KIND=4) :: aj, bj ! equation of line (x=aj.y +bj + REAL(KIND=4) :: rd, rd1, rd2 ! distance between point, between + ! vertical layers REAL(KIND=4) :: udum, vdum ! dummy velocity components for tests REAL(KIND=4) :: rau0=1000 ! density of pure water (kg/m3) REAL(KIND=4) :: rcp=4000. ! heat capacity (J/kg/K) @@ -218,6 +222,7 @@ PROGRAM cdftransport LOGICAL :: ltest = .FALSE. ! flag for test case LOGICAL :: lfull = .FALSE. ! flag for full step case + LOGICAL :: l_lonlat = .FALSE. ! flag for input lon/lat instead of i,j LOGICAL :: lheat = .TRUE. ! flag for skipping heat/salt transport computation LOGICAL :: lchk = .FALSE. ! flag for missing files LOGICAL :: lpm = .FALSE. ! flag for plus/minus transport @@ -240,7 +245,7 @@ PROGRAM cdftransport PRINT *,' ... [-vt VT-file] [-vtvar VT-var VS-var ]' PRINT *,' ... [-ut UT-file] [-utvar UT-var US-var ]' PRINT *,' ... [-test u v ] [-noheat ] [-pm ] [-obc] [-TS] ... ' - PRINT *,' ... [-full] [-time jt] [-vvl] [-self] ...' + PRINT *,' ... [-full] [-time jt] [-lonlat] [-vvl] [-self] ...' PRINT *,' ... [-zlimit dep_list] [-sfx suffix][-cumul]...' PRINT *,' ... [-nan]' PRINT *,' ' @@ -288,6 +293,7 @@ PROGRAM cdftransport PRINT *,' files. In this case use -t T-file option.' PRINT *,' [-full ] : use for full step configurations.' PRINT *,' [-time jt ]: compute transports for time index jt. Default is 1.' + PRINT *,' [-lonlat ]: Input section limits in longitude and latitude instead of (i,j).' PRINT *,' [-vvl ] : use time varying vertical metrics e3 read in the data file' PRINT *,' [-zlimit dep_list] : Specify a list of depth (meters) defining the ' PRINT *,' limits of classes for which transports will be computed.' @@ -353,6 +359,7 @@ PROGRAM cdftransport CASE ('-full' ) ; lfull = .TRUE. CASE ('-noheat') ; lheat = .FALSE. CASE ('-time' ) ; CALL getarg (ijarg, cldum ) ; ijarg = ijarg + 1 ; READ(cldum,*) itime + CASE ('-lonlat') ; l_lonlat = .TRUE. CASE ('-pm' ) ; lpm = .TRUE. ; lheat = .FALSE. CASE ('-obc' ) ; lobc = .TRUE. @@ -781,9 +788,18 @@ PROGRAM cdftransport dtim = getvar1d (cf_ufil, cn_vtimec, npt ) ierr = putvar1d (ncout, dtim(itime:itime), 1, 'T') - - PRINT *, ' Give iimin, iimax, ijmin, ijmax ' - READ(*,*) iimin, iimax, ijmin, ijmax + IF( l_lonlat ) THEN + PRINT *, ' Give lon_min, lon_max, lat_min, lat_max ' + READ(*,*) lon_min, lon_max, lat_min, lat_max + !! use cdf_findij to convert to model indices + CALL cdf_findij ( lon_min, lon_min, lat_min, lat_min, iimin, iimin, ijmin, ijmin, cd_coord=cn_fhgr, & + & cd_point="F", cd_verbose='n') + CALL cdf_findij ( lon_max, lon_max, lat_max, lat_max, iimax, iimax, ijmax, ijmax, cd_coord=cn_fhgr, & + & cd_point="F", cd_verbose='n') + ELSE + PRINT *, ' Give iimin, iimax, ijmin, ijmax ' + READ(*,*) iimin, iimax, ijmin, ijmax + ENDIF !! Find the broken line between P1 (iimin,ijmin) and P2 (iimax, ijmax) ! ... Initialization ii0 = iimin ; ij0 = ijmin ; ii1 = iimax ; ij1 = ijmax From 12fd66c04e48ce67baa3ab9d605862d1da6479c3 Mon Sep 17 00:00:00 2001 From: Dave Storkey Date: Mon, 17 Mar 2025 15:35:56 +0000 Subject: [PATCH 15/15] Fix for bug in cdfpsi.f90, issue #15. --- src/cdfio.F90 | 14 +++++++++++--- src/cdfpsi.f90 | 6 ++++-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/cdfio.F90 b/src/cdfio.F90 index bb4ae3e..9f71657 100644 --- a/src/cdfio.F90 +++ b/src/cdfio.F90 @@ -1422,7 +1422,7 @@ FUNCTION FindVarName(cdfile,cdvar,ld_verbose) END FUNCTION FindVarName - FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom) + FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom, ld_zeromask) !!--------------------------------------------------------------------- !! *** FUNCTION getvar *** !! @@ -1443,6 +1443,7 @@ FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom) INTEGER(KIND=4), OPTIONAL, INTENT(in) :: kimin, kjmin ! Optional variable. If missing 1 is assumed INTEGER(KIND=4), OPTIONAL, INTENT(in) :: ktime ! Optional variable. If missing 1 is assumed LOGICAL, OPTIONAL, INTENT(in) :: ldiom ! Optional variable. If missing false is assumed + LOGICAL, OPTIONAL, INTENT(in) :: ld_zeromask ! Optional variable. Reset field to zero at missing value points. If missing false is assumed REAL(KIND=4), DIMENSION(kpi,kpj) :: getvar ! 2D REAL 4 holding variable field at klev INTEGER(KIND=4), DIMENSION(4) :: istart, icount, inldim @@ -1454,7 +1455,7 @@ FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom) REAL(KIND=4) :: spval !: missing value REAL(KIND=4) , DIMENSION (:,:), ALLOCATABLE :: zend, zstart CHARACTER(LEN=256) :: clvar - LOGICAL :: lliom=.false., llperio=.false. + LOGICAL :: lliom=.false., llperio=.false., ll_zeromask=.false. LOGICAL :: llog=.FALSE. , lsf=.FALSE. , lao=.FALSE. !! INTEGER(KIND=4) :: ityp @@ -1499,6 +1500,12 @@ FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom) lliom=.false. ENDIF + IF (PRESENT(ld_zeromask) ) THEN + ll_zeromask=ld_zeromask + ELSE + ll_zeromask=.false. + ENDIF + ! Must reset the flags to false for every call to getvar clvar=cdvar llog = .FALSE. @@ -1603,7 +1610,8 @@ FUNCTION getvar (cdfile,cdvar,klev,kpi,kpj,kimin,kjmin, ktime, ldiom) IF (lsf ) WHERE (getvar /= spval ) getvar=getvar*sf IF (lao ) WHERE (getvar /= spval ) getvar=getvar + ao IF (llog) WHERE (getvar /= spval ) getvar=10**getvar - + IF (ll_zeromask) WHERE (getvar == spval ) getvar=0.0 + istatus=NF90_CLOSE(incid) END FUNCTION getvar diff --git a/src/cdfpsi.f90 b/src/cdfpsi.f90 index 8478b5c..809b294 100644 --- a/src/cdfpsi.f90 +++ b/src/cdfpsi.f90 @@ -243,6 +243,7 @@ PROGRAM cdfpsi CALL CreateOutput + PRINT *, 'Getting ',TRIM(cn_ve1v), TRIM(cn_ve2u) e1v(:,:) = getvar(cn_fhgr, cn_ve1v, 1, npiglo, npjglo) e2u(:,:) = getvar(cn_fhgr, cn_ve2u, 1, npiglo, npjglo) IF ( lmask) THEN @@ -270,7 +271,7 @@ PROGRAM cdfpsi DO jk = 1,npk IF ( ll_v ) THEN - zv(:,:) = getvar(cf_vfil, cn_vomecrty, jk, npiglo, npjglo, ktime=jt ) + zv(:,:) = getvar(cf_vfil, cn_vomecrty, jk, npiglo, npjglo, ktime=jt, ld_zeromask=.true. ) IF ( lfull ) THEN ; e3v(:,:) = e31d(jk) ELSE ; e3v(:,:) = getvar(cn_fe3v, cn_ve3v, jk, npiglo, npjglo, ktime=it, ldiom=.NOT.lg_vvl) ENDIF @@ -289,7 +290,8 @@ PROGRAM cdfpsi ENDIF IF ( ll_u) THEN - zu(:,:) = getvar(cf_ufil, cn_vozocrtx, jk, npiglo, npjglo, ktime=jt ) + zu(:,:) = getvar(cf_ufil, cn_vozocrtx, jk, npiglo, npjglo, ktime=jt, ld_zeromask=.true. ) + PRINT *, 'Getting ',TRIM(cn_ve3u) IF ( lfull ) THEN ; e3u(:,:) = e31d(jk) ELSE ; e3u(:,:) = getvar(cn_fe3u, cn_ve3u, jk, npiglo, npjglo, ktime=it, ldiom=.NOT.lg_vvl) ENDIF