diff --git a/Makefile.global b/Makefile.global index 795a2e8..726f075 100644 --- a/Makefile.global +++ b/Makefile.global @@ -16,55 +16,48 @@ ARG_FLAGS = -DARG_INT1=$(INPUT_INT1) -DARG_INT2=$(INPUT_INT2) OBJDIR := build ifeq ($(strip $(FF)),) - FF = $(F95COMPILER) + FF = $(F90COMPILER) endif ifeq ($(NETCDFLIBS),"none") LIB_DIR = INC_DIR = ORM_FLAGS += -Dno_netcdf -endif -ifeq ($(NETCDFLIBS),"automatic") +else ifeq ($(NETCDFLIBS),"automatic") LIB_DIR = $(shell nc-config --flibs) INC_DIR = -I$(shell nc-config --includedir) -endif -ifeq ($(NETCDFLIBS),"automatic-44") +else ifeq ($(NETCDFLIBS),"automatic-44") INC_DIR = $(shell nf-config --cflags) LIB_DIR = $(shell nf-config --flibs) -endif -ifeq ($(NETCDFLIBS),"macports") - LIB_DIR = -L/opt/local/lib - INC_DIR = -I/opt/local/include/ - LNK_FLAGS = -lnetcdf -lnetcdff -endif -ifeq ($(NETCDFLIBS),"fink") - LIB_DIR = -L/sw/lib # -L/sw/lib/netcdf-gfortran/lib - INC_DIR = -I/sw/include #-I/sw/lib/netcdf-gfortran/include - LNK_FLAGS = -lnetcdf -lnetcdff -lsz +else + LIB_DIR = -L/usr/local/lib + INC_DIR = -I/usr/local/include/ + LNK_FLAGS = -lnetcdf -lnetcdff endif -ifeq ($(F95COMPILER),"gfortran") - FF_FLAGS = -O3 -c -x f95-cpp-input -g -fbacktrace -fbounds-check -# FF_FLAGS = -O3 -c -x f95-cpp-input -g -fbacktrace -fbounds-check -ffpe-trap=denormal -static -Wall -pedantic - #FF_FLAGS = -O3 -c -x f95-cpp-input +ifeq ($(ARCH),"macports") + F90COMPILER = gfortran-mp-8 + PROD_FF_FLAGS = -O3 -c -x f95-cpp-input -g -fbacktrace -fbounds-check + DEV_FF_FLAGS = -O1 -c -x f95-cpp-input -g -fbacktrace -fbounds-check -ffpe-trap=denormal -static -Wall -pedantic F90_FLAGS = -fno-underscoring FF += $(LIB_DIR) $(INC_DIR) $(F90_FLAGS) $(ORM_FLAGS) -J$(OBJDIR) endif -ifeq ($(F95COMPILER),"g95") - FF_FLAGS = -c -cpp -fendian=big - F90_FLAGS = -O3 -C -g -fno-underscoring - FF += $(LIB_DIR) $(INC_DIR) $(F90_FLAGS) $(ORM_FLAGS) -endif -ifeq ($(F95COMPILER),"ifort") - FF_FLAGS = -O3 -cpp -convert big_endian -xcore-avx512 +ifeq ($(ARCH),"tetralith") + F90COMPILER = ifort + PROD_FF_FLAGS = -O3 -cpp -convert big_endian -xcore-avx512 + DEV_FF_FLAGS = -O1 -cpp -convert big_endian -xcore-avx512 -g -traceback -fpe0 F90_FLAGS = -free ARG_FLAGS += -Tf FF += $(LIB_DIR) $(INC_DIR) $(F90_FLAGS) $(ORM_FLAGS) TT_FLAGS = -c endif - +ifeq ($(OPT),"opt") + FF_FLAGS = $(PROD_FF_FLAGS) +else + FF_FLAGS = $(DEV_FF_FLAGS) +endif CC = gcc -O $(INC_DIR) @@ -75,12 +68,12 @@ VPATH = src:projects/$(PROJECT):src/active_particles all: runfile -objects := $(addprefix $(OBJDIR)/,modules.o interp_general.o $(ACTVEL) calendar.o \ +objects := $(addprefix $(OBJDIR)/,modules.o interp_general.o interp2.o $(ACTVEL) calendar.o \ laplacian.o savepsi.o savepsi_new.o loop_pos.o writetrajs.o \ sw_stat.o tracer_functions.o \ seed.o init_seed.o getfile.o \ vertvel.o coord.o cross.o init_par.o time_subs.o \ - pos.o interp2.o sw_seck.o sw_pres.o sw_dens0.o \ + pos.o sw_seck.o sw_pres.o sw_dens0.o \ writepsi.o writetracer.o printinfo.o loop.o main.o \ setupgrid.o readfield.o diffusion.o) diff --git a/Makefile_tmpl b/Makefile_tmpl index fb965ab..90d7d67 100644 --- a/Makefile_tmpl +++ b/Makefile_tmpl @@ -2,17 +2,23 @@ PROJECT = theoretical CASE = theoretical RUNFILE = runtracmass -#F95COMPILER = "g95" -F95COMPILER = "gfortran" -#F95COMPILER = "ifort" -#FF = "/sw/bin/gfortran-fsf-4.9" #Alt. compiler location +# Possible architectures: +# tetralith (Swedish HPC with intel) +# macports (MacPorts for Mac OS) +# generic (copy this to make your own arch) +ARCH = "macports" +#FF = "/sw/bin/gfortran-fsf-4.9" #Alt. compiler location -NETCDFLIBS = "none" -#NETCDFLIBS = "automatic" # set using nc-config -#NETCDFLIBS = "automatic-44" # set using nf-config -#NETCDFLIBS = "macports" -#NETCDFLIBS = "fink" -#NETCDFLIBS = "source" +# Possible netCDF settings +# automatic (set by nc-config) +# automatic-44 (set by nf-config, for netCDF version >4.4) +# generic (set by yourself in Makefile.global) +# none (no netCDF) +NETCDFLIBS = "automatic-44" + +# Opt = optimise, i.e. dont use so many error checks +# Else TRACMASS will compile with lots of error check which might slowdown the code +OPT = "opt" #================================================================ diff --git a/README.md b/README.md index ebdb4b2..89ab696 100644 --- a/README.md +++ b/README.md @@ -13,18 +13,42 @@ It is still in development. Quickstart ========== -Prerequisites to run on Mac OS X --------------------------------- +Download the code +----------------- -Download and install macports (www.macports.org) -Install the following ports: +```bash +git clone https://github.com/TRACMASS/tracmass.git +``` -```sh -sudo port install openmpi +gcc45 +valgrind -sudo port install netcdf +openmpi +netcdf4 -sudo port install netcdf-fortran +gcc45 +openmpi -sudo port install git-core -sudo port install git-extras +Compile +------- + +Enter the tracmass directory and copy the template Makefile + +```bash +cd tracmass +cp Makefile_tmpl Makefile +``` + +Modify the Makefile to fit your system. +You will need to set ARCH, which is the name of your system, i.e. macports, tetralith etc. +You will also need to configure how TRACMASS should find the netCDF libraries, if at all. +For most systems, we recommend the option automatic-44. +Then you can run the make command. + +```bash +make +``` + +Running a first test case +------------------------- + +We recommend testing that TRACMASS was properly compiled by letting PROJECT and CASE be "theoretical" in the Makefile (which is the default). +In this case, TRACMASS will use a simple oscillating velocity field to trace trajectories. +You can run this case by setting + +```bash +./runtracmass ``` Download test data @@ -36,26 +60,43 @@ You can find some input data for testing the code on https://www.dropbox.com/sh/b6leim7tb99i3hm/AAA7M0mxuYTwJLKjqmLnoJuca?dl=0 ``` +This test data includes data from NEMO, IFS and ROMS. Before doing any analysis we recommend to download some of the test data and make sure TRACMASS is working properly. -Compile the code ----------------- - -Copy Makefile_tmpl to Makefile. -Choose a change PROJECT and CASE in Makefile to your projectname. -If you are using netCDF version >4.4 you can set NETCDFLIBS to "automatic-44", in which case netCDF paths will be taken from nc-config and nf-config. +In order to set up TRACMASS to run trajectories on e.g. NEMO data, you will need to change PROJECT and CASE to NEMO, and then re-compile the code. -run: - -```bash +```bash +make clean make -./runtrm +./runtracmass ``` +Run your analysis +----------------- + +If you wish to run TRACMASS using fields from NEMO, ROMS, IFS or CMEMS satellite altimetry, you may use these projects as is (see projects directory). +If you wish to run another case or a very specific case of the above models, you may create your own project in the projects directory. +To run with e.g. your own NEMO data, you will need to modify the projects/nemo/nemo.in namelist to suit your needs. + Change log ========== 6.0.0 First public release -7.0.0_alpha New features for parametrising sub-grid scale motions, tracing water in the atmosphere, simplified set-up for users + code clean–up + bug-fixes. \ No newline at end of file +7.0.0_alpha New features for parametrising sub-grid scale motions, tracing water in the atmosphere, simplified set-up for users + code clean–up + bug-fixes. + + +Prerequisites to run on Mac OS X +================================ + +Download and install macports (www.macports.org) +Install the following ports: + +```sh +sudo port install openmpi +gcc45 +valgrind +sudo port install netcdf +openmpi +netcdf4 +sudo port install netcdf-fortran +gcc45 +openmpi +sudo port install git-core +sudo port install git-extras +``` diff --git a/projects/theoretical/readfield.f95 b/projects/theoretical/readfield.F90 similarity index 98% rename from projects/theoretical/readfield.f95 rename to projects/theoretical/readfield.F90 index b10a972..e5621c5 100755 --- a/projects/theoretical/readfield.f95 +++ b/projects/theoretical/readfield.F90 @@ -6,7 +6,8 @@ subroutine readfields USE mod_grid USE mod_name USE mod_vel - ! USE mod_dens + USE mod_dens + USE mod_tempsalt USE mod_stat IMPLICIT none @@ -47,11 +48,9 @@ subroutine readfields do i=1,IMT im=i-1 if(im.eq.0) im=IMT -#ifdef tempsalt tem (i,j,k,2)=20.*float(k)/float(km) sal (i,j,k,2)=30. rho (i,j,k,2)=(28.-20.)*float(km-k)/float(km) +20. -#endif ! time evolving oscilating field ! ------------------------------ diff --git a/projects/theoretical/setupgrid.f95 b/projects/theoretical/setupgrid.F90 similarity index 93% rename from projects/theoretical/setupgrid.f95 rename to projects/theoretical/setupgrid.F90 index 4909c2d..f0f4d8b 100644 --- a/projects/theoretical/setupgrid.f95 +++ b/projects/theoretical/setupgrid.F90 @@ -44,14 +44,14 @@ SUBROUTINE setupgrid integer, dimension(2) :: dim2d integer, dimension(3) :: dim3d real, allocatable, dimension(:,:) :: lon, lat - real, allocatable, dimension(:) :: depth + real, allocatable, dimension(:) :: depth_1d logical :: lwrite_nc, lread_nc lwrite_nc = .true. lread_nc = .true. allocate( lon(imt,jmt), lat(imt,jmt) ) -allocate( depth(km) ) +allocate( depth_1d(km) ) kmt=KM ! flat bottom @@ -74,9 +74,9 @@ SUBROUTINE setupgrid end do end do -depth(1) = dz(1)/2. +depth_1d(1) = dz(1)/2. do k = 2, km - depth(k) = SUM(dz(1:k-1)) + dz(k)/2. + depth_1d(k) = SUM(dz(1:k-1)) + dz(k)/2. end do if (lwrite_nc) then @@ -98,7 +98,7 @@ SUBROUTINE setupgrid call check( nf90_put_var(ncid, idlon, lon) ) call check( nf90_put_var(ncid, idlat, lat) ) - call check( nf90_put_var(ncid, iddep, depth) ) + call check( nf90_put_var(ncid, iddep, depth_1d) ) call check( nf90_put_var(ncid, iddx, dxv(1:jmt,1:imt)) ) call check( nf90_put_var(ncid, iddy, dyu(1:jmt,1:imt)) ) call check( nf90_put_var(ncid, iddz, dzt(1:imt,1:jmt,1:km,2)) ) @@ -112,13 +112,16 @@ SUBROUTINE setupgrid map3d = [2, 3, 4, 1] lon = get2DfieldNC('mesh.nc', 'lon') lat = get2DfieldNC('mesh.nc', 'lat') - depth = get1DfieldNC('mesh.nc', 'depth') + depth_1d = get1DfieldNC('mesh.nc', 'depth') dxv(1:imt,1:jmt) = get2DfieldNC('mesh.nc', 'dx') dyu(1:imt,1:jmt) = get2DfieldNC('mesh.nc', 'dy') dzt(1:imt,1:jmt,1:km,2) = get3DfieldNC('mesh.nc', 'dz') dzt(:,:,:,1) = dzt(:,:,:,2) end if +depth(:,:) = depth_1d(km) + + ! === diff --git a/projects/theoretical/theoretical.in b/projects/theoretical/theoretical.in index fd4bd0c..9855f45 100644 --- a/projects/theoretical/theoretical.in +++ b/projects/theoretical/theoretical.in @@ -37,7 +37,8 @@ / &INIT_GRID_TIME fieldsperfile = 1, - ngcm = 1, + ngcm_step = 1, + ngcm_unit = 4, iter = 1, intmax = 146000, minveljd = -1, @@ -67,7 +68,7 @@ ! 5 = write at chosen intervals ! 6 = write each spatial grid-crossing kriva = 6, - outdatadir = '/proj/bolinc/users/x_sarbe/TRACMASS/projects/theoretical/', + outdatadir = 'data_out' outdatafile = 'test_k6', ! outdatafile = 'analyt_k1', ! outdatafile = 'dt1_k1', @@ -169,3 +170,5 @@ twave = 8.0, critvel = 0.1, / +&INIT_TRACERS +/ diff --git a/src/cross.F90 b/src/cross.F90 index 12d8ed3..526035d 100755 --- a/src/cross.F90 +++ b/src/cross.F90 @@ -39,8 +39,6 @@ subroutine cross_stat(ijk,ia,ja,ka,r0,sp,sn) if(im.eq.0) im=IMT uu=(intrpg*uflux(ia,ja,ka,nsp)+intrpr*uflux(ia,ja,ka,nsm))!*ff um=(intrpg*uflux(im,ja,ka,nsp)+intrpr*uflux(im,ja,ka,nsm))!*ff - frac1 = min( abs(upr(1,1)), abs(0.99*uu/upr(1,1)) ) - frac2 = min( abs(upr(2,1)), abs(0.99*um/upr(2,1)) ) !print*,'cross x',uu,um,upr(1,1),upr(2,1) #ifdef turb if(r0.ne.dble(ii)) then @@ -66,8 +64,6 @@ subroutine cross_stat(ijk,ia,ja,ka,r0,sp,sn) ii=ja uu=(intrpg*vflux(ia,ja ,ka,nsp)+intrpr*vflux(ia,ja ,ka,nsm))!*ff um=(intrpg*vflux(ia,ja-1,ka,nsp)+intrpr*vflux(ia,ja-1,ka,nsm))!*ff - frac1 = min( abs(upr(3,1)), abs(0.99*uu/upr(3,1)) ) - frac2 = min( abs(upr(4,1)), abs(0.99*um/upr(4,1)) ) !print*,'cross y',uu,um,upr(3,1),upr(4,1) #ifdef turb if(r0.ne.dble(ja )) then diff --git a/src/init_par.F90 b/src/init_par.F90 index 2eb3aaa..697034c 100644 --- a/src/init_par.F90 +++ b/src/init_par.F90 @@ -410,14 +410,12 @@ SUBROUTINE init_params seedsubints(1) = 0 end if -#ifdef tempsalt ALLOCATE ( tem(imt,jmt,km,nst) ) ALLOCATE ( sal(imt,jmt,km,nst) ) ALLOCATE ( rho(imt,jmt,km,nst) ) tem = 0. sal = 0. rho = 0. -#endif ALLOCATE ( tracers2D(n2Dtracers), tracers3D(n3Dtracers) ) DO jt=1,n2Dtracers diff --git a/src/interp2.F90 b/src/interp2.F90 index 299b9bb..abac51e 100755 --- a/src/interp2.F90 +++ b/src/interp2.F90 @@ -1,4 +1,3 @@ -#ifdef tempsalt subroutine interp2(i,j,k,temp,salt,dens) ! === NO interpolation of the temperature, salinity, and density=== @@ -28,5 +27,5 @@ subroutine interp2(i,j,k,temp,salt,dens) return end subroutine interp2 -#endif + diff --git a/src/streamfunctions.F90 b/src/streamfunctions.F90 index 6e99244..6e19ac9 100644 --- a/src/streamfunctions.F90 +++ b/src/streamfunctions.F90 @@ -8,15 +8,19 @@ module mod_pos USE mod_streamfunctions USE mod_psi use mod_psi_new - + USE mod_interp, only: interp_gen2D, interp_gen3D + IMPLICIT none contains subroutine pos_stream_funcs_1() + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) + #if defined streamr || streamts ! call interp(ib,jb,kb,x1,y1,z1,temp,salt,dens,1) - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mrb=int((dens-rmin)/dr)+1 if(mrb.lt.1 ) mrb=1 if(mrb.gt.MR) mrb=MR @@ -31,14 +35,16 @@ subroutine pos_stream_funcs_1() #endif #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) + !call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR @@ -53,9 +59,11 @@ end subroutine pos_stream_funcs_1 subroutine pos_stream_funcs_2() + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) #if defined streamr || streamts ! call interp(ib,jb,kb,x1,y1,z1,temp,salt,dens,1) - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mrb=int((dens-rmin)/dr)+1 if(mrb.lt.1 ) mrb=1 if(mrb.gt.MR) mrb=MR @@ -70,14 +78,16 @@ subroutine pos_stream_funcs_2() #endif #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR @@ -90,9 +100,12 @@ end subroutine pos_stream_funcs_2 subroutine pos_stream_funcs_3() + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) + #if defined streamr || streamts ! call interp(ib,jb,kb,x1,y1,z1,temp,salt,dens,1) - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mrb=int((dens-rmin)/dr)+1 if(mrb.lt.1 ) mrb=1 if(mrb.gt.MR) mrb=MR @@ -107,14 +120,16 @@ subroutine pos_stream_funcs_3() #endif #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR @@ -127,9 +142,11 @@ end subroutine pos_stream_funcs_3 subroutine pos_stream_funcs_4() + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) #if defined streamr || streamts ! call interp(ib,jb,kb,x1,y1,z1,temp,salt,dens,1) - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mrb=int((dens-rmin)/dr)+1 if(mrb.lt.1 ) mrb=1 if(mrb.gt.MR) mrb=MR @@ -144,14 +161,16 @@ subroutine pos_stream_funcs_4() #endif #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR @@ -164,17 +183,20 @@ end subroutine pos_stream_funcs_4 subroutine pos_stream_funcs_5() - + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) + #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR @@ -188,17 +210,20 @@ end subroutine pos_stream_funcs_5 subroutine pos_stream_funcs_6() - + integer :: mtrb(n3Dtracers), mtra(n3Dtracers) + #if defined stream_thermohaline ! calculate the layers of temperature and salinity for both a-box and b-box - call interp2(ib,jb,kb,temp,salt,dens) +! call interp2(ib,jb,kb,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mtb=int((temp-tmin)/dtemp)+1 if(mtb.lt.1 ) mtb=1 if(mtb.gt.MR) mtb=MR msb=int((salt-smin)/dsalt)+1 if(msb.lt.1 ) msb=1 if(msb.gt.MR) msb=MR - call interp2(ia,ja,ka,temp,salt,dens) +! call interp2(ia,ja,ka,temp,salt,dens) + call interp_gen3D(ib,jb,kb,x1,y1,z1,2,trc3D,method='nearest') mta=(temp-tmin)/dtemp+1 if(mta.lt.1 ) mta=1 if(mta.gt.MR) mta=MR