diff --git a/Macrolib/macro.ukmo_gfortran b/Macrolib/macro.ukmo_gfortran
new file mode 100644
index 0000000..7ac8fb1
--- /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) -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/README.md b/README.md
index 5f92759..8e3732c 100644
--- a/README.md
+++ b/README.md
@@ -1,16 +1,12 @@
-# CDFTOOLS
- 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).
-
+# 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 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:
+#### 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 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
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
+!---------------------------------------
+/
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/cdfio.F90 b/src/cdfio.F90
index c421ac2..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
@@ -1981,13 +1989,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/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
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)
diff --git a/src/cdfsmooth.f90 b/src/cdfsmooth.f90
index 34f7339..4a4c629 100644
--- a/src/cdfsmooth.f90
+++ b/src/cdfsmooth.f90
@@ -17,9 +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
!!----------------------------------------------------------------------
- !!----------------------------------------------------------------------
- !! routines : description
+ !! routines : description
!! filterinit : initialise weight
!! filter : main routine for filter computation
!! initlanc : initialise lanczos weights
@@ -28,7 +28,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 +53,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 +67,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
@@ -76,7 +79,8 @@ 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, 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
@@ -89,8 +93,8 @@ PROGRAM cdfsmooth
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'
@@ -103,41 +107,53 @@ 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'
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 *,' '
STOP
ENDIF
!
+ cf_nam = ""
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 ( '-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
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 +163,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 +177,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 +290,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 +336,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 +553,303 @@ 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 :: 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 :: 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
+
+ CHARACTER(LEN=256) :: cn_fixed_points ! file containing list of fixed points
+ ! for ln_pass_fixed_pt_updates = T
+
+ 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 / 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
+ 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
+
+ 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)
+
+ 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 ( 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 ( 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) )
+
+ 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 ( 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) < rn_min_val ) zpx(ji,jj) = rn_min_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, 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)
+
+ 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, ln_ew_cycl, ln_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
- ! return this array
- IF( ll_cycl ) THEN
- py(:,:) = zpy(1:kpi,:)
- ELSE
- py(:,:) = zpy(: ,:)
- ENDIF
+ 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)
- END SUBROUTINE lisshapiro1d
+ ztmp(:,:) = zpy(:,:) ! update ztmp for use with the next value of jorder
+
+ ENDDO ! jorder
+
+ zpy(:,:) = zpx(:,:) - zpy(:,:) ! zpy stores k-order filter after jpass iterations
+ zpx(:,:) = zpy(:,:) ! update zpx for use with the next value of jpass
+
+ PRINT *, 'jpass, point 3 zpx = ', jpass
+ CALL prt_summary( zpx, zkiw, kpi, kpj)
+
+ END DO jpassloop
+
+ 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 rn_min_val (to within tolerance rn_tol_shallow)
+
+ IF ( ln_pass_shallow_updates ) THEN
+ jcount_shallow = 0
+ rms_int = 0.0
+ 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) < 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
+ 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
+ ENDIF ! ln_pass_shallow_updates
+
+ 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) ) > 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) ) > 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 ! 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. 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 ( 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) < 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 ( ln_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 ( 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), rn_min_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 +890,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, ln_ew_cycl, ln_npol_fold)
+ REAL(KIND=4), DIMENSION(:,:), INTENT(inout) :: pfld
+ INTEGER, INTENT(in ) :: kpi, kpj
+ LOGICAL, INTENT(in ) :: ln_ew_cycl, ln_npol_fold
+
+ INTEGER ji, ijt
+
+! cyclic points in ji
+ 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 ( ln_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
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
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