[cdftools] 134/228: JMM: modification of transport tools ( cdftransport.f90 and cdfsigtrp.f90) in order to be able to fix variable name and longname from standard input or input file. This is an option. If not provided, those tools work as before.
Alastair McKinstry
mckinstry at moszumanska.debian.org
Fri Jun 12 08:21:40 UTC 2015
This is an automated email from the git hooks/post-receive script.
mckinstry pushed a commit to branch master
in repository cdftools.
commit 116db0b70ce8d789ad160279ddf9e2b96ff3f6e9
Author: molines <molines at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date: Sun Apr 1 12:17:13 2012 +0000
JMM: modification of transport tools ( cdftransport.f90 and cdfsigtrp.f90) in
order to be able to fix variable name and longname from standard input
or input file. This is an option. If not provided, those tools work as
before.
git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@582 1055176f-818a-41d9-83e1-73fbe5b947c5
---
cdfsigtrp.f90 | 1807 +++++++++++++++++++++---------------------
cdftransport.f90 | 2294 ++++++++++++++++++++++++++++--------------------------
2 files changed, 2107 insertions(+), 1994 deletions(-)
diff --git a/cdfsigtrp.f90 b/cdfsigtrp.f90
index 86710bd..10237e9 100644
--- a/cdfsigtrp.f90
+++ b/cdfsigtrp.f90
@@ -1,885 +1,942 @@
PROGRAM cdfsigtrp
- !!======================================================================
- !! *** PROGRAM cdfsigtrp ***
- !!=====================================================================
- !! ** Purpose : Compute density class Mass transport across a section.
- !!
- !! ** Method :- The begining and end point of the section are given in
- !! term of f-points index.
- !! - The program works for zonal or meridional sections.
- !! - The section definitions are given in an ASCII FILE
- !! dens_section.dat:
- !! foreach sections, 2 lines :
- !! (i) : section name (String, no blank)
- !! (ii) : imin imax jmin jmax for the section
- !! - Only vertical slices corrsponding to the sections are
- !! read in the files.
- !! - read metrics, depth, etc
- !! - read normal velocity (either vozocrtx oy vomecrty )
- !! - read 2 rows of T and S ( i i+1 or j j+1 )
- !! - compute the mean value at velocity point
- !! - compute sigma0 (can be easily modified for sigmai )
- !! - compute the depths of isopyncal surfaces
- !! - compute the transport from surface to the isopycn
- !! - compute the transport in each class of density
- !! - compute the total transport (for information)
- !!
- !! History : 2.1 : 03/2006 : J.M. Molines : Original code
- !! : 07/2009 : R. Dussin : add cdf output
- !! 3.0 : 06/2011 : J.M. Molines : Doctor norm + Lic.
- !!----------------------------------------------------------------------
- !!----------------------------------------------------------------------
- !! routines : description
- !! section_init : initialize section names and positions
- !! print_out : routine which performs standard output if required
- !! bimg_writ : routine which performs bimg output if required
- !!----------------------------------------------------------------------
- USE cdfio
- USE eos ! for sigma0, sigmai
- USE modcdfnames ! for ReadCdfNames
- USE modutils ! for SetGlobalAtt
- !!----------------------------------------------------------------------
- !! CDFTOOLS_3.0 , MEOM 2011
- !! $Id$
- !! Copyright (c) 2011, J.-M. Molines
- !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
- !!----------------------------------------------------------------------
- IMPLICIT NONE
-
- INTEGER(KIND=4) :: ji, jk, jclass, jsec ! dummy loop index
- INTEGER(KIND=4) :: jiso, jbin, jarg ! dummy loop index
- INTEGER(KIND=4) :: nbins ! number of density classes
- INTEGER(KIND=4) :: ipos ! working variable
- INTEGER(KIND=4) :: narg, iargc ! command line
- INTEGER(KIND=4) :: ijarg, ireq ! command line
- INTEGER(KIND=4) :: npk, nk ! vertical size, number of wet layers
- INTEGER(KIND=4) :: numbimg=10 ! optional bimg logical unit
- INTEGER(KIND=4) :: numout=11 ! ascii output
- INTEGER(KIND=4) :: nsection ! number of sections (overall)
- INTEGER(KIND=4) :: iimin, iimax ! working section limits
- INTEGER(KIND=4) :: ijmin, ijmax ! working section limits
- INTEGER(KIND=4) :: npts ! number of points in section
- INTEGER(KIND=4) :: ikx=1, iky=1 ! dims of netcdf output file
- INTEGER(KIND=4) :: nboutput=2 ! number of values to write in cdf output
- INTEGER(KIND=4) :: ncout, ierr ! for netcdf output
- INTEGER(KIND=4) :: iweight ! weight of input file for further averaging
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: iimina, iimaxa ! sections limits
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ijmina, ijmaxa ! sections limits
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_varout ! variable levels and id
-
- REAL(KIND=4) :: refdep =0.e0 ! reference depth (m)
- REAL(KIND=4), DIMENSION(1) :: tim ! time counter
- REAL(KIND=4), DIMENSION(1) :: rdummy1, rdummy2 ! working variable
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdept, gdepw ! depth of T and W points
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: eu ! either e1v or e2u
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: e3t1d, e3w1d ! vertical metrics in case of full step
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rlonlat ! longitudes/latitudes if the section
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zs, zt ! salinity and temperature from file
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rdumlon, rdumlat ! dummy longitude and latitude for output
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zu ! velocity
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zmask ! mask
- REAL(KIND=4), DIMENSION(:,:,:), ALLOCATABLE :: tmpm, tmpz ! temporary arrays
-
- ! double precision for cumulative variables and densities
- REAL(KIND=8) :: dsigma_min ! minimum density for bining
- REAL(KIND=8) :: dsigma_max, dltsig ! maximum density for bining, step
- REAL(KIND=8) :: dsigma, dalfa ! working sigma, interpolation coeff.
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsigma_lev ! built array with sigma levels
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: de3 ! vertical metrics
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: ddepu ! depth of vel points
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dsig ! density
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dhiso ! depth of isopycns
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwtrp, dwtrpbin ! transport arrays
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dtrpbin ! transport arrays
-
- TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvar ! structure of output
-
- CHARACTER(LEN=256) :: cf_tfil ! temperature salinity file
- CHARACTER(LEN=256) :: cf_ufil ! zonal velocity file
- CHARACTER(LEN=256) :: cf_vfil ! meridional velocity file
- CHARACTER(LEN=256) :: cf_section='dens_section.dat' ! input section file
- CHARACTER(LEN=256) :: cf_out='trpsig.txt' ! output ascii file
- CHARACTER(LEN=256) :: cf_bimg ! output bimg file (2d)
- CHARACTER(LEN=256) :: cf_nc ! output netcdf file (2d)
- CHARACTER(LEN=256) :: cf_outnc ! output netcdf file (1d, 0d))
- CHARACTER(LEN=256) :: cv_dep ! depth variable
- CHARACTER(LEN=256) :: cldum ! dummy string
- CHARACTER(LEN=256) :: cglobal ! global attribute
- CHARACTER(LEN=80 ) :: cfmt_9000 ! format string
- CHARACTER(LEN=80 ) :: cfmt_9001 ! format string
- CHARACTER(LEN=80 ) :: cfmt_9002 ! format string
- CHARACTER(LEN=80 ) :: cfmt_9003 ! format string
- CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cv_names ! names of input variables
- CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: csection ! section name
-
- LOGICAL :: l_merid ! flag for meridional section
- LOGICAL :: ltemp =.FALSE. ! flag for extra print
- LOGICAL :: lprint =.FALSE. ! flag for extra print
- LOGICAL :: lbimg =.FALSE. ! flag for bimg output
- LOGICAL :: lncdf =.FALSE. ! flag for bimg output
- LOGICAL :: lfull =.FALSE. ! flag for bimg output
- LOGICAL :: lchk =.FALSE. ! flag for missing files
- !!----------------------------------------------------------------------
- CALL ReadCdfNames()
-
- narg= iargc()
- IF ( narg < 6 ) THEN
- PRINT *,' usage : cdfsigtrp T-file U-file V-file sigma_min sigma_max nbins ...'
- PRINT *,' ... [-print ] [-bimg ] [-full ] [ -refdep ref_depth] ...'
- PRINT *,' ... [-section file ] [-temp ]'
- PRINT *,' '
- PRINT *,' PURPOSE :'
- PRINT *,' Compute density class transports, according to the density class'
- PRINT *,' definition ( minimum, maximum and number of bins) given in arguments.'
- PRINT *,' Section position are given in ',TRIM(cf_section),', an ASCII file '
- PRINT *,' with pairs of lines giving section name and section location as'
- PRINT *,' imin imax jmin jmax. Only zonal or meridional section are allowed.'
- PRINT *,' The name of this file can be specified with the -section option, if'
- PRINT *,' it differs from the standard name.'
- PRINT *,' '
- PRINT *,' This program can also be used to compute transport by class of '
- PRINT *,' temperatures, provided the temperatures decreases monotonically '
- PRINT *,' downward. In this case, use -temp option and of course specify'
- PRINT *,' sigma_min, sigma_max as temperatures.'
- PRINT *,' '
- PRINT *,' ARGUMENTS :'
- PRINT *,' T-file : netcdf file with temperature and salinity'
- PRINT *,' U-file : netcdf file with zonal velocity component'
- PRINT *,' V-file : netcdf file with meridional velocity component'
- PRINT *,' sigma_min : minimum density for binning'
- PRINT *,' sigma_max : maximum density for binning'
- PRINT *,' nbins : number of bins. This will fix the bin ''width'' '
- PRINT *,' '
- PRINT *,' OPTIONS :'
- PRINT *,' [ -full ] : for full step configuration'
- PRINT *,' [ -bimg ] : produce extra bimg output file which shows the details'
- PRINT *,' of the sections (normal velocity, density, temperature, '
- PRINT *,' salinity, transports, isopycnal depths. (to be change to '
- PRINT *,' netcdf files for more common use.'
- PRINT *,' [ -ncdf ] : produce extra netcdf output file which shows the details'
- PRINT *,' of the sections (normal velocity, density, temperature, '
- PRINT *,' salinity, transports, isopycnal depths. '
- PRINT *,' [ -print ]: write the binned transports on standard output, for each'
- PRINT *,' sections.'
- PRINT *,' [ -refdep ref_depth ]: give a reference depths for the computation of'
- PRINT *,' potential density. Sigma_min, sigma_max must be adapted '
- PRINT *,' accordingly.'
- PRINT *,' [ -section file] : give the name of section file.'
- PRINT *,' Default is ', TRIM(cf_section)
- PRINT *,' [ -temp ] : use temperature instead of density for binning'
- PRINT *,' '
- PRINT *,' REQUIRED FILES :'
- PRINT *,' ', TRIM(cn_fhgr),', ', TRIM(cn_fzgr),' and ', TRIM(cf_section)
- PRINT *,' '
- PRINT *,' OUTPUT : '
- PRINT *,' Netcdf file : There is 1 netcdf file per section. File name is build'
- PRINT *,' from section name : Section_name_trpsig.nc'
- PRINT *,' variables : sigma_class (upper limit of the bin)'
- PRINT *,' sigtrp : transport (Sv per bin)'
- PRINT *,' '
- PRINT *,' ascii file : ', TRIM(cf_out)
- PRINT *,' '
- PRINT *,' bimg file : There are 2 bimg files whose name is build from section'
- PRINT *,' name : section_name_trpdep.bimg and section_name_trpsig.bimg.'
- PRINT *,' This file is written only if -bimg option is used.'
- PRINT *,' '
- PRINT *,' Standard output : the results are written on standard output only if '
- PRINT *,' the -print option is used.'
- PRINT *,' '
- PRINT *,' SEE ALSO :'
- PRINT *,' cdfrhoproj, cdftransport, cdfsigintegr '
- PRINT *,' '
- STOP
- ENDIF
-
- ! browse command line
- ijarg = 1 ; ireq = 0
- DO WHILE ( ijarg <= narg )
- CALL getarg(ijarg, cldum ) ; ijarg=ijarg+1
- SELECT CASE ( cldum )
- CASE ( '-full' ) ; lfull = .TRUE.
- CASE ( '-bimg' ) ; lbimg = .TRUE.
- CASE ( '-ncdf' ) ; lncdf = .TRUE.
- CASE ( '-print') ; lprint = .TRUE.
- CASE ( '-temp') ; ltemp = .TRUE.
- CASE ( '-refdep' ) ; CALL getarg(ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) refdep
- CASE ( '-section') ; CALL getarg(ijarg, cf_section ) ; ijarg=ijarg+1
- CASE DEFAULT
- ireq=ireq+1
- SELECT CASE ( ireq)
- CASE ( 1 ) ; cf_tfil = cldum
- CASE ( 2 ) ; cf_ufil = cldum
- CASE ( 3 ) ; cf_vfil = cldum
- CASE ( 4 ) ; READ(cldum,*) dsigma_min
- CASE ( 5 ) ; READ(cldum,*) dsigma_max
- CASE ( 6 ) ; READ(cldum,*) nbins
- CASE DEFAULT
- PRINT *,' Too many arguments ' ; STOP
- END SELECT
- END SELECT
- END DO
-
- ! check for file existence
- lchk = lchk .OR. chkfile( cn_fzgr )
- lchk = lchk .OR. chkfile( cn_fhgr )
- lchk = lchk .OR. chkfile( cf_section )
- lchk = lchk .OR. chkfile( cf_tfil )
- lchk = lchk .OR. chkfile( cf_ufil )
- lchk = lchk .OR. chkfile( cf_vfil )
- IF ( lchk ) STOP ! missing file
- IF ( ltemp) THEN ! temperature decrease downward. Change sign and swap min/max
- refdep = -10. ! flag value
- dltsig = dsigma_max ! use dltsig as dummy variable for swapping
- dsigma_max = -dsigma_min
- dsigma_min = -dltsig
- ENDIF
-
- ! define global attribute with command line
- CALL SetGlobalAtt( cglobal)
-
- ! get the attribute iweight from vozocrtx
- iweight = getatt(cf_ufil, cn_vozocrtx, 'iweight')
- IF ( iweight == 0 ) iweight = 1 ! if 0 means that it is not defined.
-
- ALLOCATE ( stypvar(nboutput), ipk(nboutput), id_varout(nboutput) )
- ALLOCATE ( rdumlon(ikx,iky), rdumlat(ikx,iky) )
-
- rdumlon(:,:)=0.
- rdumlat(:,:)=0.
-
- ipk(1)=nbins ! sigma for each level
- ipk(2)=nbins ! transport for each level
-
- ! define new variables for output
- stypvar%rmissing_value = 99999.
- stypvar%scale_factor = 1.
- stypvar%add_offset = 0.
- stypvar%savelog10 = 0.
- stypvar%iwght = iweight
- stypvar%conline_operation = 'N/A'
- stypvar%caxis = 'ZT'
-
- IF ( ltemp ) THEN
- stypvar(1)%cname = 'temp_class'
- stypvar(1)%cunits = '[]'
- stypvar(1)%valid_min = 0.
- stypvar(1)%valid_max = 100.
- stypvar(1)%clong_name = 'class of potential temperature'
- stypvar(1)%cshort_name = 'temp_class'
-
- stypvar(2)%cname = 'temptrp'
- stypvar(2)%cunits = 'Sv'
- stypvar(2)%valid_min = -1000.
- stypvar(2)%valid_max = 1000.
- stypvar(2)%clong_name = 'transport in temperature class'
- stypvar(2)%cshort_name = 'temptrp'
- ELSE
- stypvar(1)%cname = 'sigma_class'
- stypvar(1)%cunits = '[]'
- stypvar(1)%valid_min = 0.
- stypvar(1)%valid_max = 100.
- stypvar(1)%clong_name = 'class of potential density'
- stypvar(1)%cshort_name = 'sigma_class'
-
- stypvar(2)%cname = 'sigtrp'
- stypvar(2)%cunits = 'Sv'
- stypvar(2)%valid_min = -1000.
- stypvar(2)%valid_max = 1000.
- stypvar(2)%clong_name = 'transport in sigma class'
- stypvar(2)%cshort_name = 'sigtrp'
- ENDIF
-
- ! Initialise sections from file
- ! first call to get nsection and allocate arrays
- nsection = 0 ; CALL section_init(cf_section, csection, iimina, iimaxa, ijmina, ijmaxa, nsection)
- ALLOCATE ( csection(nsection), iimina(nsection), iimaxa(nsection), ijmina(nsection),ijmaxa(nsection) )
- CALL section_init(cf_section, csection,iimina,iimaxa,ijmina,ijmaxa, nsection)
-
- ! Allocate and build sigma levels and section array
- ALLOCATE ( dsigma_lev (nbins+1) , dtrpbin(nsection,nbins) )
-
- dsigma_lev(1)=dsigma_min
- dltsig=( dsigma_max - dsigma_min) / nbins
- DO jclass =2, nbins+1
- dsigma_lev(jclass)= dsigma_lev(1) + (jclass-1) * dltsig
- END DO
-
- ! Look for vertical size of the domain
- npk = getdim (cf_tfil,cn_z)
- ALLOCATE ( gdept(npk), gdepw(npk) )
- IF ( lfull ) ALLOCATE ( e3t1d(npk), e3w1d(npk))
-
- ! read gdept, gdepw : it is OK even in partial cells, as we never use the bottom gdep
- gdept(:) = getvare3(cn_fzgr, cn_gdept, npk)
- gdepw(:) = getvare3(cn_fzgr, cn_gdepw, npk)
-
- IF ( lfull ) THEN
- e3t1d(:) = getvare3(cn_fzgr, cn_ve3t, npk)
- e3w1d(:) = getvare3(cn_fzgr, cn_ve3w, npk)
- ENDIF
-
- !! * Main loop on sections
- DO jsec=1,nsection
- iimin=iimina(jsec) ; iimax=iimaxa(jsec)
- ijmin=ijmina(jsec) ; ijmax=ijmaxa(jsec)
-
- IF (iimin == iimax ) THEN ! meridional section
- npts = ijmax - ijmin ! number of segments
- l_merid = .TRUE.
-
- ELSE IF ( ijmin == ijmax ) THEN ! zonal section
- npts = iimax - iimin ! number of segments
- l_merid = .FALSE.
-
- ELSE
- PRINT *,' Section ',TRIM(csection(jsec)),' is neither zonal nor meridional :('
- PRINT *,' We skip this section .'
- CYCLE
- ENDIF
-
- ALLOCATE ( zu(npts,npk), zt(npts,npk), zs(npts,npk), dsig(npts,0:npk) )
- ALLOCATE ( eu(npts), de3(npts,npk), ddepu(npts, 0:npk), zmask(npts,npk) )
- ALLOCATE ( tmpm(1,npts,2), tmpz(npts,1,2) )
- ALLOCATE ( dwtrp(npts, nbins+1), dhiso(npts,nbins+1), dwtrpbin(npts,nbins) )
- ALLOCATE ( rlonlat(npts,1) )
-
- zt = 0. ; zs = 0. ; zu = 0. ; ddepu= 0. ; zmask = 0. ; dsig=0.d0
-
- IF (l_merid ) THEN ! meridional section at i=iimin=iimax
- tmpm(:,:,1) = getvar(cn_fhgr, cn_ve2u, 1, 1, npts, kimin=iimin, kjmin=ijmin+1)
- eu(:) = tmpm(1,:,1) ! metrics varies only horizontally
- tmpm(:,:,1) = getvar(cn_fhgr, cn_vlat2d, 1, 1, npts, kimin=iimin, kjmin=ijmin+1)
- rlonlat(:,1) = tmpm(1,:,1) ! latitude in this case
- DO jk = 1,npk
- ! initiliaze ddepu to gdept()
- ddepu(:,jk) = gdept(jk)
-
- IF ( lfull ) THEN
- de3(:,jk) = e3t1d(jk)
- tmpm(1,:,1) = e3w1d(jk)
- tmpm(1,:,2) = e3w1d(jk)
- ELSE
- ! vertical metrics (PS case)
- tmpm(:,:,1) = getvar(cn_fzgr, 'e3u_ps', jk, 1, npts, kimin=iimin, kjmin=ijmin+1, ldiom=.TRUE.)
- de3(:,jk) = tmpm(1,:,1)
- tmpm(:,:,1) = getvar(cn_fzgr, 'e3w_ps', jk, 1, npts, kimin=iimin, kjmin=ijmin+1, ldiom=.TRUE.)
- tmpm(:,:,2) = getvar(cn_fzgr, 'e3w_ps', jk, 1, npts, kimin=iimin+1, kjmin=ijmin+1, ldiom=.TRUE.)
- ENDIF
-
- IF (jk >= 2 ) THEN
- DO ji=1,npts
- ddepu(ji,jk)= ddepu(ji,jk-1) + MIN(tmpm(1,ji,1), tmpm(1,ji,2))
- END DO
- ENDIF
-
- ! Normal velocity
- tmpm(:,:,1) = getvar(cf_ufil,cn_vozocrtx,jk,1,npts, kimin=iimin, kjmin=ijmin+1)
- zu(:,jk) = tmpm(1,:,1)
-
- ! salinity and deduce umask for the section
- tmpm(:,:,1) = getvar(cf_tfil,cn_vosaline,jk,1,npts, kimin=iimin , kjmin=ijmin+1)
- tmpm(:,:,2) = getvar(cf_tfil,cn_vosaline,jk,1,npts, kimin=iimin+1, kjmin=ijmin+1)
- zmask(:,jk) = tmpm(1,:,1)*tmpm(1,:,2)
- WHERE ( zmask(:,jk) /= 0 ) zmask(:,jk)=1
- ! do not take special care for land value, as the corresponding velocity point is masked
- zs(:,jk) = 0.5 * ( tmpm(1,:,1) + tmpm(1,:,2) )
-
- ! limitation to 'wet' points
- IF ( SUM(zs(:,jk)) == 0 ) THEN
- nk=jk ! first vertical point of the section full on land
- EXIT ! as soon as all the points are on land
- ENDIF
-
- ! temperature
- tmpm(:,:,1) = getvar(cf_tfil, cn_votemper, jk, 1, npts, kimin=iimin, kjmin=ijmin+1)
- tmpm(:,:,2) = getvar(cf_tfil, cn_votemper, jk, 1, npts, kimin=iimin+1, kjmin=ijmin+1)
- zt(:,jk) = 0.5 * ( tmpm(1,:,1) + tmpm(1,:,2) )
- END DO
-
- ELSE ! zonal section at j=ijmin=ijmax
- tmpz(:,:,1) = getvar(cn_fhgr, cn_ve1v, 1, npts, 1, kimin=iimin, kjmin=ijmin)
- eu(:) = tmpz(:,1,1)
- tmpz(:,:,1) = getvar(cn_fhgr, cn_vlon2d, 1, npts, 1, kimin=iimin, kjmin=ijmin)
- rlonlat(:,1) = tmpz(:,1,1) ! longitude in this case
- DO jk=1,npk
- ! initiliaze ddepu to gdept()
- ddepu(:,jk) = gdept(jk)
-
- IF ( lfull ) THEN
- de3(:,jk) = e3t1d(jk)
- tmpm(:,1,1) = e3w1d(jk)
- tmpm(:,1,2) = e3w1d(jk)
- ELSE
- ! vertical metrics (PS case)
- tmpz(:,:,1)=getvar(cn_fzgr,'e3v_ps',jk, npts, 1, kimin=iimin+1, kjmin=ijmin, ldiom=.TRUE.)
- de3(:,jk) = tmpz(:,1,1)
- tmpz(:,:,1)=getvar(cn_fzgr,'e3w_ps',jk,npts,1, kimin=iimin+1, kjmin=ijmin, ldiom=.TRUE.)
- tmpz(:,:,2)=getvar(cn_fzgr,'e3w_ps',jk,npts,1, kimin=iimin+1, kjmin=ijmin+1, ldiom=.TRUE.)
- ENDIF
-
- IF (jk >= 2 ) THEN
- DO ji=1,npts
- ddepu(ji,jk)= ddepu(ji,jk-1) + MIN(tmpz(ji,1,1), tmpz(ji,1,2))
- END DO
- ENDIF
-
- ! Normal velocity
- tmpz(:,:,1)=getvar(cf_vfil,cn_vomecrty,jk,npts,1, kimin=iimin+1, kjmin=ijmin)
- zu(:,jk)=tmpz(:,1,1)
-
- ! salinity and mask
- tmpz(:,:,1)=getvar(cf_tfil,cn_vosaline,jk, npts, 1, kimin=iimin+1, kjmin=ijmin)
- tmpz(:,:,2)=getvar(cf_tfil,cn_vosaline,jk, npts, 1, kimin=iimin+1, kjmin=ijmin+1)
- zmask(:,jk)=tmpz(:,1,1)*tmpz(:,1,2)
- WHERE ( zmask(:,jk) /= 0 ) zmask(:,jk)=1
- ! do not take special care for land value, as the corresponding velocity point is masked
- zs(:,jk) = 0.5 * ( tmpz(:,1,1) + tmpz(:,1,2) )
-
- ! limitation to 'wet' points
- IF ( SUM(zs(:,jk)) == 0 ) THEN
- nk=jk ! first vertical point of the section full on land
- EXIT ! as soon as all the points are on land
- ENDIF
-
- ! temperature
- tmpz(:,:,1)=getvar(cf_tfil,cn_votemper,jk, npts, 1, kimin=iimin+1, kjmin=ijmin)
- tmpz(:,:,2)=getvar(cf_tfil,cn_votemper,jk, npts, 1, kimin=iimin+1, kjmin=ijmin+1)
- zt(:,jk) = 0.5 * ( tmpz(:,1,1) + tmpz(:,1,2) )
- END DO
-
- ENDIF
-
- ! compute density only for wet points
- IF ( refdep == -10. ) THEN
- dsig(:,1:nk)= -zt(:,:) ! change sign
- ELSEIF ( refdep == 0. ) THEN
- dsig(:,1:nk)=sigma0( zt, zs, npts, nk)*zmask(:,:)
- ELSE
- dsig(:,1:nk)=sigmai( zt, zs, refdep, npts, nk)*zmask(:,:)
- ENDIF
-
- dsig(:,0)=dsig(:,1)-1.e-4 ! dummy layer for easy interpolation
-
- ! compute depth of isopynals (nbins+1 )
- DO jiso =1, nbins+1
- dsigma=dsigma_lev(jiso)
+ !!======================================================================
+ !! *** PROGRAM cdfsigtrp ***
+ !!======================================================================
+ !! ** Purpose : Compute density class Mass transport across a section.
+ !!
+ !! ** Method :- The begining and end point of the section are given in
+ !! term of f-points index.
+ !! - The program works for zonal or meridional sections.
+ !! - The section definitions are given in an ASCII FILE
+ !! dens_section.dat:
+ !! foreach sections, 2 lines :
+ !! (i) : section name (String, no blank)
+ !! (ii) : imin imax jmin jmax for the section
+ !! - Only vertical slices corrsponding to the sections are
+ !! read in the files.
+ !! - read metrics, depth, etc
+ !! - read normal velocity (either vozocrtx oy vomecrty )
+ !! - read 2 rows of T and S ( i i+1 or j j+1 )
+ !! - compute the mean value at velocity point
+ !! - compute sigma0 (can be easily modified for sigmai )
+ !! - compute the depths of isopyncal surfaces
+ !! - compute the transport from surface to the isopycn
+ !! - compute the transport in each class of density
+ !! - compute the total transport (for information)
+ !!
+ !! History : 2.1 : 03/2006 : J.M. Molines : Original code
+ !! : 07/2009 : R. Dussin : add cdf output
+ !! 3.0 : 06/2011 : J.M. Molines : Doctor norm + Lic.
+ !!----------------------------------------------------------------------
+ !!----------------------------------------------------------------------
+ !! routines : description
+ !! section_init : initialize section names and positions
+ !! print_out : routine which performs standard output if required
+ !! bimg_writ : routine which performs bimg output if required
+ !!----------------------------------------------------------------------
+ USE cdfio
+ USE eos ! for sigma0, sigmai
+ USE modcdfnames ! for ReadCdfNames
+ USE modutils ! for SetGlobalAtt
+ !!----------------------------------------------------------------------
+ !! CDFTOOLS_3.0 , MEOM 2011
+ !! $Id$
+ !! Copyright (c) 2011, J.-M. Molines
+ !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
+ !!----------------------------------------------------------------------
+ IMPLICIT NONE
+
+ INTEGER(KIND=4) :: ji, jk, jclass, jsec ! dummy loop index
+ INTEGER(KIND=4) :: jiso, jbin, jarg ! dummy loop index
+ INTEGER(KIND=4) :: nbins ! number of density classes
+ INTEGER(KIND=4) :: ipos ! working variable
+ INTEGER(KIND=4) :: narg, iargc ! command line
+ INTEGER(KIND=4) :: ijarg, ireq ! command line
+ INTEGER(KIND=4) :: npk, nk ! vertical size, number of wet layers
+ INTEGER(KIND=4) :: numbimg=10 ! optional bimg logical unit
+ INTEGER(KIND=4) :: numout=11 ! ascii output
+ INTEGER(KIND=4) :: nsection ! number of sections (overall)
+ INTEGER(KIND=4) :: iimin, iimax ! working section limits
+ INTEGER(KIND=4) :: ijmin, ijmax ! working section limits
+ INTEGER(KIND=4) :: npts ! number of points in section
+ INTEGER(KIND=4) :: ikx=1, iky=1 ! dims of netcdf output file
+ INTEGER(KIND=4) :: nboutput=2 ! number of values to write in cdf output
+ INTEGER(KIND=4) :: ncout, ierr ! for netcdf output
+ INTEGER(KIND=4) :: iweight ! weight of input file for further averaging
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: iimina, iimaxa ! sections limits
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ijmina, ijmaxa ! sections limits
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_varout ! variable levels and id
+
+ REAL(KIND=4) :: refdep =0.e0 ! reference depth (m)
+ REAL(KIND=4), DIMENSION(1) :: tim ! time counter
+ REAL(KIND=4), DIMENSION(1) :: rdummy1, rdummy2 ! working variable
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdept, gdepw ! depth of T and W points
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: eu ! either e1v or e2u
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: e3t1d, e3w1d ! vertical metrics in case of full step
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rlonlat ! longitudes/latitudes if the section
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zs, zt ! salinity and temperature from file
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rdumlon, rdumlat ! dummy longitude and latitude for output
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zu ! velocity
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zmask ! mask
+ REAL(KIND=4), DIMENSION(:,:,:), ALLOCATABLE :: tmpm, tmpz ! temporary arrays
+
+ ! double precision for cumulative variables and densities
+ REAL(KIND=8) :: dsigma_min ! minimum density for bining
+ REAL(KIND=8) :: dsigma_max, dltsig ! maximum density for bining, step
+ REAL(KIND=8) :: dsigma, dalfa ! working sigma, interpolation coeff.
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsigma_lev ! built array with sigma levels
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: de3 ! vertical metrics
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: ddepu ! depth of vel points
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dsig ! density
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dhiso ! depth of isopycns
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwtrp, dwtrpbin ! transport arrays
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dtrpbin ! transport arrays
+
+ TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvar ! structure of output
+
+ CHARACTER(LEN=256) :: cf_tfil ! temperature salinity file
+ CHARACTER(LEN=256) :: cf_ufil ! zonal velocity file
+ CHARACTER(LEN=256) :: cf_vfil ! meridional velocity file
+ CHARACTER(LEN=256) :: cf_section='dens_section.dat' ! input section file
+ CHARACTER(LEN=256) :: cf_out='trpsig.txt' ! output ascii file
+ CHARACTER(LEN=256) :: cf_bimg ! output bimg file (2d)
+ CHARACTER(LEN=256) :: cf_nc ! output netcdf file (2d)
+ CHARACTER(LEN=256) :: cf_outnc ! output netcdf file (1d, 0d))
+ CHARACTER(LEN=256) :: cv_dep ! depth variable
+ CHARACTER(LEN=256) :: cldum ! dummy string
+ CHARACTER(LEN=256) :: cglobal ! global attribute
+ CHARACTER(LEN=80 ) :: cfmt_9000 ! format string
+ CHARACTER(LEN=80 ) :: cfmt_9001 ! format string
+ CHARACTER(LEN=80 ) :: cfmt_9002 ! format string
+ CHARACTER(LEN=80 ) :: cfmt_9003 ! format string
+ CHARACTER(LEN=256) :: cl_vnam, cl_lname ! working variables
+ CHARACTER(LEN=256) :: csuffixvarname !
+ CHARACTER(LEN=256) :: cprefixlongname !
+ CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cv_names ! names of input variables
+ CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: csection ! section name
+ CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cvarname ! output variable name (root)
+ CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: clongname ! output long name (root)
+
+ LOGICAL :: l_merid ! flag for meridional section
+ LOGICAL :: ltemp =.FALSE. ! flag for extra print
+ LOGICAL :: lprint =.FALSE. ! flag for extra print
+ LOGICAL :: lbimg =.FALSE. ! flag for bimg output
+ LOGICAL :: lncdf =.FALSE. ! flag for bimg output
+ LOGICAL :: lfull =.FALSE. ! flag for bimg output
+ LOGICAL :: lchk =.FALSE. ! flag for missing files
+ !!----------------------------------------------------------------------
+ CALL ReadCdfNames()
+
+ narg= iargc()
+ IF ( narg < 6 ) THEN
+ PRINT *,' usage : cdfsigtrp T-file U-file V-file sigma_min sigma_max nbins ...'
+ PRINT *,' ... [-print ] [-bimg ] [-full ] [ -refdep ref_depth] ...'
+ PRINT *,' ... [-section file ] [-temp ]'
+ PRINT *,' '
+ PRINT *,' PURPOSE :'
+ PRINT *,' Compute density class transports, according to the density class'
+ PRINT *,' definition ( minimum, maximum and number of bins) given in arguments.'
+ PRINT *,' Section position are given in ',TRIM(cf_section),', an ASCII file '
+ PRINT *,' with pairs of lines giving section name and section location as'
+ PRINT *,' imin imax jmin jmax. Only zonal or meridional section are allowed.'
+ PRINT *,' The name of this file can be specified with the -section option, if'
+ PRINT *,' it differs from the standard name. Optionaly, a netcdf root variable '
+ PRINT *,' name and a netcdf root long-name can be provided on the line giving '
+ PRINT *,' the section name.'
+ PRINT *,' '
+ PRINT *,' This program can also be used to compute transport by class of '
+ PRINT *,' temperatures, provided the temperatures decreases monotonically '
+ PRINT *,' downward. In this case, use -temp option and of course specify'
+ PRINT *,' sigma_min, sigma_max as temperatures.'
+ PRINT *,' '
+ PRINT *,' ARGUMENTS :'
+ PRINT *,' T-file : netcdf file with temperature and salinity'
+ PRINT *,' U-file : netcdf file with zonal velocity component'
+ PRINT *,' V-file : netcdf file with meridional velocity component'
+ PRINT *,' sigma_min : minimum density for binning'
+ PRINT *,' sigma_max : maximum density for binning'
+ PRINT *,' nbins : number of bins. This will fix the bin ''width'' '
+ PRINT *,' '
+ PRINT *,' OPTIONS :'
+ PRINT *,' [ -full ] : for full step configuration'
+ PRINT *,' [ -bimg ] : produce extra bimg output file which shows the details'
+ PRINT *,' of the sections (normal velocity, density, temperature, '
+ PRINT *,' salinity, transports, isopycnal depths. (to be change to '
+ PRINT *,' netcdf files for more common use.'
+ PRINT *,' [ -ncdf ] : produce extra netcdf output file which shows the details'
+ PRINT *,' of the sections (normal velocity, density, temperature, '
+ PRINT *,' salinity, transports, isopycnal depths. '
+ PRINT *,' [ -print ]: write the binned transports on standard output, for each'
+ PRINT *,' sections.'
+ PRINT *,' [ -refdep ref_depth ]: give a reference depths for the computation of'
+ PRINT *,' potential density. Sigma_min, sigma_max must be adapted '
+ PRINT *,' accordingly.'
+ PRINT *,' [ -section file] : give the name of section file.'
+ PRINT *,' Default is ', TRIM(cf_section)
+ PRINT *,' [ -temp ] : use temperature instead of density for binning'
+ PRINT *,' '
+ PRINT *,' REQUIRED FILES :'
+ PRINT *,' ', TRIM(cn_fhgr),', ', TRIM(cn_fzgr),' and ', TRIM(cf_section)
+ PRINT *,' '
+ PRINT *,' OUTPUT : '
+ PRINT *,' Netcdf file : There is 1 netcdf file per section. File name is build'
+ PRINT *,' from section name : Section_name_trpsig.nc'
+ PRINT *,' variables : sigma_class (upper limit of the bin)'
+ PRINT *,' sigtrp : transport (Sv per bin)'
+ PRINT *,' '
+ PRINT *,' ascii file : ', TRIM(cf_out)
+ PRINT *,' '
+ PRINT *,' bimg file : There are 2 bimg files whose name is build from section'
+ PRINT *,' name : section_name_trpdep.bimg and section_name_trpsig.bimg.'
+ PRINT *,' This file is written only if -bimg option is used.'
+ PRINT *,' '
+ PRINT *,' Standard output : the results are written on standard output only if '
+ PRINT *,' the -print option is used.'
+ PRINT *,' '
+ PRINT *,' SEE ALSO :'
+ PRINT *,' cdfrhoproj, cdftransport, cdfsigintegr '
+ PRINT *,' '
+ STOP
+ ENDIF
+
+ ! browse command line
+ ijarg = 1 ; ireq = 0
+ DO WHILE ( ijarg <= narg )
+ CALL getarg(ijarg, cldum ) ; ijarg=ijarg+1
+ SELECT CASE ( cldum )
+ CASE ( '-full' ) ; lfull = .TRUE.
+ CASE ( '-bimg' ) ; lbimg = .TRUE.
+ CASE ( '-ncdf' ) ; lncdf = .TRUE.
+ CASE ( '-print') ; lprint = .TRUE.
+ CASE ( '-temp') ; ltemp = .TRUE.
+ CASE ( '-refdep' ) ; CALL getarg(ijarg, cldum ) ; ijarg=ijarg+1 ; READ(cldum,*) refdep
+ CASE ( '-section') ; CALL getarg(ijarg, cf_section ) ; ijarg=ijarg+1
+ CASE DEFAULT
+ ireq=ireq+1
+ SELECT CASE ( ireq)
+ CASE ( 1 ) ; cf_tfil = cldum
+ CASE ( 2 ) ; cf_ufil = cldum
+ CASE ( 3 ) ; cf_vfil = cldum
+ CASE ( 4 ) ; READ(cldum,*) dsigma_min
+ CASE ( 5 ) ; READ(cldum,*) dsigma_max
+ CASE ( 6 ) ; READ(cldum,*) nbins
+ CASE DEFAULT
+ PRINT *,' Too many arguments ' ; STOP
+ END SELECT
+ END SELECT
+ END DO
+
+ ! check for file existence
+ lchk = lchk .OR. chkfile( cn_fzgr )
+ lchk = lchk .OR. chkfile( cn_fhgr )
+ lchk = lchk .OR. chkfile( cf_section )
+ lchk = lchk .OR. chkfile( cf_tfil )
+ lchk = lchk .OR. chkfile( cf_ufil )
+ lchk = lchk .OR. chkfile( cf_vfil )
+ IF ( lchk ) STOP ! missing file
+ IF ( ltemp) THEN ! temperature decrease downward. Change sign and swap min/max
+ refdep = -10. ! flag value
+ dltsig = dsigma_max ! use dltsig as dummy variable for swapping
+ dsigma_max = -dsigma_min
+ dsigma_min = -dltsig
+ ENDIF
+
+ ! define global attribute with command line
+ CALL SetGlobalAtt( cglobal)
+
+ ! get the attribute iweight from vozocrtx
+ iweight = getatt(cf_ufil, cn_vozocrtx, 'iweight')
+ IF ( iweight == 0 ) iweight = 1 ! if 0 means that it is not defined.
+
+ ALLOCATE ( stypvar(nboutput), ipk(nboutput), id_varout(nboutput) )
+ ALLOCATE ( rdumlon(ikx,iky), rdumlat(ikx,iky) )
+
+ rdumlon(:,:)=0.
+ rdumlat(:,:)=0.
+
+ ipk(1)=nbins ! sigma for each level
+ ipk(2)=nbins ! transport for each level
+ ! initialisation of variable names etc... is done according to section name
+
+ ! Initialise sections from file
+ ! first call to get nsection and allocate arrays
+ nsection = 0 ; CALL section_init(cf_section, csection,cvarname,clongname,iimina, iimaxa, ijmina, ijmaxa, nsection)
+ ALLOCATE ( csection(nsection), cvarname(nsection), clongname(nsection) )
+ ALLOCATE ( iimina(nsection), iimaxa(nsection), ijmina(nsection),ijmaxa(nsection) )
+ CALL section_init(cf_section, csection,cvarname,clongname, iimina,iimaxa,ijmina,ijmaxa, nsection)
+
+ ! Allocate and build sigma levels and section array
+ ALLOCATE ( dsigma_lev (nbins+1) , dtrpbin(nsection,nbins) )
+
+ dsigma_lev(1)=dsigma_min
+ dltsig=( dsigma_max - dsigma_min) / nbins
+ DO jclass =2, nbins+1
+ dsigma_lev(jclass)= dsigma_lev(1) + (jclass-1) * dltsig
+ END DO
+
+ ! Look for vertical size of the domain
+ npk = getdim (cf_tfil,cn_z)
+ ALLOCATE ( gdept(npk), gdepw(npk) )
+ IF ( lfull ) ALLOCATE ( e3t1d(npk), e3w1d(npk))
+
+ ! read gdept, gdepw : it is OK even in partial cells, as we never use the bottom gdep
+ gdept(:) = getvare3(cn_fzgr, cn_gdept, npk)
+ gdepw(:) = getvare3(cn_fzgr, cn_gdepw, npk)
+
+ IF ( lfull ) THEN
+ e3t1d(:) = getvare3(cn_fzgr, cn_ve3t, npk)
+ e3w1d(:) = getvare3(cn_fzgr, cn_ve3w, npk)
+ ENDIF
+
+ !! * Main loop on sections
+ DO jsec=1,nsection
+ iimin=iimina(jsec) ; iimax=iimaxa(jsec)
+ ijmin=ijmina(jsec) ; ijmax=ijmaxa(jsec)
+
+ IF (iimin == iimax ) THEN ! meridional section
+ npts = ijmax - ijmin ! number of segments
+ l_merid = .TRUE.
+
+ ELSE IF ( ijmin == ijmax ) THEN ! zonal section
+ npts = iimax - iimin ! number of segments
+ l_merid = .FALSE.
+
+ ELSE
+ PRINT *,' Section ',TRIM(csection(jsec)),' is neither zonal nor meridional :('
+ PRINT *,' We skip this section .'
+ CYCLE
+ ENDIF
+
+ ALLOCATE ( zu(npts,npk), zt(npts,npk), zs(npts,npk), dsig(npts,0:npk) )
+ ALLOCATE ( eu(npts), de3(npts,npk), ddepu(npts, 0:npk), zmask(npts,npk) )
+ ALLOCATE ( tmpm(1,npts,2), tmpz(npts,1,2) )
+ ALLOCATE ( dwtrp(npts, nbins+1), dhiso(npts,nbins+1), dwtrpbin(npts,nbins) )
+ ALLOCATE ( rlonlat(npts,1) )
+
+ zt = 0. ; zs = 0. ; zu = 0. ; ddepu= 0. ; zmask = 0. ; dsig=0.d0
+
+ IF (l_merid ) THEN ! meridional section at i=iimin=iimax
+ tmpm(:,:,1) = getvar(cn_fhgr, cn_ve2u, 1, 1, npts, kimin=iimin, kjmin=ijmin+1)
+ eu(:) = tmpm(1,:,1) ! metrics varies only horizontally
+ tmpm(:,:,1) = getvar(cn_fhgr, cn_vlat2d, 1, 1, npts, kimin=iimin, kjmin=ijmin+1)
+ rlonlat(:,1) = tmpm(1,:,1) ! latitude in this case
+ DO jk = 1,npk
+ ! initiliaze ddepu to gdept()
+ ddepu(:,jk) = gdept(jk)
+
+ IF ( lfull ) THEN
+ de3(:,jk) = e3t1d(jk)
+ tmpm(1,:,1) = e3w1d(jk)
+ tmpm(1,:,2) = e3w1d(jk)
+ ELSE
+ ! vertical metrics (PS case)
+ tmpm(:,:,1) = getvar(cn_fzgr, 'e3u_ps', jk, 1, npts, kimin=iimin, kjmin=ijmin+1, ldiom=.TRUE.)
+ de3(:,jk) = tmpm(1,:,1)
+ tmpm(:,:,1) = getvar(cn_fzgr, 'e3w_ps', jk, 1, npts, kimin=iimin, kjmin=ijmin+1, ldiom=.TRUE.)
+ tmpm(:,:,2) = getvar(cn_fzgr, 'e3w_ps', jk, 1, npts, kimin=iimin+1, kjmin=ijmin+1, ldiom=.TRUE.)
+ ENDIF
+
+ IF (jk >= 2 ) THEN
+ DO ji=1,npts
+ ddepu(ji,jk)= ddepu(ji,jk-1) + MIN(tmpm(1,ji,1), tmpm(1,ji,2))
+ END DO
+ ENDIF
+
+ ! Normal velocity
+ tmpm(:,:,1) = getvar(cf_ufil,cn_vozocrtx,jk,1,npts, kimin=iimin, kjmin=ijmin+1)
+ zu(:,jk) = tmpm(1,:,1)
+
+ ! salinity and deduce umask for the section
+ tmpm(:,:,1) = getvar(cf_tfil,cn_vosaline,jk,1,npts, kimin=iimin , kjmin=ijmin+1)
+ tmpm(:,:,2) = getvar(cf_tfil,cn_vosaline,jk,1,npts, kimin=iimin+1, kjmin=ijmin+1)
+ zmask(:,jk) = tmpm(1,:,1)*tmpm(1,:,2)
+ WHERE ( zmask(:,jk) /= 0 ) zmask(:,jk)=1
+ ! do not take special care for land value, as the corresponding velocity point is masked
+ zs(:,jk) = 0.5 * ( tmpm(1,:,1) + tmpm(1,:,2) )
+
+ ! limitation to 'wet' points
+ IF ( SUM(zs(:,jk)) == 0 ) THEN
+ nk=jk ! first vertical point of the section full on land
+ EXIT ! as soon as all the points are on land
+ ENDIF
+
+ ! temperature
+ tmpm(:,:,1) = getvar(cf_tfil, cn_votemper, jk, 1, npts, kimin=iimin, kjmin=ijmin+1)
+ tmpm(:,:,2) = getvar(cf_tfil, cn_votemper, jk, 1, npts, kimin=iimin+1, kjmin=ijmin+1)
+ zt(:,jk) = 0.5 * ( tmpm(1,:,1) + tmpm(1,:,2) )
+ END DO
+
+ ELSE ! zonal section at j=ijmin=ijmax
+ tmpz(:,:,1) = getvar(cn_fhgr, cn_ve1v, 1, npts, 1, kimin=iimin, kjmin=ijmin)
+ eu(:) = tmpz(:,1,1)
+ tmpz(:,:,1) = getvar(cn_fhgr, cn_vlon2d, 1, npts, 1, kimin=iimin, kjmin=ijmin)
+ rlonlat(:,1) = tmpz(:,1,1) ! longitude in this case
+ DO jk=1,npk
+ ! initiliaze ddepu to gdept()
+ ddepu(:,jk) = gdept(jk)
+
+ IF ( lfull ) THEN
+ de3(:,jk) = e3t1d(jk)
+ tmpm(:,1,1) = e3w1d(jk)
+ tmpm(:,1,2) = e3w1d(jk)
+ ELSE
+ ! vertical metrics (PS case)
+ tmpz(:,:,1)=getvar(cn_fzgr,'e3v_ps',jk, npts, 1, kimin=iimin+1, kjmin=ijmin, ldiom=.TRUE.)
+ de3(:,jk) = tmpz(:,1,1)
+ tmpz(:,:,1)=getvar(cn_fzgr,'e3w_ps',jk,npts,1, kimin=iimin+1, kjmin=ijmin, ldiom=.TRUE.)
+ tmpz(:,:,2)=getvar(cn_fzgr,'e3w_ps',jk,npts,1, kimin=iimin+1, kjmin=ijmin+1, ldiom=.TRUE.)
+ ENDIF
+
+ IF (jk >= 2 ) THEN
+ DO ji=1,npts
+ ddepu(ji,jk)= ddepu(ji,jk-1) + MIN(tmpz(ji,1,1), tmpz(ji,1,2))
+ END DO
+ ENDIF
+
+ ! Normal velocity
+ tmpz(:,:,1)=getvar(cf_vfil,cn_vomecrty,jk,npts,1, kimin=iimin+1, kjmin=ijmin)
+ zu(:,jk)=tmpz(:,1,1)
+
+ ! salinity and mask
+ tmpz(:,:,1)=getvar(cf_tfil,cn_vosaline,jk, npts, 1, kimin=iimin+1, kjmin=ijmin)
+ tmpz(:,:,2)=getvar(cf_tfil,cn_vosaline,jk, npts, 1, kimin=iimin+1, kjmin=ijmin+1)
+ zmask(:,jk)=tmpz(:,1,1)*tmpz(:,1,2)
+ WHERE ( zmask(:,jk) /= 0 ) zmask(:,jk)=1
+ ! do not take special care for land value, as the corresponding velocity point is masked
+ zs(:,jk) = 0.5 * ( tmpz(:,1,1) + tmpz(:,1,2) )
+
+ ! limitation to 'wet' points
+ IF ( SUM(zs(:,jk)) == 0 ) THEN
+ nk=jk ! first vertical point of the section full on land
+ EXIT ! as soon as all the points are on land
+ ENDIF
+
+ ! temperature
+ tmpz(:,:,1)=getvar(cf_tfil,cn_votemper,jk, npts, 1, kimin=iimin+1, kjmin=ijmin)
+ tmpz(:,:,2)=getvar(cf_tfil,cn_votemper,jk, npts, 1, kimin=iimin+1, kjmin=ijmin+1)
+ zt(:,jk) = 0.5 * ( tmpz(:,1,1) + tmpz(:,1,2) )
+ END DO
+
+ ENDIF
+
+ ! compute density only for wet points
+ IF ( refdep == -10. ) THEN
+ dsig(:,1:nk)= -zt(:,:) ! change sign
+ ELSEIF ( refdep == 0. ) THEN
+ dsig(:,1:nk)=sigma0( zt, zs, npts, nk)*zmask(:,:)
+ ELSE
+ dsig(:,1:nk)=sigmai( zt, zs, refdep, npts, nk)*zmask(:,:)
+ ENDIF
+
+ dsig(:,0)=dsig(:,1)-1.e-4 ! dummy layer for easy interpolation
+
+ ! compute depth of isopynals (nbins+1 )
+ DO jiso =1, nbins+1
+ dsigma=dsigma_lev(jiso)
!!! REM : I and K loop can be inverted if necessary
- DO ji=1,npts
- dhiso(ji,jiso) = gdept(npk)
- DO jk=1,nk
- IF ( dsig(ji,jk) < dsigma ) THEN
- ELSE
- ! interpolate between jk-1 and jk
- dalfa=(dsigma - dsig(ji,jk-1)) / ( dsig(ji,jk) -dsig(ji,jk-1) )
- IF (ABS(dalfa) > 1.1 .OR. dalfa < 0 ) THEN ! case dsig(0) = dsig(1)-1.e-4
- dhiso(ji,jiso)= 0.d0
- ELSE
- dhiso(ji,jiso)= ddepu(ji,jk)*dalfa + (1.d0-dalfa)* ddepu(ji,jk-1)
- ENDIF
- EXIT
- ENDIF
- END DO
- END DO
- END DO
-
- ! compute transport between surface and isopycn
- DO jiso = 1, nbins + 1
- dsigma=dsigma_lev(jiso)
- DO ji=1,npts
- dwtrp(ji,jiso) = 0.d0
- DO jk=1, nk-1
- IF ( gdepw(jk+1) < dhiso(ji,jiso) ) THEN
- dwtrp(ji,jiso)= dwtrp(ji,jiso) + eu(ji)*de3(ji,jk)*zu(ji,jk)*1.d0
- ELSE ! last box ( fraction)
- dwtrp(ji,jiso)= dwtrp(ji,jiso) + eu(ji)*(dhiso(ji,jiso)-gdepw(jk))*zu(ji,jk)*1.d0
- EXIT ! jk loop
- ENDIF
- END DO
- END DO
- END DO
-
- ! binned transport : difference between 2 isopycns
- DO jbin=1, nbins
- dsigma=dsigma_lev(jbin)
- DO ji=1, npts
- dwtrpbin(ji,jbin) = dwtrp(ji,jbin+1) - dwtrp(ji,jbin)
- END DO
- dtrpbin(jsec,jbin)=SUM(dwtrpbin(:,jbin) )
- END DO
-
- ! output of the code for 1 section
- IF (lprint) CALL print_out(jsec)
- IF (lbimg ) CALL bimg_writ(jsec)
- IF (lncdf ) CALL cdf_writ(jsec)
- PRINT *,' Total transport in all bins :',TRIM(csection(jsec)),' ',SUM(dtrpbin(jsec,:) )/1.d6
-
- ! free memory for the next section
- DEALLOCATE ( zu, zt, zs, dsig, ddepu, dhiso, dwtrp, dwtrpbin )
- DEALLOCATE ( eu, de3, tmpm, tmpz, zmask, rlonlat )
-
- END DO ! next section
-
- !! Global Output
- OPEN( numout, FILE=cf_out)
- ipos=INDEX(cf_tfil,'_gridT.nc')
- WRITE(numout,9006) TRIM(cf_tfil(1:ipos-1))
- WRITE(numout,9005) ' sigma ', (csection(jsec),jsec=1,nsection)
- DO jiso=1,nbins
- WRITE(numout,9004) dsigma_lev(jiso), (dtrpbin(jsec,jiso),jsec=1,nsection)
- ENDDO
- CLOSE(numout)
-
- cv_dep='levels'
-
- DO jsec=1,nsection
- ! create output fileset
- IF (ltemp) THEN
- cf_outnc = TRIM(csection(jsec))//'_trptemp.nc'
- ELSE
- cf_outnc = TRIM(csection(jsec))//'_trpsig.nc'
- ENDIF
-
- ncout = create (cf_outnc, 'none', ikx, iky, nbins, cdep=cv_dep )
- ierr = createvar (ncout, stypvar, nboutput, ipk, id_varout, cdglobal=TRIM(cglobal))
- ierr = putheadervar(ncout, cf_tfil, ikx, iky, nbins, &
- & pnavlon=rdumlon, pnavlat=rdumlat, pdep=REAL(dsigma_lev), cdep=cv_dep )
-
- tim = getvar1d(cf_tfil, cn_vtimec, 1 )
- ierr = putvar1d(ncout, tim, 1, 'T')
-
- DO jiso=1,nbins
- rdummy1 = dsigma_lev(jiso)
- rdummy2 = dtrpbin(jsec,jiso)/1.d6 ! Sv
- ierr = putvar(ncout, id_varout(1), rdummy1, jiso, ikx, iky )
- ierr = putvar(ncout, id_varout(2), rdummy2, jiso, ikx, iky )
- END DO
-
- ierr = closeout(ncout)
-
- END DO
+ DO ji=1,npts
+ dhiso(ji,jiso) = gdept(npk)
+ DO jk=1,nk
+ IF ( dsig(ji,jk) < dsigma ) THEN
+ ELSE
+ ! interpolate between jk-1 and jk
+ dalfa=(dsigma - dsig(ji,jk-1)) / ( dsig(ji,jk) -dsig(ji,jk-1) )
+ IF (ABS(dalfa) > 1.1 .OR. dalfa < 0 ) THEN ! case dsig(0) = dsig(1)-1.e-4
+ dhiso(ji,jiso)= 0.d0
+ ELSE
+ dhiso(ji,jiso)= ddepu(ji,jk)*dalfa + (1.d0-dalfa)* ddepu(ji,jk-1)
+ ENDIF
+ EXIT
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ ! compute transport between surface and isopycn
+ DO jiso = 1, nbins + 1
+ dsigma=dsigma_lev(jiso)
+ DO ji=1,npts
+ dwtrp(ji,jiso) = 0.d0
+ DO jk=1, nk-1
+ IF ( gdepw(jk+1) < dhiso(ji,jiso) ) THEN
+ dwtrp(ji,jiso)= dwtrp(ji,jiso) + eu(ji)*de3(ji,jk)*zu(ji,jk)*1.d0
+ ELSE ! last box ( fraction)
+ dwtrp(ji,jiso)= dwtrp(ji,jiso) + eu(ji)*(dhiso(ji,jiso)-gdepw(jk))*zu(ji,jk)*1.d0
+ EXIT ! jk loop
+ ENDIF
+ END DO
+ END DO
+ END DO
+
+ ! binned transport : difference between 2 isopycns
+ DO jbin=1, nbins
+ dsigma=dsigma_lev(jbin)
+ DO ji=1, npts
+ dwtrpbin(ji,jbin) = dwtrp(ji,jbin+1) - dwtrp(ji,jbin)
+ END DO
+ dtrpbin(jsec,jbin)=SUM(dwtrpbin(:,jbin) )
+ END DO
+
+ ! output of the code for 1 section
+ IF (lprint) CALL print_out(jsec)
+ IF (lbimg ) CALL bimg_writ(jsec)
+ IF (lncdf ) CALL cdf_writ(jsec)
+ PRINT *,' Total transport in all bins :',TRIM(csection(jsec)),' ',SUM(dtrpbin(jsec,:) )/1.d6
+
+ ! free memory for the next section
+ DEALLOCATE ( zu, zt, zs, dsig, ddepu, dhiso, dwtrp, dwtrpbin )
+ DEALLOCATE ( eu, de3, tmpm, tmpz, zmask, rlonlat )
+
+ END DO ! next section
+
+ !! Global Output
+ OPEN( numout, FILE=cf_out)
+ ipos=INDEX(cf_tfil,'_gridT.nc')
+ WRITE(numout,9006) TRIM(cf_tfil(1:ipos-1))
+ WRITE(numout,9005) ' sigma ', (csection(jsec),jsec=1,nsection)
+ DO jiso=1,nbins
+ WRITE(numout,9004) dsigma_lev(jiso), (dtrpbin(jsec,jiso),jsec=1,nsection)
+ ENDDO
+ CLOSE(numout)
+
+ cv_dep='levels'
+ ! need to call section_init again in order to reset cvarname, clongname if they where modified
+ ! previously in cdf_writ( case lncdf=true )
+ IF (lncdf) THEN
+ CALL section_init(cf_section, csection,cvarname,clongname, iimina,iimaxa,ijmina,ijmaxa, nsection)
+ ENDIF
+
+ DO jsec=1,nsection
+ ! setup output variables (section dependant for adaptative variable name (if possible)
+ ! define new variables for output
+ IF ( cvarname(jsec) /= 'none' ) THEN
+ csuffixvarname='_'//TRIM(cvarname(jsec))
+ ELSE
+ csuffixvarname=''
+ ENDIF
+ IF ( clongname(jsec) /= 'none' ) THEN
+ cprefixlongname=TRIM(clongname(jsec))//'_'
+ ELSE
+ cprefixlongname=''
+ ENDIF
+
+ stypvar%rmissing_value = 99999.
+ stypvar%scale_factor = 1.
+ stypvar%add_offset = 0.
+ stypvar%savelog10 = 0.
+ stypvar%iwght = iweight
+ stypvar%conline_operation = 'N/A'
+ stypvar%caxis = 'ZT'
+
+ IF ( ltemp ) THEN
+ stypvar(1)%cname = 'temp_class'
+ stypvar(1)%cunits = '[]'
+ stypvar(1)%valid_min = 0.
+ stypvar(1)%valid_max = 100.
+ stypvar(1)%clong_name = 'class of potential temperature'
+ stypvar(1)%cshort_name = 'temp_class'
+
+ stypvar(2)%cname = 'temptrp'//TRIM(csuffixvarname)
+ stypvar(2)%cunits = 'Sv'
+ stypvar(2)%valid_min = -1000.
+ stypvar(2)%valid_max = 1000.
+ stypvar(2)%clong_name = TRIM(cprefixlongname)//'transport in temperature class'
+ stypvar(2)%cshort_name = 'temptrp'
+ ELSE
+ stypvar(1)%cname = 'sigma_class'
+ stypvar(1)%cunits = '[]'
+ stypvar(1)%valid_min = 0.
+ stypvar(1)%valid_max = 100.
+ stypvar(1)%clong_name = 'class of potential density'
+ stypvar(1)%cshort_name = 'sigma_class'
+
+ stypvar(2)%cname = 'sigtrp'//TRIM(csuffixvarname)
+ stypvar(2)%cunits = 'Sv'
+ stypvar(2)%valid_min = -1000.
+ stypvar(2)%valid_max = 1000.
+ stypvar(2)%clong_name = TRIM(cprefixlongname)//'transport in sigma class'
+ stypvar(2)%cshort_name = 'sigtrp'
+ ENDIF
+
+
+ ! create output fileset
+ IF (ltemp) THEN
+ cf_outnc = TRIM(csection(jsec))//'_trptemp.nc'
+ ELSE
+ cf_outnc = TRIM(csection(jsec))//'_trpsig.nc'
+ ENDIF
+
+ ncout = create (cf_outnc, 'none', ikx, iky, nbins, cdep=cv_dep )
+ ierr = createvar (ncout, stypvar, nboutput, ipk, id_varout, cdglobal=TRIM(cglobal))
+ ierr = putheadervar(ncout, cf_tfil, ikx, iky, nbins, &
+ & pnavlon=rdumlon, pnavlat=rdumlat, pdep=REAL(dsigma_lev), cdep=cv_dep )
+
+ tim = getvar1d(cf_tfil, cn_vtimec, 1 )
+ ierr = putvar1d(ncout, tim, 1, 'T')
+
+ DO jiso=1,nbins
+ rdummy1 = dsigma_lev(jiso)
+ rdummy2 = dtrpbin(jsec,jiso)/1.d6 ! Sv
+ ierr = putvar(ncout, id_varout(1), rdummy1, jiso, ikx, iky )
+ ierr = putvar(ncout, id_varout(2), rdummy2, jiso, ikx, iky )
+ END DO
+
+ ierr = closeout(ncout)
+
+ END DO
9004 FORMAT(f9.4, 20e16.7)
9005 FORMAT('#',a9, 20(2x,a12,2x) )
9006 FORMAT('# ',a)
CONTAINS
- SUBROUTINE section_init(cdfile, cdsection, kimin, kimax, kjmin, kjmax, knumber)
- !!---------------------------------------------------------------------
- !! *** ROUTINE section_init ***
- !!
- !! ** Purpose : Read input ASCII file that defines section names and limit of
- !! sections.
- !!
- !! ** Method : At fisrt call only return the number of sections for further
- !! allocation.
- !!
- !!----------------------------------------------------------------------
- CHARACTER(LEN=*), INTENT(in ) :: cdfile
- CHARACTER(LEN=256), DIMENSION(knumber), INTENT(out ) :: cdsection
- INTEGER(KIND=4), INTENT(inout) :: knumber
- INTEGER(KIND=4), DIMENSION(knumber), INTENT(out ) :: kimin, kimax, kjmin, kjmax
-
- ! Local variables
- INTEGER(KIND=4) :: jsec
- INTEGER(KIND=4) :: ii, inum=10
- CHARACTER(LEN=256) :: cline
- LOGICAL :: llfirst
- !!----------------------------------------------------------------------
- llfirst=.FALSE.
- IF ( knumber == 0 ) llfirst=.TRUE.
-
- OPEN(inum, FILE=cdfile)
- REWIND(inum)
- ii = 0
-
- ! read the file just to count the number of sections
- DO
- READ(inum,'(a)') cline
- IF (INDEX(cline,'EOF') == 0 ) THEN
- READ(inum,*) ! skip one line
- ii = ii + 1
- ELSE
- EXIT
- ENDIF
- END DO
-
- knumber=ii
- IF ( llfirst ) RETURN
-
- REWIND(inum)
- DO jsec=1,knumber
- READ(inum,'(a)') cdsection(jsec)
- READ(inum,* ) kimin(jsec), kimax(jsec), kjmin(jsec), kjmax(jsec)
- END DO
-
- CLOSE(inum)
-
- END SUBROUTINE section_init
-
- SUBROUTINE bimg_writ( ksec)
- !!---------------------------------------------------------------------
- !! *** ROUTINE bimg_writ ***
- !!
- !! ** Purpose : Write output bimg files if required
- !!
- !! ** Method : Most of the variables are global
- !!
- !!----------------------------------------------------------------------
- INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
-
- INTEGER(KIND=4) :: ji, jk
- !!----------------------------------------------------------------------
- ! (along section, depth ) 2D variables
- cf_bimg=TRIM(csection(ksec))//'_trpdep.bimg'
- OPEN(numbimg,FILE=cf_bimg,FORM='UNFORMATTED')
- cldum=' 4 dimensions in this isopycnal file '
- WRITE(numbimg) cldum
-
- cldum=' 1: T ; 2: S ; 3: sigma ; 4: Velocity '
- WRITE(numbimg) cldum
-
- WRITE(cldum,'(a,4i5.4)') ' from '//TRIM(csection(ksec)), iimin,iimax,ijmin,ijmax
- WRITE(numbimg) cldum
-
- cldum=' file '//TRIM(cf_tfil)
- WRITE(numbimg) cldum
-
- WRITE(numbimg) npts,nk,1,1,4,0
- WRITE(numbimg) 1.,-float(nk),1.,1., 0.
- WRITE(numbimg) 0.
- WRITE(numbimg) 0.
-
- WRITE(numbimg) (( REAL(zt(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! temperature
- WRITE(numbimg) (( REAL(zs(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! salinity
- WRITE(numbimg) (( REAL(dsig(ji,jk)), ji=1,npts), jk=nk,1,-1 ) ! density
- WRITE(numbimg) (( REAL(zu(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! normal velocity
- CLOSE(numbimg)
-
- ! (along section, sigma ) 2D variables
- cf_bimg=TRIM(csection(ksec))//'_trpsig.bimg'
- OPEN(numbimg,FILE=cf_bimg,FORM='UNFORMATTED')
- cldum=' 3 dimensions in this isopycnal file '
- WRITE(numbimg) cldum
- cldum=' 1: hiso ; 2: bin trp ; 3: cumulated trp '
- WRITE(numbimg) cldum
- WRITE(cldum,'(a,4i5.4)') ' from '//TRIM(csection(ksec)), iimin,iimax,ijmin,ijmax
- WRITE(numbimg) cldum
- cldum=' file '//TRIM(cf_tfil)
- WRITE(numbimg) cldum
- WRITE(numbimg) npts,nbins,1,1,3,0
- WRITE(numbimg) 1.,-REAL(dsigma_lev(nbins)),1.,REAL(dltsig), 0.
- WRITE(numbimg) 0.
- WRITE(numbimg) 0.
- WRITE(numbimg) (( REAL(dhiso(ji,jiso) ), ji=1,npts), jiso=nbins,1,-1) ! isopyc depth
- WRITE(numbimg) (( REAL(dwtrpbin(ji,jiso))/1.e6, ji=1,npts), jiso=nbins,1,-1) ! binned transport
- WRITE(numbimg) (( REAL(dwtrp(ji,jiso) )/1.e6, ji=1,npts), jiso=nbins,1,-1) ! cumulated transport
- CLOSE(numbimg)
-
- END SUBROUTINE bimg_writ
-
- SUBROUTINE cdf_writ( ksec)
- !!---------------------------------------------------------------------
- !! *** ROUTINE cdf_writ ***
- !!
- !! ** Purpose : Write output cdf files if required
- !!
- !! ** Method : Most of the variables are global
- !!
- !!----------------------------------------------------------------------
- INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
-
- INTEGER(KIND=4) :: ji, jk
- INTEGER(KIND=4) :: ivar
- INTEGER(KIND=4) :: icout
- INTEGER(KIND=4), DIMENSION(4) :: ipk, id_varout
-
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zdum
- TYPE(variable), DIMENSION(4) :: sl_typvar
- !!----------------------------------------------------------------------
- ALLOCATE ( zdum(npts,1))
- ! (along section, depth ) 2D variables
- cf_nc=TRIM(csection(ksec))//'_secdep.nc'
- ! define variables
- ipk(:)=nk
- sl_typvar%rmissing_value = 0.
- sl_typvar%rmissing_value = 0.
- sl_typvar%scale_factor = 1.
- sl_typvar%add_offset = 0.
- sl_typvar%savelog10 = 0.
- sl_typvar%iwght = iweight
- sl_typvar%conline_operation = 'N/A'
- sl_typvar%caxis = 'XZT'
-
- ivar=1
- sl_typvar(ivar)%cname = 'temperature'
- sl_typvar(ivar)%cunits = 'Celsius'
- sl_typvar(ivar)%valid_min = -2.
- sl_typvar(ivar)%valid_max = 45.
- sl_typvar(ivar)%clong_name = 'Potential_temperature'
- sl_typvar(ivar)%cshort_name = 'temperature'
-
- ivar=ivar+1
- sl_typvar(ivar)%cname = 'salinity'
- sl_typvar(ivar)%cunits = 'PSU'
- sl_typvar(ivar)%valid_min = 0.
- sl_typvar(ivar)%valid_max = 45.
- sl_typvar(ivar)%clong_name = 'Salinity'
- sl_typvar(ivar)%cshort_name = 'salinity'
-
- ivar=ivar+1
- sl_typvar(ivar)%cname = 'density'
- sl_typvar(ivar)%cunits = 'kg/m3 -1000'
- sl_typvar(ivar)%valid_min = 0.
- sl_typvar(ivar)%valid_max = 45.
- sl_typvar(ivar)%clong_name = 'potential_density'
- sl_typvar(ivar)%cshort_name = 'density'
-
- ivar=ivar+1
- sl_typvar(ivar)%cname = 'velocity'
- sl_typvar(ivar)%cunits = 'm/s'
- sl_typvar(ivar)%valid_min = -3.
- sl_typvar(ivar)%valid_max = 3.
- sl_typvar(ivar)%clong_name = 'Normal_velocity'
- sl_typvar(ivar)%cshort_name = 'velocity'
-
- icout = create (cf_nc, 'none', npts, 1, nk, cdep=cn_vdeptht )
- ierr = createvar (icout, sl_typvar, ivar, ipk, id_varout, cdglobal=TRIM(cglobal) )
- ierr = putheadervar(icout, cf_tfil, npts, 1, nk, &
- & pnavlon=rlonlat, pnavlat=rlonlat, cdep=cn_vdeptht )
-
-! tim = getvar1d(cf_tfil, cn_vtimec, 1 )
-! ierr = putvar1d(icout, tim, 1, 'T')
-
- DO jk = 1, nk
- zdum(:,1)=zt(:,jk) ; ierr = putvar ( icout, id_varout(1), zdum, jk, npts, 1 )
- zdum(:,1)=zs(:,jk) ; ierr = putvar ( icout, id_varout(2), zdum, jk, npts, 1 )
- zdum(:,1)=dsig(:,jk) ; ierr = putvar ( icout, id_varout(3), zdum, jk, npts, 1 )
- zdum(:,1)=zu(:,jk) ; ierr = putvar ( icout, id_varout(4), zdum, jk, npts, 1 )
- END DO
-
- ierr = closeout(icout)
-
- ! (along section, sigma ) 2D variables
- cf_nc=TRIM(csection(ksec))//'_secsig.nc'
- ! define variables
- ipk(:)=nbins
- sl_typvar%rmissing_value = 99999.
- sl_typvar%rmissing_value = 99999.
- sl_typvar%scale_factor = 1.
- sl_typvar%add_offset = 0.
- sl_typvar%savelog10 = 0.
- sl_typvar%iwght = iweight
- sl_typvar%conline_operation = 'N/A'
- sl_typvar%caxis = 'XST'
-
- ivar=1
- ipk(ivar)=nbins-1
- sl_typvar(ivar)%cname = 'isodep'
- sl_typvar(ivar)%cunits = 'm'
- sl_typvar(ivar)%valid_min = 0.
- sl_typvar(ivar)%valid_max = 6000.
- sl_typvar(ivar)%clong_name = 'isopycnal_depth'
- sl_typvar(ivar)%cshort_name = 'isodep'
-
- ivar=ivar+1
- sl_typvar(ivar)%cname = 'bintrp'
- sl_typvar(ivar)%cunits = 'SV'
- sl_typvar(ivar)%valid_min = -5.
- sl_typvar(ivar)%valid_max = 5.
- sl_typvar(ivar)%clong_name = 'Binned_transport'
- sl_typvar(ivar)%cshort_name = 'bintrp'
-
- ivar=ivar+1
- sl_typvar(ivar)%cname = 'sumtrp'
- sl_typvar(ivar)%cunits = 'SV'
- sl_typvar(ivar)%valid_min = -20.
- sl_typvar(ivar)%valid_max = 20.
- sl_typvar(ivar)%clong_name = 'cumulated_transport'
- sl_typvar(ivar)%cshort_name = 'sumtrp'
-
- icout = create (cf_nc, 'none', npts, 1, nbins, cdep='levels' )
- ierr = createvar (icout, sl_typvar, ivar, ipk, id_varout, cdglobal=TRIM(cglobal) )
- ierr = putheadervar(icout, cf_tfil, npts, 1, nbins, &
- & pnavlon=rlonlat, pnavlat=rlonlat, pdep=REAL(dsigma_lev), cdep='levels' )
-
- PRINT *, 'NBINS', nbins, npts
- DO jk = 1, nbins-1
- zdum(:,1)=dhiso (:,jk) ; ierr = putvar ( icout, id_varout(1), zdum, jk, npts, 1 )
- END DO
- DO jk = 1, nbins
- zdum(:,1)=dwtrpbin(:,jk)/1.e6 ; ierr = putvar ( icout, id_varout(2), zdum, jk, npts, 1 )
- zdum(:,1)=dwtrp (:,jk)/1.e6 ; ierr = putvar ( icout, id_varout(3), zdum, jk, npts, 1 )
- END DO
- ierr = closeout(icout)
-
- DEALLOCATE ( zdum )
-
- END SUBROUTINE cdf_writ
-
- SUBROUTINE print_out(ksec)
- !!---------------------------------------------------------------------
- !! *** ROUTINE print_out ***
- !!
- !! ** Purpose : Print results on standard output
- !!
- !! ** Method : Most of the variables are global and already known
- !!
- !!----------------------------------------------------------------------
- INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
-
- INTEGER(KIND=4) :: ji, jk, jiso, jbin
- !!----------------------------------------------------------------------
- WRITE(cfmt_9000,'(a,i3,a)') '(i7, ',npts,'f8.3)'
- WRITE(cfmt_9001,'(a,i3,a)') '(i7, ',npts,'f8.0)'
- WRITE(cfmt_9002,'(a,i3,a)') '(f7.3,',npts,'f8.0)'
- WRITE(cfmt_9003,'(a,i3,a)') '(f7.3,',npts,'f8.3)'
- PRINT *,' T (deg C)'
- DO jk=1,nk
- PRINT cfmt_9000, jk, (zt(ji,jk),ji=1,npts)
- END DO
-
- PRINT *,' S (PSU)'
- DO jk=1,nk
- PRINT cfmt_9000, jk, (zs(ji,jk),ji=1,npts)
- END DO
-
- PRINT *,' SIG (kg/m3 - 1000 )'
- DO jk=1,nk
- PRINT cfmt_9000, jk, (dsig(ji,jk),ji=1,npts)
- END DO
-
- PRINT *,' VELOCITY (cm/s ) '
- DO jk=1,nk
- PRINT cfmt_9000, jk, (zu(ji,jk)*100,ji=1,npts)
- END DO
-
- PRINT *,' GDEPU (m) '
- DO jk=1,nk
- PRINT cfmt_9001,jk, (ddepu(ji,jk)*zmask(ji,jk),ji=1,npts)
- END DO
-
- PRINT *, 'E3 (m)'
- DO jk=1,nk
- PRINT cfmt_9001,jk, (de3(ji,jk)*zmask(ji,jk),ji=1,npts)
- END DO
-
- PRINT *,' DEP ISO ( m )'
- DO jiso =1, nbins+1
- PRINT cfmt_9002, dsigma_lev(jiso),(dhiso(ji,jiso),ji=1,npts)
- END DO
-
- PRINT *,' TRP SURF --> ISO (SV)'
- DO jiso =1, nbins+1
- PRINT cfmt_9003, dsigma_lev(jiso),(dwtrp(ji,jiso)/1.d6,ji=1,npts)
- END DO
-
- PRINT *,' TRP bins (SV)'
- DO jbin =1, nbins
- PRINT cfmt_9003, dsigma_lev(jbin),(dwtrpbin(ji,jbin)/1.d6,ji=1,npts), dtrpbin(ksec,jbin)/1.d6
- END DO
-
-
-
- END SUBROUTINE print_out
+ SUBROUTINE section_init(cdfile, cdsection, cdvarname, cdlongname, kimin, kimax, kjmin, kjmax, knumber)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE section_init ***
+ !!
+ !! ** Purpose : Read input ASCII file that defines section names and limit of
+ !! sections.
+ !!
+ !! ** Method : At fisrt call only return the number of sections for further
+ !! allocation.
+ !!
+ !!----------------------------------------------------------------------
+ CHARACTER(LEN=*), INTENT(in ) :: cdfile
+ CHARACTER(LEN=256), DIMENSION(knumber), INTENT(out ) :: cdsection
+ CHARACTER(LEN=256), DIMENSION(knumber), INTENT(out ) :: cdvarname
+ CHARACTER(LEN=256), DIMENSION(knumber), INTENT(out ) :: cdlongname
+ INTEGER(KIND=4), INTENT(inout) :: knumber
+ INTEGER(KIND=4), DIMENSION(knumber), INTENT(out ) :: kimin, kimax, kjmin, kjmax
+
+ ! Local variables
+ INTEGER(KIND=4) :: jsec
+ INTEGER(KIND=4) :: ii, inum=10
+ INTEGER(KIND=4) :: ipos
+ CHARACTER(LEN=256) :: cline
+ CHARACTER(LEN=80), DIMENSION(3) :: cldum
+ LOGICAL :: llfirst
+ !!----------------------------------------------------------------------
+ llfirst=.FALSE.
+ IF ( knumber == 0 ) llfirst=.TRUE.
+
+ OPEN(inum, FILE=cdfile)
+ REWIND(inum)
+ ii = 0
+
+ ! read the file just to count the number of sections
+ DO
+ READ(inum,'(a)') cline
+ IF (INDEX(cline,'EOF') == 0 ) THEN
+ READ(inum,*) ! skip one line
+ ii = ii + 1
+ ELSE
+ EXIT
+ ENDIF
+ END DO
+
+ knumber=ii
+ IF ( llfirst ) RETURN
+
+ REWIND(inum)
+ DO jsec=1,knumber
+ READ(inum,'(a)') cline
+ ii = 0
+ cldum(:) = 'none'
+ ipos = index(cline,' ')
+ DO WHILE ( ipos > 1 )
+ ii = ii + 1
+ cldum(ii) = cline(1:ipos - 1 )
+ cline = TRIM ( cline(ipos+1:) )
+ ipos = index( cline,' ' )
+ IF ( ii >= 3 ) EXIT
+ END DO
+ cdsection(jsec) = TRIM(cldum(1) )
+ cdvarname(jsec) = TRIM(cldum(2) )
+ cdlongname(jsec) = TRIM(cldum(3) )
+ READ(inum,* ) kimin(jsec), kimax(jsec), kjmin(jsec), kjmax(jsec)
+ END DO
+
+ CLOSE(inum)
+
+ END SUBROUTINE section_init
+
+ SUBROUTINE bimg_writ( ksec)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE bimg_writ ***
+ !!
+ !! ** Purpose : Write output bimg files if required
+ !!
+ !! ** Method : Most of the variables are global
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
+
+ INTEGER(KIND=4) :: ji, jk
+ !!----------------------------------------------------------------------
+ ! (along section, depth ) 2D variables
+ cf_bimg=TRIM(csection(ksec))//'_trpdep.bimg'
+ OPEN(numbimg,FILE=cf_bimg,FORM='UNFORMATTED')
+ cldum=' 4 dimensions in this isopycnal file '
+ WRITE(numbimg) cldum
+
+ cldum=' 1: T ; 2: S ; 3: sigma ; 4: Velocity '
+ WRITE(numbimg) cldum
+
+ WRITE(cldum,'(a,4i5.4)') ' from '//TRIM(csection(ksec)), iimin,iimax,ijmin,ijmax
+ WRITE(numbimg) cldum
+
+ cldum=' file '//TRIM(cf_tfil)
+ WRITE(numbimg) cldum
+
+ WRITE(numbimg) npts,nk,1,1,4,0
+ WRITE(numbimg) 1.,-float(nk),1.,1., 0.
+ WRITE(numbimg) 0.
+ WRITE(numbimg) 0.
+
+ WRITE(numbimg) (( REAL(zt(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! temperature
+ WRITE(numbimg) (( REAL(zs(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! salinity
+ WRITE(numbimg) (( REAL(dsig(ji,jk)), ji=1,npts), jk=nk,1,-1 ) ! density
+ WRITE(numbimg) (( REAL(zu(ji,jk) ), ji=1,npts), jk=nk,1,-1 ) ! normal velocity
+ CLOSE(numbimg)
+
+ ! (along section, sigma ) 2D variables
+ cf_bimg=TRIM(csection(ksec))//'_trpsig.bimg'
+ OPEN(numbimg,FILE=cf_bimg,FORM='UNFORMATTED')
+ cldum=' 3 dimensions in this isopycnal file '
+ WRITE(numbimg) cldum
+ cldum=' 1: hiso ; 2: bin trp ; 3: cumulated trp '
+ WRITE(numbimg) cldum
+ WRITE(cldum,'(a,4i5.4)') ' from '//TRIM(csection(ksec)), iimin,iimax,ijmin,ijmax
+ WRITE(numbimg) cldum
+ cldum=' file '//TRIM(cf_tfil)
+ WRITE(numbimg) cldum
+ WRITE(numbimg) npts,nbins,1,1,3,0
+ WRITE(numbimg) 1.,-REAL(dsigma_lev(nbins)),1.,REAL(dltsig), 0.
+ WRITE(numbimg) 0.
+ WRITE(numbimg) 0.
+ WRITE(numbimg) (( REAL(dhiso(ji,jiso) ), ji=1,npts), jiso=nbins,1,-1) ! isopyc depth
+ WRITE(numbimg) (( REAL(dwtrpbin(ji,jiso))/1.e6, ji=1,npts), jiso=nbins,1,-1) ! binned transport
+ WRITE(numbimg) (( REAL(dwtrp(ji,jiso) )/1.e6, ji=1,npts), jiso=nbins,1,-1) ! cumulated transport
+ CLOSE(numbimg)
+
+ END SUBROUTINE bimg_writ
+
+ SUBROUTINE cdf_writ( ksec)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE cdf_writ ***
+ !!
+ !! ** Purpose : Write output cdf files if required
+ !!
+ !! ** Method : Most of the variables are global
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
+
+ INTEGER(KIND=4) :: ji, jk
+ INTEGER(KIND=4) :: ivar
+ INTEGER(KIND=4) :: icout
+ INTEGER(KIND=4), DIMENSION(4) :: ipk, id_varout
+
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zdum
+ TYPE(variable), DIMENSION(4) :: sl_typvar
+ CHARACTER(LEN=255) :: csuffixvarnam
+ CHARACTER(LEN=255) :: cprefixlongnam
+ !!----------------------------------------------------------------------
+ IF ( cvarname(ksec) /= 'none' ) THEN
+ csuffixvarnam = '_'//TRIM(cvarname(ksec))
+ ELSE
+ csuffixvarnam = ''
+ ENDIF
+ IF ( clongname(ksec) /= 'none' ) THEN
+ cprefixlongnam = TRIM(clongname(ksec))//'_'
+ ELSE
+ cprefixlongnam = ''
+ ENDIF
+
+ ALLOCATE ( zdum(npts,1))
+ ! (along section, depth ) 2D variables
+ cf_nc=TRIM(csection(ksec))//'_secdep.nc'
+ ! define variables
+ ipk(:)=nk
+ sl_typvar%rmissing_value = 0.
+ sl_typvar%rmissing_value = 0.
+ sl_typvar%scale_factor = 1.
+ sl_typvar%add_offset = 0.
+ sl_typvar%savelog10 = 0.
+ sl_typvar%iwght = iweight
+ sl_typvar%conline_operation = 'N/A'
+ sl_typvar%caxis = 'XZT'
+
+ ivar=1
+ sl_typvar(ivar)%cname = 'temperature'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'Celsius'
+ sl_typvar(ivar)%valid_min = -2.
+ sl_typvar(ivar)%valid_max = 45.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Potential_temperature'
+ sl_typvar(ivar)%cshort_name = 'temperature'
+
+ ivar=ivar+1
+ sl_typvar(ivar)%cname = 'salinity'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'PSU'
+ sl_typvar(ivar)%valid_min = 0.
+ sl_typvar(ivar)%valid_max = 45.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Salinity'
+ sl_typvar(ivar)%cshort_name = 'salinity'
+
+ ivar=ivar+1
+ sl_typvar(ivar)%cname = 'density'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'kg/m3 -1000'
+ sl_typvar(ivar)%valid_min = 0.
+ sl_typvar(ivar)%valid_max = 45.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'potential_density'
+ sl_typvar(ivar)%cshort_name = 'density'
+
+ ivar=ivar+1
+ sl_typvar(ivar)%cname = 'velocity'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'm/s'
+ sl_typvar(ivar)%valid_min = -3.
+ sl_typvar(ivar)%valid_max = 3.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Normal_velocity'
+ sl_typvar(ivar)%cshort_name = 'velocity'
+
+ icout = create (cf_nc, 'none', npts, 1, nk, cdep=cn_vdeptht )
+ ierr = createvar (icout, sl_typvar, ivar, ipk, id_varout, cdglobal=TRIM(cglobal) )
+ ierr = putheadervar(icout, cf_tfil, npts, 1, nk, &
+ & pnavlon=rlonlat, pnavlat=rlonlat, cdep=cn_vdeptht )
+
+ ! tim = getvar1d(cf_tfil, cn_vtimec, 1 )
+ ! ierr = putvar1d(icout, tim, 1, 'T')
+
+ DO jk = 1, nk
+ zdum(:,1)=zt(:,jk) ; ierr = putvar ( icout, id_varout(1), zdum, jk, npts, 1 )
+ zdum(:,1)=zs(:,jk) ; ierr = putvar ( icout, id_varout(2), zdum, jk, npts, 1 )
+ zdum(:,1)=dsig(:,jk) ; ierr = putvar ( icout, id_varout(3), zdum, jk, npts, 1 )
+ zdum(:,1)=zu(:,jk) ; ierr = putvar ( icout, id_varout(4), zdum, jk, npts, 1 )
+ END DO
+
+ ierr = closeout(icout)
+
+ ! (along section, sigma ) 2D variables
+ cf_nc=TRIM(csection(ksec))//'_secsig.nc'
+ ! define variables
+ ipk(:)=nbins
+ sl_typvar%rmissing_value = 99999.
+ sl_typvar%rmissing_value = 99999.
+ sl_typvar%scale_factor = 1.
+ sl_typvar%add_offset = 0.
+ sl_typvar%savelog10 = 0.
+ sl_typvar%iwght = iweight
+ sl_typvar%conline_operation = 'N/A'
+ sl_typvar%caxis = 'XST'
+
+ ivar=1
+ ipk(ivar)=nbins-1
+ sl_typvar(ivar)%cname = 'isodep'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'm'
+ sl_typvar(ivar)%valid_min = 0.
+ sl_typvar(ivar)%valid_max = 6000.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'isopycnal_depth'
+ sl_typvar(ivar)%cshort_name = 'isodep'
+
+ ivar=ivar+1
+ sl_typvar(ivar)%cname = 'bintrp'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'SV'
+ sl_typvar(ivar)%valid_min = -5.
+ sl_typvar(ivar)%valid_max = 5.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Binned_transport'
+ sl_typvar(ivar)%cshort_name = 'bintrp'
+
+ ivar=ivar+1
+ sl_typvar(ivar)%cname = 'sumtrp'//TRIM(csuffixvarnam)
+ sl_typvar(ivar)%cunits = 'SV'
+ sl_typvar(ivar)%valid_min = -20.
+ sl_typvar(ivar)%valid_max = 20.
+ sl_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'cumulated_transport'
+ sl_typvar(ivar)%cshort_name = 'sumtrp'
+
+ icout = create (cf_nc, 'none', npts, 1, nbins, cdep='levels' )
+ ierr = createvar (icout, sl_typvar, ivar, ipk, id_varout, cdglobal=TRIM(cglobal) )
+ ierr = putheadervar(icout, cf_tfil, npts, 1, nbins, &
+ & pnavlon=rlonlat, pnavlat=rlonlat, pdep=REAL(dsigma_lev), cdep='levels' )
+
+ PRINT *, 'NBINS', nbins, npts
+ DO jk = 1, nbins-1
+ zdum(:,1)=dhiso (:,jk) ; ierr = putvar ( icout, id_varout(1), zdum, jk, npts, 1 )
+ END DO
+ DO jk = 1, nbins
+ zdum(:,1)=dwtrpbin(:,jk)/1.e6 ; ierr = putvar ( icout, id_varout(2), zdum, jk, npts, 1 )
+ zdum(:,1)=dwtrp (:,jk)/1.e6 ; ierr = putvar ( icout, id_varout(3), zdum, jk, npts, 1 )
+ END DO
+ ierr = closeout(icout)
+
+ DEALLOCATE ( zdum )
+
+ END SUBROUTINE cdf_writ
+
+ SUBROUTINE print_out(ksec)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE print_out ***
+ !!
+ !! ** Purpose : Print results on standard output
+ !!
+ !! ** Method : Most of the variables are global and already known
+ !!
+ !!----------------------------------------------------------------------
+ INTEGER(KIND=4), INTENT(in) :: ksec ! number of the section
+
+ INTEGER(KIND=4) :: ji, jk, jiso, jbin
+ !!----------------------------------------------------------------------
+ WRITE(cfmt_9000,'(a,i3,a)') '(i7, ',npts,'f8.3)'
+ WRITE(cfmt_9001,'(a,i3,a)') '(i7, ',npts,'f8.0)'
+ WRITE(cfmt_9002,'(a,i3,a)') '(f7.3,',npts,'f8.0)'
+ WRITE(cfmt_9003,'(a,i3,a)') '(f7.3,',npts,'f8.3)'
+ PRINT *,' T (deg C)'
+ DO jk=1,nk
+ PRINT cfmt_9000, jk, (zt(ji,jk),ji=1,npts)
+ END DO
+
+ PRINT *,' S (PSU)'
+ DO jk=1,nk
+ PRINT cfmt_9000, jk, (zs(ji,jk),ji=1,npts)
+ END DO
+
+ PRINT *,' SIG (kg/m3 - 1000 )'
+ DO jk=1,nk
+ PRINT cfmt_9000, jk, (dsig(ji,jk),ji=1,npts)
+ END DO
+
+ PRINT *,' VELOCITY (cm/s ) '
+ DO jk=1,nk
+ PRINT cfmt_9000, jk, (zu(ji,jk)*100,ji=1,npts)
+ END DO
+
+ PRINT *,' GDEPU (m) '
+ DO jk=1,nk
+ PRINT cfmt_9001,jk, (ddepu(ji,jk)*zmask(ji,jk),ji=1,npts)
+ END DO
+
+ PRINT *, 'E3 (m)'
+ DO jk=1,nk
+ PRINT cfmt_9001,jk, (de3(ji,jk)*zmask(ji,jk),ji=1,npts)
+ END DO
+
+ PRINT *,' DEP ISO ( m )'
+ DO jiso =1, nbins+1
+ PRINT cfmt_9002, dsigma_lev(jiso),(dhiso(ji,jiso),ji=1,npts)
+ END DO
+
+ PRINT *,' TRP SURF --> ISO (SV)'
+ DO jiso =1, nbins+1
+ PRINT cfmt_9003, dsigma_lev(jiso),(dwtrp(ji,jiso)/1.d6,ji=1,npts)
+ END DO
+
+ PRINT *,' TRP bins (SV)'
+ DO jbin =1, nbins
+ PRINT cfmt_9003, dsigma_lev(jbin),(dwtrpbin(ji,jbin)/1.d6,ji=1,npts), dtrpbin(ksec,jbin)/1.d6
+ END DO
+
+
+
+ END SUBROUTINE print_out
END PROGRAM cdfsigtrp
diff --git a/cdftransport.f90 b/cdftransport.f90
index 8a129f0..de36953 100644
--- a/cdftransport.f90
+++ b/cdftransport.f90
@@ -1,1120 +1,1176 @@
PROGRAM cdftransport
- !!======================================================================
- !! *** PROGRAM cdftransport ***
- !!=====================================================================
- !! ** Purpose : Compute Transports across a section.
- !! By default, mass (Sv) and heat(PW)/salt(kT/s) transports
- !! are computed unless -noheat option is used (mass
- !! transport only).
- !!
- !! ** Method : The begining and end point of the section are given in
- !! term of F-points index. A broken line joining successive
- !! F-points is defined between the begining and end point
- !! of the section. Therefore each segment between F-points
- !! is either a zonal or meridional segment corresponding to
- !! V or U velocity component. Doing so, the volume conservation
- !! is ensured as velocities are not interpolated, and stay
- !! on the native model grid.
- !! The section name and the begin/end point of a section are
- !! read from standard input, till 'EOF' is given as section
- !! name. This make possible to give a bunch of sections in
- !! an ASCII files and use the < redirection.
- !! SIGN CONVENTION : The transport is positive when the flow cross
- !! the section to the right, negative otherwise. This depends
- !! on the sense the section is described. With this convention
- !! The algebric sum of transports accross sections forming a
- !! closed area is 0.
- !! OPTIONS :
- !! -full : full step case
- !! -noheat : only mass transport is computed.
- !! -time : specify the time frame to be used
- !! -zlimit : transports can be computed in different depth layers
- !! defined by their depth limit
- !! REQUIREMENT :
- !! mesh-mask file are required in the current directory.
- !!
- !!
- !! History : 2.1 : 01/2005 : J.M. Molines : Original code
- !! 2.1 : 07/2009 : R. Dussin : add cdf output
- !! 2.1 : 01/2010 : M.A. Balmaseda : Change integration signs
- !! so that the transport across a segment is
- !! independent of the chosen trajectory.
- !! 3.0 : 04/2011 : J.M. Molines : Doctor norm + Lic.
- !!----------------------------------------------------------------------
- !!----------------------------------------------------------------------
- !! routines : description
- !! interm_pt : choose intermediate points on a broken line.
- !!----------------------------------------------------------------------
- USE cdfio
- USE modcdfnames
- USE modutils ! for global attribute
- !!----------------------------------------------------------------------
- !! CDFTOOLS_3.0 , MEOM 2011
- !! $Id$
- !! Copyright (c) 2011, J.-M. Molines
- !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
- !!----------------------------------------------------------------------
- IMPLICIT NONE
-
- INTEGER(KIND=4) :: jclass, jseg ! dummy loop index
- INTEGER(KIND=4) :: ji, jj, jk ! dummy loop index
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: imeter ! limit beetween depth level, in m (nclass -1)
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ilev0, ilev1 ! limit in levels (nclass)
- INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_varout ! Netcdf output
- INTEGER(KIND=4) :: ncout, ierr ! for netcdf output
- INTEGER(KIND=4) :: nvarout=12 ! number of values to write in cdf output
- INTEGER(KIND=4) :: ivtrp ! var index of volume transport (barotrope)
- INTEGER(KIND=4) :: iptrp ! var index of volume transport (barotrope)
- INTEGER(KIND=4) :: imtrp ! var index of volume transport (barotrope)
- INTEGER(KIND=4) :: ihtrp ! var index of heat transport (barotrope)
- INTEGER(KIND=4) :: istrp ! var index of sal transport (barotrope)
- INTEGER(KIND=4) :: ivtrpcl ! var index of volume transport (p. class)
- INTEGER(KIND=4) :: iptrpcl ! var index of volume transport (p. class)
- INTEGER(KIND=4) :: imtrpcl ! var index of volume transport (p. class)
- INTEGER(KIND=4) :: ihtrpcl ! var index of heat transport (p. class)
- INTEGER(KIND=4) :: istrpcl ! var index of sal transport (p. class)
- INTEGER(KIND=4) :: ilonmin ! var index of starting section longitude
- INTEGER(KIND=4) :: ilonmax ! var index of ending section longitude
- INTEGER(KIND=4) :: ilatmin ! var index of starting section latitude
- INTEGER(KIND=4) :: ilatmax ! var index of ending section latitude
- INTEGER(KIND=4) :: itop ! var index of top depth class
- INTEGER(KIND=4) :: ibot ! var index of bottom depth class
- INTEGER(KIND=4) :: ikx=1, iky=1 ! dims of netcdf output file
- INTEGER(KIND=4) :: numout = 10 ! logical unit for output file (overall)
- INTEGER(KIND=4) :: numvtrp = 11 ! logical unit for volume transport file
- INTEGER(KIND=4) :: numhtrp = 12 ! logical unit for heat transport file
- INTEGER(KIND=4) :: numstrp = 14 ! logical unit for salt trp file
- INTEGER(KIND=4) :: nclass ! number of depth class
- INTEGER(KIND=4) :: narg, iargc ! command line
- INTEGER(KIND=4) :: ijarg, nxtarg ! " "
- INTEGER(KIND=4) :: npiglo, npjglo ! size of the domain
- INTEGER(KIND=4) :: npk, npt ! size of the domain
- INTEGER(KIND=4) :: iimin, iimax ! i-limit of the section
- INTEGER(KIND=4) :: ijmin, ijmax ! j-limit of the section
- INTEGER(KIND=4) :: ivar, itime ! working integer
- INTEGER(KIND=4) :: ii, ij, ik ! working integer
- INTEGER(KIND=4), PARAMETER :: jpseg=10000 ! used for broken line algorithm
- INTEGER(KIND=4) :: ii0, ij0 ! " " "
- INTEGER(KIND=4) :: ii1, ij1 ! " " "
- INTEGER(KIND=4) :: iitmp, ijtmp ! " " "
- INTEGER(KIND=4) :: np, nn ! segment counters,
- INTEGER(KIND=4) :: iist, ijst ! local point offset for velocity
- INTEGER(KIND=4) :: norm_u, norm_v ! normalization factor (sign of normal transport)
- INTEGER(KIND=4) :: idirx, idiry ! sense of description of the section
-
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e1v, e2u ! horizontal metric
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e3u, e3v ! vertical metric
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamf ! longitudes of F points
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: gphif ! latitudes of F points
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zu, zut, zus ! Zonal velocities and uT uS
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zv, zvt, zvs ! Meridional velocities and uT uS
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rdum ! dummy (1x1) array for ncdf output
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zuobc, zvobc ! arrays for OBC files (vertical slice)
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tim ! time counter
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdepw ! depth at layer interface
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: e31d ! vertical metric in case of full step
- REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: rclass ! vertical metric in case of full step
- REAL(KIND=4), DIMENSION(2) :: gla, gphi ! lon/lat of the begining/end of section (f point)
- REAL(KIND=4), DIMENSION(jpseg) :: rxx, ryy ! working variables
- 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)
-
- ! at every model point
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwku, dwkv ! volume transport at each cell boundary
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkut, dwkvt ! heat transport at each cell boundary
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkus, dwkvs ! salt transport at each cell boundary
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkup, dwkvp ! volume transport in the positive direction
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkum, dwkvm ! volume transport in the negatibe direction
- REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpu, dtrpv ! volume transport integrated in depth class
- REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrput, dtrpvt ! heat transport integrated in depth class
- REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpus, dtrpvs ! salt transport integrated in depth class
- REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpup, dtrpvp ! volume transport integrated in depth class (positive)
- REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpum, dtrpvm ! volume transport integrated in depth class (negative)
- ! for a given section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsum ! volume transport by depth class across section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsump ! volume transport by depth class across section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsumm ! volume transport by depth class across section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dheatrpsum ! heat transport by depth class across section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsaltrpsum ! salt transport by depth class across section
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegcl ! over all leg volume transport by depth class
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegclp ! over all leg volume transport by depth class +
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegclm ! over all leg volume transport by depth class -
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dheatallegcl ! over all leg heat transport by depth class
- REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsaltallegcl ! over all leg salt transport by depth class
- REAL(KIND=8), DIMENSION(jpseg) :: dvoltrp ! volume transport across each segment of a section
- REAL(KIND=8), DIMENSION(jpseg) :: dvoltrpp ! volume transport across each segment of a section
- REAL(KIND=8), DIMENSION(jpseg) :: dvoltrpm ! volume transport across each segment of a section
- REAL(KIND=8), DIMENSION(jpseg) :: dheatrp ! heat transport across each segment of a section
- REAL(KIND=8), DIMENSION(jpseg) :: dsaltrp ! salt transport across each segment of a section
- REAL(KIND=8) :: dvoltrpbrtp ! volume transport integrated over the whole depth
- REAL(KIND=8) :: dvoltrpbrtpp ! volume transport integrated over the whole depth
- REAL(KIND=8) :: dvoltrpbrtpm ! volume transport integrated over the whole depth
- REAL(KIND=8) :: dheatrpbrtp ! heat transport integrated over the whole depth
- REAL(KIND=8) :: dsaltrpbrtp ! salt transport integrated over the whole depth
- REAL(KIND=8) :: dvolalleg ! over all leg sum of volume transport
- REAL(KIND=8) :: dvolallegp ! over all leg sum of volume transport +
- REAL(KIND=8) :: dvolallegm ! over all leg sum of volume transport -
- REAL(KIND=8) :: dheatalleg ! over all leg sum of heat transport
- REAL(KIND=8) :: dsaltalleg ! over all leg sum of salt transport
-
- COMPLEX, DIMENSION(jpseg) :: yypt ! array of points coordinates in a section
- COMPLEX :: yypti ! working point
-
- TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvar ! structure of output
-
- CHARACTER(LEN=256) :: cf_tfil ! VT file (in)
- CHARACTER(LEN=256) :: cf_ufil ! U file (in)
- CHARACTER(LEN=256) :: cf_vfil ! V file (in)
- CHARACTER(LEN=256) :: cf_out='section_trp.dat' ! output file name (ASCII)
- CHARACTER(LEN=256) :: cf_outnc ! output netcdf file
- CHARACTER(LEN=256) :: cf_vtrp='vtrp.txt' ! output volume transport file
- CHARACTER(LEN=256) :: cf_htrp='htrp.txt' ! output heat transport file
- CHARACTER(LEN=256) :: cf_strp='strp.txt' ! output salt transport file
- CHARACTER(LEN=256) :: csection ! section names
- CHARACTER(LEN=512) :: cglobal ! global attribute
- CHARACTER(LEN=256) :: cldum ! dummy char variable
-
- LOGICAL :: ltest = .FALSE. ! flag for test case
- LOGICAL :: lfull = .FALSE. ! flag for full step case
- 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
- LOGICAL :: lobc = .FALSE. ! flag for obc input files
- LOGICAL :: l_merid = .FALSE. ! flag for meridional obc
- LOGICAL :: l_zonal = .FALSE. ! flag for zonal obc
- !!----------------------------------------------------------------------
- CALL ReadCdfNames()
-
- narg= iargc()
- ! Print usage if no argument
- IF ( narg == 0 ) THEN
- PRINT *,' usage : cdftransport [-test u v ] [-noheat ] [-plus_minus ] [-obc]...'
- PRINT *,' ... [VT-file] U-file V-file [-full] |-time jt] ...'
- PRINT *,' ... [-time jt ] [-zlimit limits of level]'
- PRINT *,' '
- PRINT *,' PURPOSE :'
- PRINT *,' Compute the transports accross a section.'
- PRINT *,' The name of the section and the imin, imax, jmin, jmax for the section '
- PRINT *,' is read from the standard input. To finish the program use the key name'
- PRINT *,' ''EOF'' for the section name.'
- PRINT *,' OBC U,V files can be used if -obc option is specified.'
- PRINT *,' '
- PRINT *,' ARGUMENTS :'
- PRINT *,' [VT-file ] : netcdf file with mean values of vt, vs, ut, us for heat and'
- PRINT *,' salt transport. If options -noheat or -plus_minus are used'
- PRINT *,' this file name must be omitted.'
- PRINT *,' [U-file ] : netcdf file with the zonal velocity component.'
- PRINT *,' [V-file ] : netcdf file with the meridional velocity component.'
- PRINT *,' '
- PRINT *,' OPTIONS :'
- PRINT *,' [-test u v ]: use constant the u and v velocity components for sign '
- PRINT *,' test purpose.'
- PRINT *,' [-noheat ] : use when heat and salt transport are not requested.'
- PRINT *,' This option must come before the file names, and if used'
- PRINT *,' VT file must not be given.'
- PRINT *,' [ -plus_minus or -pm ] : separate positive and negative contribution to'
- PRINT *,' the volume transport. This option implicitly set -noheat,'
- PRINT *,' and must be used before the file names.'
- PRINT *,' [-obc ] : indicates that input files are obc files (vertical slices)'
- PRINT *,' Take care that for this case, mesh files must be adapted.'
- PRINT *,' This option implicitly set -noheat, and must be used before'
- PRINT *,' the file names.'
- PRINT *,' [-full ] : use for full step configurations.'
- PRINT *,' [-time jt ]: compute transports for time index jt. Default is 1.'
- PRINT *,' [-zlimit list of depth] : Specify depths limits defining layers where the'
- PRINT *,' transports will be computed. If not used, the transports '
- PRINT *,' are computed for the whole water column. If used, this '
- PRINT *,' option must be the last on the command line.'
- PRINT *,' '
- PRINT *,' REQUIRED FILES :'
- PRINT *,' Files ',TRIM(cn_fhgr),', ',TRIM(cn_fzgr),' must be in the current directory.'
- PRINT *,' '
- PRINT *,' OUTPUT : '
- PRINT *,' - Standard output '
- PRINT *,' - ASCII file reflecting the standard output: section_trp.dat'
- PRINT *,' - ASCII files for volume, heat and salt transport: vtrp.txt, htrp.txt '
- PRINT *,' and strp.txt.'
- PRINT *,' - Netcdf files for each section. name of the file is buildt'
- PRINT *,' from section name.'
- PRINT *,' '
- PRINT *,' SEE ALSO :'
- PRINT *,' cdfsigtrp'
- PRINT *,' '
- STOP
- ENDIF
-
- itime = 1
- nclass = 1
- ijarg = 1
- CALL SetGlobalAtt(cglobal)
-
- ! Browse command line for arguments and/or options
- DO WHILE (ijarg <= narg )
- CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1
- SELECT CASE ( cldum )
- CASE ('-test ')
- CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) udum
- CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) vdum
- ltest = .TRUE.
-
- CASE ('-full' )
- lfull = .TRUE.
-
- CASE ('-noheat' ) ! it must be called before the list of files
- lheat = .FALSE.
-
- CASE ('-time' )
- CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) itime
-
- CASE ('-plus_minus', '-pm' )
- lpm = .TRUE.
- lheat = .FALSE.
-
- CASE ('-obc' )
- lobc = .TRUE.
- lheat = .FALSE.
-
- CASE ('-zlimit' ) ! this should be the last option on the line
- nxtarg = ijarg - 1
- nclass = narg - nxtarg + 1
- ALLOCATE ( imeter(nclass -1) ) ! if no zlimit option, this array is never used
- DO jclass =1, nclass -1
- CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) imeter(jclass)
- END DO
-
- CASE DEFAULT
- ijarg = ijarg -1 ! re-read argument in this case
- IF ( lheat) THEN
- CALL getarg (ijarg, cf_tfil) ; ijarg = ijarg + 1
- ENDIF
- CALL getarg (ijarg, cf_ufil) ; ijarg = ijarg + 1
- CALL getarg (ijarg, cf_vfil) ; ijarg = ijarg + 1
- END SELECT
- END DO
-
- ! checking if all required files are available
- lchk = lchk .OR. chkfile(cn_fzgr)
- lchk = lchk .OR. chkfile(cn_fhgr)
- IF ( ltest ) THEN
- ! OK
- ELSE
- lchk = lchk .OR. chkfile(cf_ufil)
- lchk = lchk .OR. chkfile(cf_vfil)
- IF (lheat) THEN
- lchk = lchk .OR. chkfile(cf_tfil)
- ENDIF
- ENDIF
- IF ( lchk ) STOP ! missing files
-
- ! adjust the number of output variables according to options
- IF ( nclass > 1 ) THEN
- IF ( lheat ) THEN
- nvarout = 12
- ELSE
- nvarout = 8
- ENDIF
- IF ( lpm ) nvarout=nvarout+4
- ELSE
- IF ( lheat ) THEN
- nvarout = 9
- ELSE
- nvarout = 7
- ENDIF
- IF ( lpm ) nvarout=nvarout+2
- ENDIF
-
- ALLOCATE ( ilev0(nclass), ilev1(nclass), rclass(nclass) )
- rclass=(/(jclass, jclass=1,nclass)/)
-
- npiglo = getdim (cf_ufil,cn_x)
- npjglo = getdim (cf_ufil,cn_y)
- npk = getdim (cf_ufil,cn_z)
- npt = getdim (cf_ufil,cn_t)
-
- PRINT *, 'npiglo =', npiglo
- PRINT *, 'npjglo =', npjglo
- PRINT *, 'npk =', npk
- PRINT *, 'npt =', npt
-
- IF ( lobc ) THEN ! if lobc false, l_merid and l_zonal are false (default)
- IF ( npiglo == 1 ) THEN
- l_merid=.TRUE.
- ALLOCATE (zuobc(npjglo,npk), zvobc(npjglo,npk) )
- PRINT *,' Meridional OBC'
- ENDIF
-
- IF ( npjglo == 1 ) THEN
- l_zonal=.TRUE.
- ALLOCATE (zuobc(npiglo,npk), zvobc(npiglo,npk) )
- PRINT *,' Zonal OBC'
- ENDIF
- ENDIF
-
- ALLOCATE ( e31d(npk) )
-
- ! define new variables for output
- ALLOCATE ( stypvar(nvarout), ipk(nvarout), id_varout(nvarout) )
- ALLOCATE ( rdum(1,1) )
-
- rdum(:,:)=0.e0
-
- stypvar%rmissing_value=99999.
- stypvar%scale_factor= 1.
- stypvar%add_offset= 0.
- stypvar%savelog10= 0.
- stypvar%conline_operation='N/A'
- stypvar%caxis='T'
-
- ivar = 1 ; ivtrp = ivar
- ipk(ivar) = 1
- stypvar(ivar)%cname = 'vtrp'
- stypvar(ivar)%cunits = 'Sverdrup'
- stypvar(ivar)%valid_min = -500.
- stypvar(ivar)%valid_max = 500.
- stypvar(ivar)%clong_name = 'Volume_Transport'
- stypvar(ivar)%cshort_name = 'vtrp'
-
- IF ( lpm ) THEN
- ivar = ivar + 1 ; iptrp = ivar ; imtrp = ivar+1
- ipk(ivar) = 1 ; ipk(ivar+1) = 1
- stypvar(ivar)%cname = 'ptrp' ; stypvar(ivar+1)%cname = 'mtrp'
- stypvar(ivar)%cunits = 'Sverdrup' ; stypvar(ivar+1)%cunits = 'Sverdrup'
- stypvar(ivar)%valid_min = -500. ; stypvar(ivar+1)%valid_min = -500.
- stypvar(ivar)%valid_max = 500. ; stypvar(ivar+1)%valid_max = 500.
- stypvar(ivar)%clong_name = 'Positive_volume_transport' ; stypvar(ivar+1)%clong_name = 'Negative_volume_transport'
- stypvar(ivar)%cshort_name = 'ptrp' ; stypvar(ivar+1)%cshort_name = 'mtrp'
- ivar = ivar + 1
- ENDIF
-
- IF ( lheat ) THEN
- ivar = ivar + 1 ; ihtrp = ivar ; istrp = ivar+1
- ipk(ivar) = 1 ; ipk(ivar+1) = 1
- stypvar(ivar)%cname = 'htrp' ; stypvar(ivar+1)%cname = 'strp'
- stypvar(ivar)%cunits = 'PW' ; stypvar(ivar+1)%cunits = 'kt/s'
- stypvar(ivar)%valid_min = -1000. ; stypvar(ivar+1)%valid_min = -1000.
- stypvar(ivar)%valid_max = 1000. ; stypvar(ivar+1)%valid_max = 1000.
- stypvar(ivar)%clong_name = 'Heat_Transport' ; stypvar(ivar+1)%clong_name = 'Salt_Transport'
- stypvar(ivar)%cshort_name = 'htrp' ; stypvar(ivar+1)%cshort_name = 'strp'
- ivar = ivar + 1
- ENDIF
-
- ivar = ivar + 1 ; ilonmin = ivar ; ilonmax = ivar+1
- ipk(ivar) = 1 ; ipk(ivar+1) = 1
- stypvar(ivar)%cname = 'lonmin' ; stypvar(ivar+1)%cname = 'lonmax'
- stypvar(ivar)%cunits = 'deg' ; stypvar(ivar+1)%cunits = 'deg'
- stypvar(ivar)%valid_min = -180. ; stypvar(ivar+1)%valid_min = -180.
- stypvar(ivar)%valid_max = 180. ; stypvar(ivar+1)%valid_max = 180.
- stypvar(ivar)%clong_name = 'begin_longitude' ; stypvar(ivar+1)%clong_name = 'end_longitude'
- stypvar(ivar)%cshort_name = 'lonmin' ; stypvar(ivar+1)%cshort_name = 'lonmax'
- ivar = ivar + 1
-
- ivar = ivar + 1 ; ilatmin = ivar ; ilatmax = ivar+1
- ipk(ivar) = 1 ; ipk(ivar+1) = 1
- stypvar(ivar)%cname = 'latmin' ; stypvar(ivar+1)%cname = 'latmax'
- stypvar(ivar)%cunits = 'deg' ; stypvar(ivar+1)%cunits = 'deg'
- stypvar(ivar)%valid_min = -90. ; stypvar(ivar+1)%valid_min = -90.
- stypvar(ivar)%valid_max = 90. ; stypvar(ivar+1)%valid_max = 90.
- stypvar(ivar)%clong_name = 'begin_latitude' ; stypvar(ivar+1)%clong_name = 'end_latitude'
- stypvar(ivar)%cshort_name = 'latmin' ; stypvar(ivar+1)%cshort_name = 'latmax'
- ivar = ivar + 1
-
- ivar = ivar + 1 ; itop = ivar ; ibot = ivar+1
- ipk(ivar) = nclass ; ipk(ivar+1) = nclass
- stypvar(ivar)%cname = 'top' ; stypvar(ivar+1)%cname = 'bottom'
- stypvar(ivar)%cunits = 'meters' ; stypvar(ivar+1)%cunits = 'meters'
- stypvar(ivar)%valid_min = 0. ; stypvar(ivar+1)%valid_min = 0.
- stypvar(ivar)%valid_max = 10000. ; stypvar(ivar+1)%valid_max = 10000.
- stypvar(ivar)%clong_name = 'class_min_depth' ; stypvar(ivar+1)%clong_name = 'class_max_depth'
- stypvar(ivar)%cshort_name = 'top' ; stypvar(ivar+1)%cshort_name = 'bottom'
- ivar = ivar + 1
-
- ivtrpcl = -1 ; ihtrpcl = -1 ; istrpcl = -1
- IF ( nclass > 1 ) THEN ! define additional variable for vertical profile of transport (per class)
- ivar = ivar + 1 ; ivtrpcl = ivar
- ipk(ivar) = nclass
- stypvar(ivar)%cname = 'vtrp_dep'
- stypvar(ivar)%cunits = 'SV'
- stypvar(ivar)%valid_min = 0.
- stypvar(ivar)%valid_max = 10000.
- stypvar(ivar)%clong_name = 'Volume_Transport_per_class'
- stypvar(ivar)%cshort_name = 'vtrp_dep'
-
- IF ( lpm ) THEN
- ivar = ivar + 1 ; iptrpcl = ivar ; imtrpcl = ivar+1
- ipk(ivar) = nclass ; ipk(ivar+1) = nclass
- stypvar(ivar)%cname = 'ptrp_dep' ; stypvar(ivar+1)%cname = 'mtrp_dep'
- stypvar(ivar)%cunits = 'SV' ; stypvar(ivar+1)%cunits = 'SV'
- stypvar(ivar)%valid_min = -500. ; stypvar(ivar+1)%valid_min = -500.
- stypvar(ivar)%valid_max = 500. ; stypvar(ivar+1)%valid_max = 500.
- stypvar(ivar)%clong_name = 'Positive_trp_per_class' ; stypvar(ivar+1)%clong_name = 'Negative_trp_per_class'
- stypvar(ivar)%cshort_name = 'ptrp_dep' ; stypvar(ivar+1)%cshort_name = 'mtrp_dep'
- ivar = ivar + 1
- ENDIF
-
- IF ( lheat ) THEN
- ivar = ivar + 1 ; ihtrpcl = ivar ; istrpcl = ivar+1
- ipk(ivar) = nclass ; ipk(ivar+1) = nclass
- stypvar(ivar)%cname = 'htrp_dep' ; stypvar(ivar+1)%cname = 'strp_dep'
- stypvar(ivar)%cunits = 'PW' ; stypvar(ivar+1)%cunits = 'kt/s'
- stypvar(ivar)%valid_min = -1000. ; stypvar(ivar+1)%valid_min = -1000.
- stypvar(ivar)%valid_max = 1000. ; stypvar(ivar+1)%valid_max = 1000.
- stypvar(ivar)%clong_name = 'Heat_Transport_per_class' ; stypvar(ivar+1)%clong_name = 'Salt_Transport_per_class'
- stypvar(ivar)%cshort_name = 'htrp_dep' ; stypvar(ivar+1)%cshort_name = 'strp_dep'
- ivar = ivar + 1
- ENDIF
- ENDIF
-
- ! Allocate arrays
- ALLOCATE ( zu(npiglo,npjglo), zv(npiglo,npjglo) )
- ALLOCATE ( dwku(npiglo,npjglo), dwkv(npiglo,npjglo) )
- ALLOCATE ( dtrpu(npiglo,npjglo,nclass), dtrpv(npiglo,npjglo,nclass))
- ALLOCATE ( dvoltrpsum(nclass), dvolallegcl(nclass) )
-
- IF ( lpm ) THEN
- ALLOCATE ( dwkup(npiglo,npjglo), dwkvp(npiglo,npjglo) )
- ALLOCATE ( dwkum(npiglo,npjglo), dwkvm(npiglo,npjglo) )
- ALLOCATE ( dtrpup(npiglo,npjglo,nclass), dtrpvp(npiglo,npjglo,nclass))
- ALLOCATE ( dtrpum(npiglo,npjglo,nclass), dtrpvm(npiglo,npjglo,nclass))
- ALLOCATE ( dvoltrpsump(nclass), dvoltrpsumm(nclass) )
- ALLOCATE ( dvolallegclp(nclass), dvolallegclm(nclass) )
- ENDIF
-
- IF ( lheat ) THEN
- ALLOCATE ( zut(npiglo,npjglo), zus(npiglo,npjglo) )
- ALLOCATE ( zvt(npiglo,npjglo), zvs(npiglo,npjglo) )
- ALLOCATE ( dwkut(npiglo,npjglo), dwkus(npiglo,npjglo) )
- ALLOCATE ( dwkvt(npiglo,npjglo), dwkvs(npiglo,npjglo) )
- ALLOCATE ( dtrput(npiglo,npjglo,nclass), dtrpvt(npiglo,npjglo,nclass))
- ALLOCATE ( dtrpus(npiglo,npjglo,nclass), dtrpvs(npiglo,npjglo,nclass))
- ALLOCATE ( dheatrpsum(nclass), dsaltrpsum(nclass) )
- ALLOCATE ( dheatallegcl(nclass), dsaltallegcl(nclass) )
- ENDIF
- !
- ALLOCATE ( e1v(npiglo,npjglo),e3v(npiglo,npjglo) )
- ALLOCATE ( e2u(npiglo,npjglo),e3u(npiglo,npjglo) )
- !
- ALLOCATE ( gphif(npiglo,npjglo) )
- ALLOCATE ( glamf(npiglo,npjglo) )
- ALLOCATE ( gdepw(npk) , tim(npt) )
- !
- ! read metrics and grid position
- e1v(:,:) = getvar(cn_fhgr, cn_ve1v, 1, npiglo, npjglo)
- e2u(:,:) = getvar(cn_fhgr, cn_ve2u, 1, npiglo, npjglo)
-
- glamf(:,:) = getvar(cn_fhgr, cn_glamf, 1,npiglo, npjglo)
- gphif(:,:) = getvar(cn_fhgr, cn_gphif, 1,npiglo, npjglo)
-
- gdepw(:) = getvare3(cn_fzgr, cn_gdepw, npk)
- e31d(:) = getvare3(cn_fzgr, cn_ve3t, npk) ! used only for full step
-
- ! look for nearest level to imeter and setup ilev0 and ilev1 (t-index of class limit)
- ik = 1
- ilev0(1) = 1 ; ilev1(nclass) = npk-1 ! default value if nclass=1
-
- IF ( lobc ) THEN
- ! read u, v on OBC
- IF ( l_zonal ) THEN ! (jpiglo,jpk)
- zuobc(:,:)= getvarxz(cf_ufil, cn_vozocrtx, 1, npiglo, npk)
- zvobc(:,:)= getvarxz(cf_vfil, cn_vomecrty, 1, npiglo, npk)
- ENDIF
- IF ( l_merid ) THEN ! (jpjglo,jpk)
- zuobc(:,:)= getvaryz(cf_ufil, cn_vozocrtx, 1, npjglo, npk)
- zvobc(:,:)= getvaryz(cf_vfil, cn_vomecrty, 1, npjglo, npk)
- ENDIF
- ENDIF
-
- DO jclass = 1, nclass -1
- DO WHILE ( gdepw(ik) < imeter(jclass) )
- ik = ik +1
- END DO
-
- rd1 = ABS(gdepw(ik-1) - imeter(jclass) )
- rd2 = ABS(gdepw(ik ) - imeter(jclass) )
- IF ( rd2 < rd1 ) THEN
- ilev1(jclass ) = ik - 1 ! t-levels index
- ilev0(jclass+1) = ik
- ELSE
- ilev1(jclass ) = ik - 2 ! t-levels index
- ilev0(jclass+1) = ik - 1
- END IF
- END DO
-
- PRINT *, 'Limits : '
- DO jclass = 1, nclass
- PRINT *, ilev0(jclass),ilev1(jclass), gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
- END DO
-
- ! compute the transports at each grid cell
- dtrpu (:,:,:)= 0.d0 ; dtrpv (:,:,:)= 0.d0 ! initialization to 0
-
- IF ( lpm ) THEN
- dtrpup(:,:,:)= 0.d0 ; dtrpvp(:,:,:)= 0.d0
- dtrpum(:,:,:)= 0.d0 ; dtrpvm(:,:,:)= 0.d0
- ENDIF
- IF ( lheat ) THEN
- dtrput(:,:,:)= 0.d0 ; dtrpvt(:,:,:)= 0.d0
- dtrpus(:,:,:)= 0.d0 ; dtrpvs(:,:,:)= 0.d0
- ENDIF
-
- DO jclass = 1, nclass
- DO jk = ilev0(jclass),ilev1(jclass)
- PRINT *,'level ',jk
- ! Get velocities, temperature and salinity fluxes at jk
- IF ( ltest ) THEN
- zu (:,:) = udum ; zv (:,:) = vdum
- IF (lheat) THEN
- zut(:,:) = udum ; zvt(:,:) = vdum
- zus(:,:) = udum ; zvs(:,:) = vdum
- ENDIF
- ELSEIF ( lobc ) THEN
- IF ( l_zonal ) THEN ; zu(:,1)=zuobc(:,jk) ; zv(:,1)=zvobc(:,jk)
- ELSE IF ( l_merid ) THEN ; zu(1,:)=zuobc(:,jk) ; zv(1,:)=zvobc(:,jk) ; ENDIF
- ELSE
- zu (:,:) = getvar(cf_ufil, cn_vozocrtx, jk, npiglo, npjglo, ktime=itime)
- zv (:,:) = getvar(cf_vfil, cn_vomecrty, jk, npiglo, npjglo, ktime=itime)
- IF (lheat) THEN
- zut(:,:) = getvar(cf_tfil, cn_vozout, jk, npiglo, npjglo, ktime=itime)
- zvt(:,:) = getvar(cf_tfil, cn_vomevt, jk, npiglo, npjglo, ktime=itime)
- zus(:,:) = getvar(cf_tfil, cn_vozous, jk, npiglo, npjglo, ktime=itime)
- zvs(:,:) = getvar(cf_tfil, cn_vomevs, jk, npiglo, npjglo, ktime=itime)
- ENDIF
- ENDIF
-
- ! get e3u, e3v at level jk
- IF ( lfull ) THEN
- e3v(:,:) = e31d(jk)
- e3u(:,:) = e31d(jk)
- ELSE
- e3v(:,:) = getvar(cn_fzgr, 'e3v_ps', jk, npiglo, npjglo, ldiom=.TRUE.)
- e3u(:,:) = getvar(cn_fzgr, 'e3u_ps', jk, npiglo, npjglo, ldiom=.TRUE.)
- ENDIF
-
- dwku (:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
- dwkv (:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
-
- IF ( lpm ) THEN
- dwkup = 0.d0 ; dwkum = 0.d0
- WHERE ( zu >= 0. )
- dwkup(:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
- ELSEWHERE
- dwkum(:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
- END WHERE
-
- dwkvp = 0.d0 ; dwkvm = 0.d0
- WHERE ( zv >= 0. )
- dwkvp(:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
- ELSEWHERE
- dwkvm(:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
- END WHERE
- ENDIF
-
- IF ( lheat ) THEN
- dwkut(:,:) = zut(:,:)*e2u(:,:)*e3u(:,:)*1.d0
- dwkvt(:,:) = zvt(:,:)*e1v(:,:)*e3v(:,:)*1.d0
- dwkus(:,:) = zus(:,:)*e2u(:,:)*e3u(:,:)*1.d0
- dwkvs(:,:) = zvs(:,:)*e1v(:,:)*e3v(:,:)*1.d0
- ENDIF
-
- ! integrates vertically
- dtrpu (:,:,jclass) = dtrpu (:,:,jclass) + dwku (:,:)
- dtrpv (:,:,jclass) = dtrpv (:,:,jclass) + dwkv (:,:)
-
- IF ( lpm ) THEN
- dtrpup(:,:,jclass) = dtrpup(:,:,jclass) + dwkup(:,:)
- dtrpvp(:,:,jclass) = dtrpvp(:,:,jclass) + dwkvp(:,:)
- dtrpum(:,:,jclass) = dtrpum(:,:,jclass) + dwkum(:,:)
- dtrpvm(:,:,jclass) = dtrpvm(:,:,jclass) + dwkvm(:,:)
- ENDIF
-
- IF ( lheat ) THEN
- dtrput(:,:,jclass) = dtrput(:,:,jclass) + dwkut(:,:) * rau0*rcp
- dtrpvt(:,:,jclass) = dtrpvt(:,:,jclass) + dwkvt(:,:) * rau0*rcp
- dtrpus(:,:,jclass) = dtrpus(:,:,jclass) + dwkus(:,:)
- dtrpvs(:,:,jclass) = dtrpvs(:,:,jclass) + dwkvs(:,:)
- ENDIF
-
- END DO ! loop to next level
- END DO ! next class
-
- OPEN(numout,FILE=cf_out)
- ! also dump the results on txt files without any comments, some users like it !
- OPEN(numvtrp,FILE=cf_vtrp)
- IF ( lheat ) THEN
- OPEN(numhtrp,FILE=cf_htrp) ; OPEN(numstrp,FILE=cf_strp)
- ENDIF
-
- !################################################################################
- ! enter interactive part
- !################################################################################
- ! initialize all legs arrays and variable to 0
- dvolalleg = 0.d0 ; dvolallegcl(:) = 0.d0
- IF ( lpm ) THEN
- dvolallegp = 0.d0 ; dvolallegclp(:) = 0.d0
- dvolallegm = 0.d0 ; dvolallegclm(:) = 0.d0
- ENDIF
- IF ( lheat ) THEN
- dheatalleg = 0.d0 ; dheatallegcl(:) = 0.d0
- dsaltalleg = 0.d0 ; dsaltallegcl(:) = 0.d0
- ENDIF
- DO
- PRINT *, ' Give name of section (EOF to finish)'
- READ(*,'(a)') csection
- IF (TRIM(csection) == 'EOF' ) THEN
- CLOSE(numout) ; CLOSE(numvtrp)
- IF ( lheat ) THEN
- CLOSE(numhtrp) ; CLOSE(numstrp)
- ENDIF
- EXIT ! infinite DO loop
- ENDIF
-
- ! create output fileset
- cf_outnc = TRIM(csection)//'_transports.nc'
- ncout = create (cf_outnc, 'none', ikx, iky, nclass, cdep='depth_class')
- ierr = createvar (ncout, stypvar, nvarout, ipk, id_varout, cdglobal=TRIM(cglobal) )
- ierr = putheadervar(ncout, cf_ufil, ikx, iky, nclass, pnavlon=rdum, pnavlat=rdum, pdep=rclass )
- tim = getvar1d (cf_ufil, cn_vtimec, npt )
- ierr = putvar1d (ncout, tim, npt, 'T')
-
- PRINT *, ' Give iimin, iimax, ijmin, ijmax '
- READ(*,*) iimin, iimax, ijmin, ijmax
- !! Find the broken line between P1 (iimin,ijmin) and P2 (iimax, ijmax)
- ! ... Initialization
- ii0 = iimin ; ij0 = ijmin ; ii1 = iimax ; ij1 = ijmax
- rxi0 = ii0 ; ryj0 = ij0 ; rxi1 = ii1 ; ryj1 = ij1
-
- ! compute direction of integrations and signs
- !The transport across the section is the dot product of
- !integral(line){(Mx,My)*dS}
- !Mx=integral(u*dz) My=integral(v*dz)) and dS=(dy,-dx)}
-
- !By defining the direction of the integration as
- idirx = SIGN(1,ii1-ii0) !positive to the east or if ii1=ii0
- idiry = SIGN(1,ij1-ij0) !positive to the north or if ij1=ij0
-
- !Then dS=(e2u*idiry,-e1v*idirx)
- !This will produce the following sign convention:
- ! West-to-est line (dx>0, dy=0)=> -My*dx (-ve for a northward flow)
- ! South-to-north (dy>0, dx=0)=> Mx*dy (+ve for an eastward flow)
- norm_u = idiry
- norm_v = -idirx
-
- ! .. Compute equation: ryj = aj rxi + bj [valid in the (i,j) plane]
- IF ( (rxi1 -rxi0) /= 0 ) THEN
- aj = (ryj1 - ryj0 ) / (rxi1 -rxi0)
- bj = ryj0 - aj * rxi0
- ELSE
- aj = 10000. ! flag value
- bj = 0.
- END IF
-
- ! .. Compute equation: rxi = ai ryj + bi [valid in the (i,j) plane]
- IF ( (ryj1 -ryj0) /= 0 ) THEN
- ai = (rxi1 - rxi0 ) / ( ryj1 -ryj0 )
- bi = rxi0 - ai * ryj0
- ELSE
- ai = 10000. ! flag value
- bi = 0.
- END IF
-
- ! .. Compute the integer pathway: a succession of F points
- np=0
- ! .. Chose the strait line with the smallest slope
- IF (ABS(aj) <= 1 ) THEN
- ! ... Here, the best line is y(x)
- ! ... If ii1 < ii0 swap points [ always describe section from left to right ]
- IF (ii1 < ii0 ) THEN
- iitmp = ii0 ; ijtmp = ij0
- ii0 = ii1 ; ij0 = ij1
- ii1 = iitmp ; ij1 = ijtmp
- END IF
-
- ! iist,ijst is the grid offset to pass from F point to either U/V point
- IF ( ij1 >= ij0 ) THEN ! line heading NE
- iist = 1 ; ijst = 1
- ELSE ! line heading SE
- iist = 1 ; ijst = 0
- END IF
-
- ! ... compute the nearest ji point on the line crossing at ji
- DO ji=ii0, ii1
- np=np+1
- IF (np > jpseg) STOP 'np > jpseg !'
- ij=NINT(aj*ji + bj )
- yypt(np) = CMPLX(ji,ij)
- END DO
- ELSE
- ! ... Here, the best line is x(y)
- ! ... If ij1 < ij0 swap points [ always describe section from bottom to top ]
- IF (ij1 < ij0 ) THEN
- iitmp = ii0 ; ijtmp = ij0
- ii0 = ii1 ; ij0 = ij1
- ii1 = iitmp ; ij1 = ijtmp
- END IF
-
- ! iist,ijst is the grid offset to pass from F point to either U/V point
- IF ( ii1 >= ii0 ) THEN
- iist = 1 ; ijst = 1
- ELSE
- iist = 0 ; ijst = 1
- END IF
-
- ! ... compute the nearest ji point on the line crossing at jj
- DO jj=ij0,ij1
- np=np+1
- IF (np > jpseg) STOP 'np > jpseg !'
- ii=NINT(ai*jj + bi)
- yypt(np) = CMPLX(ii,jj)
- END DO
- END IF
-
- !!
- !! Look for intermediate points to be added.
- ! .. The final positions are saved in rxx,ryy
- rxx(1) = REAL(yypt(1))
- ryy(1) = IMAG(yypt(1))
- nn = 1
-
- DO jk=2,np
- ! .. distance between 2 neighbour points
- rd=ABS(yypt(jk)-yypt(jk-1))
- ! .. intermediate points required if rd > 1
- IF ( rd > 1 ) THEN
- CALL interm_pt(yypt, jk, ai, bi, aj, bj, yypti)
- nn=nn+1
- IF (nn > jpseg) STOP 'nn>jpseg !'
- rxx(nn) = REAL(yypti)
- ryy(nn) = IMAG(yypti)
- END IF
- nn=nn+1
- IF (nn > jpseg) STOP 'nn>jpseg !'
- rxx(nn) = REAL(yypt(jk))
- ryy(nn) = IMAG(yypt(jk))
- END DO
- ! record longitude and latitude of initial en endind point of the section
- gla (1) = glamf( INT(rxx(1)), INT(ryy(1)) )
- gphi(1) = gphif( INT(rxx(1)), INT(ryy(1)) )
- gla (2) = glamf( INT(rxx(nn)), INT(ryy(nn)) )
- gphi(2) = gphif( INT(rxx(nn)), INT(ryy(nn)) )
-
- ! Now extract the transport through a section
- ! ... Check whether we need a u velocity or a v velocity
- ! Think that the points are f-points and delimit either a U segment
- ! or a V segment (iist and ijst are set in order to look for the correct
- ! velocity point on the C-grid
- PRINT *, TRIM(csection)
- PRINT *, 'IMIN IMAX JMIN JMAX', iimin, iimax, ijmin, ijmax
- WRITE(numout,*) '% Transport along a section by levels' ,TRIM(csection)
- WRITE(numout,*) '% ---- IMIN IMAX JMIN JMAX'
-
- dvoltrpbrtp = 0.d0
- dvoltrpbrtpp = 0.d0
- dvoltrpbrtpm = 0.d0
- dheatrpbrtp = 0.d0
- dsaltrpbrtp = 0.d0
- DO jclass=1,nclass
- dvoltrpsum(jclass) = 0.d0
- IF ( lpm ) THEN
- dvoltrpsump(jclass) = 0.d0
- dvoltrpsumm(jclass) = 0.d0
- ENDIF
- IF ( lheat ) THEN
- dheatrpsum(jclass) = 0.d0
- dsaltrpsum(jclass) = 0.d0
- ENDIF
-
- ! segment jseg is a line between (rxx(jseg),ryy(jseg)) and (rxx(jseg+1),ryy(jseg+1))
- DO jseg = 1, nn-1
- ii0=rxx(jseg)
- ij0=ryy(jseg)
- IF ( rxx(jseg) == rxx(jseg+1) ) THEN ! meridional segment, use U velocity
- dvoltrp(jseg)= dtrpu (ii0,ij0+ijst,jclass)*norm_u
-
- IF ( lpm ) THEN
- dvoltrpp(jseg)= dtrpup(ii0,ij0+ijst,jclass)*norm_u
- dvoltrpm(jseg)= dtrpum(ii0,ij0+ijst,jclass)*norm_u
- ENDIF
-
- IF ( lheat ) THEN
- dheatrp(jseg)= dtrput(ii0,ij0+ijst,jclass)*norm_u
- dsaltrp(jseg)= dtrpus(ii0,ij0+ijst,jclass)*norm_u
- ENDIF
- ELSE IF ( ryy(jseg) == ryy(jseg+1) ) THEN ! zonal segment, use V velocity
- dvoltrp(jseg)=dtrpv (ii0+iist,ij0,jclass)*norm_v
-
- IF ( lpm ) THEN
- dvoltrpp(jseg)=dtrpvp(ii0+iist,ij0,jclass)*norm_v
- dvoltrpm(jseg)=dtrpvm(ii0+iist,ij0,jclass)*norm_v
- ENDIF
-
- IF ( lheat ) THEN
- dheatrp(jseg)=dtrpvt(ii0+iist,ij0,jclass)*norm_v
- dsaltrp(jseg)=dtrpvs(ii0+iist,ij0,jclass)*norm_v
- ENDIF
- ELSE
- PRINT *,' ERROR :', rxx(jseg),ryy(jseg),rxx(jseg+1),ryy(jseg+1) ! likely to never happen !
- END IF
- dvoltrpsum(jclass) = dvoltrpsum(jclass) + dvoltrp(jseg)
- IF ( lpm ) THEN
- dvoltrpsump(jclass) = dvoltrpsump(jclass) + dvoltrpp(jseg)
- dvoltrpsumm(jclass) = dvoltrpsumm(jclass) + dvoltrpm(jseg)
- ENDIF
- IF ( lheat ) THEN
- dheatrpsum(jclass) = dheatrpsum(jclass) + dheatrp(jseg)
- dsaltrpsum(jclass) = dsaltrpsum(jclass) + dsaltrp(jseg)
- ENDIF
- END DO ! next segment
-
- ! Ascii outputs :
- IF (jclass == 1 ) THEN ! print header when it is the first class
- PRINT '(a,2f8.2,a,2f8.2)', 'FROM (LON LAT): ', gla(1),gphi(1),' TO (LON LAT) ', gla(2), gphi(2)
- WRITE(numout,*) '% ---- LONmin LATmin LONmax LATmax'
- WRITE(numout,*) '% Top(m) Bottom(m) MassTrans(Sv) HeatTrans(PW) SaltTrans(kt/s)'
- WRITE(numout,*) 0 ,iimin, iimax, ijmin, ijmax
- WRITE(numout,9003) 0. ,gla(1),gphi(1), gla(2), gphi(2)
- ENDIF
-
- PRINT *, gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
- PRINT *, ' Mass transport : ', dvoltrpsum(jclass)/1.e6,' SV'
- WRITE(numvtrp,'(e12.6)') dvoltrpsum(jclass)
- IF ( lpm ) THEN
- PRINT *, ' Positive Mass transport : ', dvoltrpsump(jclass)/1.e6,' SV'
- PRINT *, ' Negative Mass transport : ', dvoltrpsumm(jclass)/1.e6,' SV'
- WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), &
- & dvoltrpsum(jclass)/1.e6, dvoltrpsump(jclass)/1.e6, dvoltrpsumm(jclass)/1.e6
- WRITE(numvtrp,'(e12.6)') dvoltrpsump(jclass)
- WRITE(numvtrp,'(e12.6)') dvoltrpsumm(jclass)
- ENDIF
-
- IF ( lheat ) THEN
- PRINT *, ' Heat transport : ', dheatrpsum(jclass)/1.e15,' PW'
- PRINT *, ' Salt transport : ', dsaltrpsum(jclass)/1.e6,' kT/s'
- WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), &
- & dvoltrpsum(jclass)/1.e6, dheatrpsum(jclass)/1.e15, dsaltrpsum(jclass)/1.e6
- WRITE(numhtrp,'(e12.6)') dheatrpsum(jclass)
- WRITE(numstrp,'(e12.6)') dsaltrpsum(jclass)
- ELSE
- IF ( .NOT. lpm ) WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), dvoltrpsum(jclass)/1.e6
- ENDIF
-
- ! netcdf output
- IF ( nclass > 1 ) THEN
- rdum(1,1) = REAL(dvoltrpsum(jclass)/1.e6)
- ierr = putvar(ncout,id_varout(ivtrpcl), rdum, jclass, 1, 1, ktime=itime )
- IF ( lpm ) THEN
- rdum(1,1) = REAL(dvoltrpsump(jclass)/1.e6)
- ierr = putvar(ncout,id_varout(iptrpcl), rdum, jclass, 1, 1, ktime=itime )
- rdum(1,1) = REAL(dvoltrpsumm(jclass)/1.e6)
- ierr = putvar(ncout,id_varout(imtrpcl), rdum, jclass, 1, 1, ktime=itime )
- ENDIF
- IF ( lheat ) THEN
- rdum(1,1) = REAL(dheatrpsum(jclass)/1.e15)
- ierr = putvar(ncout,id_varout(ihtrpcl), rdum, jclass, 1, 1, ktime=itime )
- rdum(1,1) = REAL(dsaltrpsum(jclass)/1.e6)
- ierr = putvar(ncout,id_varout(istrpcl), rdum, jclass, 1, 1, ktime=itime )
- ENDIF
- ENDIF
- rdum(1,1) = REAL(gdepw(ilev0(jclass)))
- ierr = putvar(ncout,id_varout(itop), rdum, jclass, 1, 1, ktime=itime )
- rdum(1,1) = REAL(gdepw(ilev1(jclass)+1))
- ierr = putvar(ncout,id_varout(ibot), rdum, jclass, 1, 1, ktime=itime )
-
- dvoltrpbrtp = dvoltrpbrtp + dvoltrpsum(jclass)
- IF ( lpm ) THEN
- dvoltrpbrtpp = dvoltrpbrtpp + dvoltrpsump(jclass)
- dvoltrpbrtpm = dvoltrpbrtpm + dvoltrpsumm(jclass)
- ENDIF
- IF ( lheat) THEN
- dheatrpbrtp = dheatrpbrtp + dheatrpsum(jclass)
- dsaltrpbrtp = dsaltrpbrtp + dsaltrpsum(jclass)
- ENDIF
- ! save sum over legs
- dvolallegcl(jclass) = dvolallegcl(jclass) + dvoltrpsum(jclass)
- IF ( lpm ) THEN
- dvolallegclp(jclass) = dvolallegclp(jclass) + dvoltrpsump(jclass)
- dvolallegclm(jclass) = dvolallegclm(jclass) + dvoltrpsumm(jclass)
- ENDIF
- IF ( lheat ) THEN
- dheatallegcl(jclass) = dheatallegcl(jclass) + dheatrpsum(jclass)
- dsaltallegcl(jclass) = dsaltallegcl(jclass) + dsaltrpsum(jclass)
- ENDIF
- END DO ! next class
- ! save sum over legs
- dvolalleg = dvolalleg + dvoltrpbrtp
- IF ( lpm ) THEN
- dvolallegp = dvolallegp + dvoltrpbrtpp
- dvolallegm = dvolallegm + dvoltrpbrtpm
- ENDIF
- IF ( lheat ) THEN
- dheatalleg = dheatalleg + dheatrpbrtp
- dsaltalleg = dsaltalleg + dsaltrpbrtp
- ENDIF
-
- IF ( nclass > 1 ) THEN
- PRINT *, ' ====================================================='
- PRINT *, ' total Mass transport : ', dvoltrpbrtp/1.e6,' SV'
- IF ( lpm ) THEN
- PRINT *, ' total positive transport : ', dvoltrpbrtpp/1.e6,' SV'
- PRINT *, ' total negative transport : ', dvoltrpbrtpm/1.e6,' SV'
- ENDIF
- IF ( lheat ) THEN
- PRINT *, ' total Heat transport : ', dheatrpbrtp/1.e15,' PW'
- PRINT *, ' total Salt transport : ', dsaltrpbrtp/1.e6,' kT/s'
- ENDIF
- ENDIF
- ierr = putvar0d(ncout,id_varout(ivtrp), REAL(dvoltrpbrtp/1.e6) )
- IF ( lpm ) THEN
- ierr = putvar0d(ncout,id_varout(iptrp), REAL(dvoltrpbrtpp/1.e6) )
- ierr = putvar0d(ncout,id_varout(imtrp), REAL(dvoltrpbrtpm/1.e6) )
- ENDIF
- IF ( lheat ) THEN
- ierr = putvar0d(ncout,id_varout(ihtrp), REAL(dheatrpbrtp/1.e15) )
- ierr = putvar0d(ncout,id_varout(istrp), REAL(dsaltrpbrtp/1.e6 ) )
- ENDIF
- ierr = putvar0d(ncout,id_varout(ilonmin), REAL(gla(1)) )
- ierr = putvar0d(ncout,id_varout(ilonmax), REAL(gla(2)) )
- ierr = putvar0d(ncout,id_varout(ilatmin), REAL(gphi(1)) )
- ierr = putvar0d(ncout,id_varout(ilatmax), REAL(gphi(2)) )
- ierr = closeout(ncout)
- END DO ! infinite loop : gets out when input is EOF
-
-
- PRINT *,' '
- PRINT *,' Overall transports (sum of all legs done so far)'
- DO jclass = 1, nclass
- PRINT *, gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
- PRINT *, ' Mass transport : ', dvolallegcl(jclass)/1.e6,' SV'
- IF ( lpm ) THEN
- PRINT *, ' Positive Mass transport : ', dvolallegclp(jclass)/1.e6,' SV'
- PRINT *, ' Negative Mass transport : ', dvolallegclm(jclass)/1.e6,' SV'
- ENDIF
-
- IF ( lheat ) THEN
- PRINT *, ' Heat transport : ', dheatallegcl(jclass)/1.e15,' PW'
- PRINT *, ' Salt transport : ', dsaltallegcl(jclass)/1.e6,' kT/s'
- ENDIF
- ENDDO
-
- IF ( nclass > 1 ) THEN
- PRINT *, ' ====================================================='
- PRINT *, ' Mass transport : ', dvolalleg/1.e6,' SV'
- IF ( lpm ) THEN
- PRINT *, ' positive transport : ', dvolallegp/1.e6,' SV'
- PRINT *, ' negative transport : ', dvolallegm/1.e6,' SV'
- ENDIF
- IF ( lheat ) THEN
- PRINT *, ' heat transport : ', dheatalleg/1.e15,' PW'
- PRINT *, ' salt transport : ', dsaltalleg/1.e6,' kT/s'
- ENDIF
- ENDIF
-
-
-9000 FORMAT(I4,6(f9.3,f8.4))
-9001 FORMAT(I4,6(f9.2,f9.3))
-9002 FORMAT(f9.0,f9.0,f9.2,f9.2,f9.2)
-9003 FORMAT(f9.2,f9.2,f9.2,f9.2,f9.2)
-
-CONTAINS
- SUBROUTINE interm_pt (ydpt, kk, pai, pbi, paj, pbj, ydpti)
- !!---------------------------------------------------------------------
- !! *** ROUTINE nterm_pt ***
- !!
- !! ** Purpose : Find the best intermediate points on a pathway.
- !!
- !! ** Method : ydpt : complex vector of the positions of the nearest points
- !! kk : current working index
- !! pai, pbi : slope and original ordinate of x(y)
- !! paj, pbj : slope and original ordinate of y(x)
- !! ydpti : Complex holding the position of intermediate point
- !!
- !! ** Reference : 19/07/1999 : J.M. Molines in Clipper
- !!----------------------------------------------------------------------
- COMPLEX, DIMENSION(:), INTENT(in ) :: ydpt
- COMPLEX, INTENT(out) :: ydpti
- REAL(KIND=4), INTENT(in ) :: pai, pbi, paj, pbj
- INTEGER(KIND=4), INTENT(in ) :: kk
- ! ... local
- COMPLEX :: ylptmp1, ylptmp2
- REAL(KIND=4) :: za0, zb0
- REAL(KIND=4) :: za1, zb1
- REAL(KIND=4) :: zd1, zd2
- REAL(KIND=4) :: zxm, zym
- REAL(KIND=4) :: zxp, zyp
- !!----------------------------------------------------------------------
- ! ... Determines whether we use y(x) or x(y):
- IF (ABS(paj) <= 1) THEN
- ! ..... use y(x)
- ! ... possible intermediate points:
- ylptmp1=ydpt(kk-1)+(1.,0.) ! M1
- ylptmp2=ydpt(kk-1)+CMPLX(0.,SIGN(1.,paj)) ! M2
- !
- ! ... M1 is the candidate point:
- zxm=REAL(ylptmp1)
- zym=IMAG(ylptmp1)
- za0=paj
- zb0=pbj
- !
- za1=-1./za0
- zb1=zym - za1*zxm
- ! ... P1 is the projection of M1 on the strait line
- zxp=-(zb1-zb0)/(za1-za0)
- zyp=za0*zxp + zb0
- ! ... zd1 is the distance M1P1
- zd1=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
- !
- ! ... M2 is the candidate point:
- zxm=REAL(ylptmp2)
- zym=IMAG(ylptmp2)
- za1=-1./za0
- zb1=zym - za1*zxm
- ! ... P2 is the projection of M2 on the strait line
- zxp=-(zb1-zb0)/(za1-za0)
- zyp=za0*zxp + zb0
- ! ... zd2 is the distance M2P2
- zd2=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
- ! ... chose the smallest (zd1,zd2)
- IF (zd2 <= zd1) THEN
- ydpti=ylptmp2 ! use M2
- ELSE
- ydpti=ylptmp1 ! use M1
- END IF
- !
- ELSE
- ! ... use x(y)
- ! ... possible intermediate points:
- ylptmp1=ydpt(kk-1)+CMPLX(SIGN(1.,pai),0.) ! M1
- ylptmp2=ydpt(kk-1)+(0.,1.) ! M2
- !
- ! ... M1 is the candidate point:
- zxm=REAL(ylptmp1)
- zym=IMAG(ylptmp1)
- za0=pai
- zb0=pbi
- !
- za1=-1./za0
- zb1=zxm - za1*zym
- zyp=-(zb1-zb0)/(za1-za0)
- zxp=za0*zyp + zb0
- zd1=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
- !
- zxm=REAL(ylptmp2)
- zym=IMAG(ylptmp2)
- za1=-1./za0
- zb1=zxm - za1*zym
- zyp=-(zb1-zb0)/(za1-za0)
- zxp=za0*zyp + zb0
- zd2=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
- IF (zd2 <= zd1) THEN
- ydpti=ylptmp2
- ELSE
- ydpti=ylptmp1
- END IF
- END IF
- END SUBROUTINE interm_pt
-
-END PROGRAM cdftransport
+ !!======================================================================
+ !! *** PROGRAM cdftransport ***
+ !!=====================================================================
+ !! ** Purpose : Compute Transports across a section.
+ !! By default, mass (Sv) and heat(PW)/salt(kT/s) transports
+ !! are computed unless -noheat option is used (mass
+ !! transport only).
+ !!
+ !! ** Method : The begining and end point of the section are given in
+ !! term of F-points index. A broken line joining successive
+ !! F-points is defined between the begining and end point
+ !! of the section. Therefore each segment between F-points
+ !! is either a zonal or meridional segment corresponding to
+ !! V or U velocity component. Doing so, the volume conservation
+ !! is ensured as velocities are not interpolated, and stay
+ !! on the native model grid.
+ !! The section name and the begin/end point of a section are
+ !! read from standard input, till 'EOF' is given as section
+ !! name. This make possible to give a bunch of sections in
+ !! an ASCII files and use the < redirection.
+ !! SIGN CONVENTION : The transport is positive when the flow cross
+ !! the section to the right, negative otherwise. This depends
+ !! on the sense the section is described. With this convention
+ !! The algebric sum of transports accross sections forming a
+ !! closed area is 0.
+ !! OPTIONS :
+ !! -full : full step case
+ !! -noheat : only mass transport is computed.
+ !! -time : specify the time frame to be used
+ !! -zlimit : transports can be computed in different depth layers
+ !! defined by their depth limit
+ !! REQUIREMENT :
+ !! mesh-mask file are required in the current directory.
+ !!
+ !!
+ !! History : 2.1 : 01/2005 : J.M. Molines : Original code
+ !! 2.1 : 07/2009 : R. Dussin : add cdf output
+ !! 2.1 : 01/2010 : M.A. Balmaseda : Change integration signs
+ !! so that the transport across a segment is
+ !! independent of the chosen trajectory.
+ !! 3.0 : 04/2011 : J.M. Molines : Doctor norm + Lic.
+ !!----------------------------------------------------------------------
+ !!----------------------------------------------------------------------
+ !! routines : description
+ !! interm_pt : choose intermediate points on a broken line.
+ !!----------------------------------------------------------------------
+ USE cdfio
+ USE modcdfnames
+ USE modutils ! for global attribute
+ !!----------------------------------------------------------------------
+ !! CDFTOOLS_3.0 , MEOM 2011
+ !! $Id$
+ !! Copyright (c) 2011, J.-M. Molines
+ !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
+ !!----------------------------------------------------------------------
+ IMPLICIT NONE
+
+ INTEGER(KIND=4) :: jclass, jseg ! dummy loop index
+ INTEGER(KIND=4) :: ji, jj, jk ! dummy loop index
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: imeter ! limit beetween depth level, in m (nclass -1)
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ilev0, ilev1 ! limit in levels (nclass)
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_varout ! Netcdf output
+ INTEGER(KIND=4) :: ipos ! working integer (position of ' ' in strings)
+ INTEGER(KIND=4) :: ncout, ierr ! for netcdf output
+ INTEGER(KIND=4) :: nvarout=12 ! number of values to write in cdf output
+ INTEGER(KIND=4) :: ivtrp ! var index of volume transport (barotrope)
+ INTEGER(KIND=4) :: iptrp ! var index of volume transport (barotrope)
+ INTEGER(KIND=4) :: imtrp ! var index of volume transport (barotrope)
+ INTEGER(KIND=4) :: ihtrp ! var index of heat transport (barotrope)
+ INTEGER(KIND=4) :: istrp ! var index of sal transport (barotrope)
+ INTEGER(KIND=4) :: ivtrpcl ! var index of volume transport (p. class)
+ INTEGER(KIND=4) :: iptrpcl ! var index of volume transport (p. class)
+ INTEGER(KIND=4) :: imtrpcl ! var index of volume transport (p. class)
+ INTEGER(KIND=4) :: ihtrpcl ! var index of heat transport (p. class)
+ INTEGER(KIND=4) :: istrpcl ! var index of sal transport (p. class)
+ INTEGER(KIND=4) :: ilonmin ! var index of starting section longitude
+ INTEGER(KIND=4) :: ilonmax ! var index of ending section longitude
+ INTEGER(KIND=4) :: ilatmin ! var index of starting section latitude
+ INTEGER(KIND=4) :: ilatmax ! var index of ending section latitude
+ INTEGER(KIND=4) :: itop ! var index of top depth class
+ INTEGER(KIND=4) :: ibot ! var index of bottom depth class
+ INTEGER(KIND=4) :: ikx=1, iky=1 ! dims of netcdf output file
+ INTEGER(KIND=4) :: numout = 10 ! logical unit for output file (overall)
+ INTEGER(KIND=4) :: numvtrp = 11 ! logical unit for volume transport file
+ INTEGER(KIND=4) :: numhtrp = 12 ! logical unit for heat transport file
+ INTEGER(KIND=4) :: numstrp = 14 ! logical unit for salt trp file
+ INTEGER(KIND=4) :: nclass ! number of depth class
+ INTEGER(KIND=4) :: narg, iargc ! command line
+ INTEGER(KIND=4) :: ijarg, nxtarg ! " "
+ INTEGER(KIND=4) :: npiglo, npjglo ! size of the domain
+ INTEGER(KIND=4) :: npk, npt ! size of the domain
+ INTEGER(KIND=4) :: iimin, iimax ! i-limit of the section
+ INTEGER(KIND=4) :: ijmin, ijmax ! j-limit of the section
+ INTEGER(KIND=4) :: ivar, itime ! working integer
+ INTEGER(KIND=4) :: ii, ij, ik ! working integer
+ INTEGER(KIND=4), PARAMETER :: jpseg=10000 ! used for broken line algorithm
+ INTEGER(KIND=4) :: ii0, ij0 ! " " "
+ INTEGER(KIND=4) :: ii1, ij1 ! " " "
+ INTEGER(KIND=4) :: iitmp, ijtmp ! " " "
+ INTEGER(KIND=4) :: np, nn ! segment counters,
+ INTEGER(KIND=4) :: iist, ijst ! local point offset for velocity
+ INTEGER(KIND=4) :: norm_u, norm_v ! normalization factor (sign of normal transport)
+ INTEGER(KIND=4) :: idirx, idiry ! sense of description of the section
+
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e1v, e2u ! horizontal metric
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e3u, e3v ! vertical metric
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamf ! longitudes of F points
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: gphif ! latitudes of F points
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zu, zut, zus ! Zonal velocities and uT uS
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zv, zvt, zvs ! Meridional velocities and uT uS
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rdum ! dummy (1x1) array for ncdf output
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zuobc, zvobc ! arrays for OBC files (vertical slice)
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tim ! time counter
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdepw ! depth at layer interface
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: e31d ! vertical metric in case of full step
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: rclass ! vertical metric in case of full step
+ REAL(KIND=4), DIMENSION(2) :: gla, gphi ! lon/lat of the begining/end of section (f point)
+ REAL(KIND=4), DIMENSION(jpseg) :: rxx, ryy ! working variables
+ 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)
+
+ ! at every model point
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwku, dwkv ! volume transport at each cell boundary
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkut, dwkvt ! heat transport at each cell boundary
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkus, dwkvs ! salt transport at each cell boundary
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkup, dwkvp ! volume transport in the positive direction
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: dwkum, dwkvm ! volume transport in the negatibe direction
+ REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpu, dtrpv ! volume transport integrated in depth class
+ REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrput, dtrpvt ! heat transport integrated in depth class
+ REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpus, dtrpvs ! salt transport integrated in depth class
+ REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpup, dtrpvp ! volume transport integrated in depth class (positive)
+ REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: dtrpum, dtrpvm ! volume transport integrated in depth class (negative)
+ ! for a given section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsum ! volume transport by depth class across section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsump ! volume transport by depth class across section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvoltrpsumm ! volume transport by depth class across section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dheatrpsum ! heat transport by depth class across section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsaltrpsum ! salt transport by depth class across section
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegcl ! over all leg volume transport by depth class
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegclp ! over all leg volume transport by depth class +
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dvolallegclm ! over all leg volume transport by depth class -
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dheatallegcl ! over all leg heat transport by depth class
+ REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: dsaltallegcl ! over all leg salt transport by depth class
+ REAL(KIND=8), DIMENSION(jpseg) :: dvoltrp ! volume transport across each segment of a section
+ REAL(KIND=8), DIMENSION(jpseg) :: dvoltrpp ! volume transport across each segment of a section
+ REAL(KIND=8), DIMENSION(jpseg) :: dvoltrpm ! volume transport across each segment of a section
+ REAL(KIND=8), DIMENSION(jpseg) :: dheatrp ! heat transport across each segment of a section
+ REAL(KIND=8), DIMENSION(jpseg) :: dsaltrp ! salt transport across each segment of a section
+ REAL(KIND=8) :: dvoltrpbrtp ! volume transport integrated over the whole depth
+ REAL(KIND=8) :: dvoltrpbrtpp ! volume transport integrated over the whole depth
+ REAL(KIND=8) :: dvoltrpbrtpm ! volume transport integrated over the whole depth
+ REAL(KIND=8) :: dheatrpbrtp ! heat transport integrated over the whole depth
+ REAL(KIND=8) :: dsaltrpbrtp ! salt transport integrated over the whole depth
+ REAL(KIND=8) :: dvolalleg ! over all leg sum of volume transport
+ REAL(KIND=8) :: dvolallegp ! over all leg sum of volume transport +
+ REAL(KIND=8) :: dvolallegm ! over all leg sum of volume transport -
+ REAL(KIND=8) :: dheatalleg ! over all leg sum of heat transport
+ REAL(KIND=8) :: dsaltalleg ! over all leg sum of salt transport
+
+ COMPLEX, DIMENSION(jpseg) :: yypt ! array of points coordinates in a section
+ COMPLEX :: yypti ! working point
+
+ TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvar ! structure of output
+
+ CHARACTER(LEN=256) :: cf_tfil ! VT file (in)
+ CHARACTER(LEN=256) :: cf_ufil ! U file (in)
+ CHARACTER(LEN=256) :: cf_vfil ! V file (in)
+ CHARACTER(LEN=256) :: cf_out='section_trp.dat' ! output file name (ASCII)
+ CHARACTER(LEN=256) :: cf_outnc ! output netcdf file
+ CHARACTER(LEN=256) :: cf_vtrp='vtrp.txt' ! output volume transport file
+ CHARACTER(LEN=256) :: cf_htrp='htrp.txt' ! output heat transport file
+ CHARACTER(LEN=256) :: cf_strp='strp.txt' ! output salt transport file
+ CHARACTER(LEN=256) :: csection ! section names
+ CHARACTER(LEN=256) :: cvarname ! variable names (root)
+ CHARACTER(LEN=256) :: clongname ! variable longname (root)
+ CHARACTER(LEN=512) :: cglobal ! global attribute
+ CHARACTER(LEN=256) :: cldum ! dummy char variable
+ CHARACTER(LEN=256) :: cline ! dummy char variable
+ CHARACTER(LEN=256), DIMENSION(3) :: cldumt ! dummy char variable
+
+ LOGICAL :: ltest = .FALSE. ! flag for test case
+ LOGICAL :: lfull = .FALSE. ! flag for full step case
+ 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
+ LOGICAL :: lobc = .FALSE. ! flag for obc input files
+ LOGICAL :: l_merid = .FALSE. ! flag for meridional obc
+ LOGICAL :: l_zonal = .FALSE. ! flag for zonal obc
+ !!----------------------------------------------------------------------
+ CALL ReadCdfNames()
+
+ narg= iargc()
+ ! Print usage if no argument
+ IF ( narg == 0 ) THEN
+ PRINT *,' usage : cdftransport [-test u v ] [-noheat ] [-plus_minus ] [-obc]...'
+ PRINT *,' ... [VT-file] U-file V-file [-full] |-time jt] ...'
+ PRINT *,' ... [-time jt ] [-zlimit limits of level]'
+ PRINT *,' '
+ PRINT *,' PURPOSE :'
+ PRINT *,' Compute the transports accross a section.'
+ PRINT *,' The name of the section and the imin, imax, jmin, jmax for the section '
+ PRINT *,' is read from the standard input. To finish the program use the key name'
+ PRINT *,' ''EOF'' for the section name.'
+ PRINT *,' OBC U,V files can be used if -obc option is specified.'
+ PRINT *,' '
+ PRINT *,' ARGUMENTS :'
+ PRINT *,' [VT-file ] : netcdf file with mean values of vt, vs, ut, us for heat and'
+ PRINT *,' salt transport. If options -noheat or -plus_minus are used'
+ PRINT *,' this file name must be omitted.'
+ PRINT *,' [U-file ] : netcdf file with the zonal velocity component.'
+ PRINT *,' [V-file ] : netcdf file with the meridional velocity component.'
+ PRINT *,' '
+ PRINT *,' OPTIONS :'
+ PRINT *,' [-test u v ]: use constant the u and v velocity components for sign '
+ PRINT *,' test purpose.'
+ PRINT *,' [-noheat ] : use when heat and salt transport are not requested.'
+ PRINT *,' This option must come before the file names, and if used'
+ PRINT *,' VT file must not be given.'
+ PRINT *,' [ -plus_minus or -pm ] : separate positive and negative contribution to'
+ PRINT *,' the volume transport. This option implicitly set -noheat,'
+ PRINT *,' and must be used before the file names.'
+ PRINT *,' [-obc ] : indicates that input files are obc files (vertical slices)'
+ PRINT *,' Take care that for this case, mesh files must be adapted.'
+ PRINT *,' This option implicitly set -noheat, and must be used before'
+ PRINT *,' the file names.'
+ PRINT *,' [-full ] : use for full step configurations.'
+ PRINT *,' [-time jt ]: compute transports for time index jt. Default is 1.'
+ PRINT *,' [-zlimit list of depth] : Specify depths limits defining layers where the'
+ PRINT *,' transports will be computed. If not used, the transports '
+ PRINT *,' are computed for the whole water column. If used, this '
+ PRINT *,' option must be the last on the command line.'
+ PRINT *,' '
+ PRINT *,' REQUIRED FILES :'
+ PRINT *,' Files ',TRIM(cn_fhgr),', ',TRIM(cn_fzgr),' must be in the current directory.'
+ PRINT *,' '
+ PRINT *,' OUTPUT : '
+ PRINT *,' - Standard output '
+ PRINT *,' - ASCII file reflecting the standard output: section_trp.dat'
+ PRINT *,' - ASCII files for volume, heat and salt transport: vtrp.txt, htrp.txt '
+ PRINT *,' and strp.txt.'
+ PRINT *,' - Netcdf files for each section. name of the file is buildt'
+ PRINT *,' from section name.'
+ PRINT *,' '
+ PRINT *,' SEE ALSO :'
+ PRINT *,' cdfsigtrp'
+ PRINT *,' '
+ STOP
+ ENDIF
+
+ itime = 1
+ nclass = 1
+ ijarg = 1
+ CALL SetGlobalAtt(cglobal)
+
+ ! Browse command line for arguments and/or options
+ DO WHILE (ijarg <= narg )
+ CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1
+ SELECT CASE ( cldum )
+ CASE ('-test ')
+ CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) udum
+ CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) vdum
+ ltest = .TRUE.
+
+ CASE ('-full' )
+ lfull = .TRUE.
+
+ CASE ('-noheat' ) ! it must be called before the list of files
+ lheat = .FALSE.
+
+ CASE ('-time' )
+ CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) itime
+
+ CASE ('-plus_minus', '-pm' )
+ lpm = .TRUE.
+ lheat = .FALSE.
+
+ CASE ('-obc' )
+ lobc = .TRUE.
+ lheat = .FALSE.
+
+ CASE ('-zlimit' ) ! this should be the last option on the line
+ nxtarg = ijarg - 1
+ nclass = narg - nxtarg + 1
+ ALLOCATE ( imeter(nclass -1) ) ! if no zlimit option, this array is never used
+ DO jclass =1, nclass -1
+ CALL getarg (ijarg, cldum) ; ijarg = ijarg + 1 ; READ(cldum,*) imeter(jclass)
+ END DO
+
+ CASE DEFAULT
+ ijarg = ijarg -1 ! re-read argument in this case
+ IF ( lheat) THEN
+ CALL getarg (ijarg, cf_tfil) ; ijarg = ijarg + 1
+ ENDIF
+ CALL getarg (ijarg, cf_ufil) ; ijarg = ijarg + 1
+ CALL getarg (ijarg, cf_vfil) ; ijarg = ijarg + 1
+ END SELECT
+ END DO
+
+ ! checking if all required files are available
+ lchk = lchk .OR. chkfile(cn_fzgr)
+ lchk = lchk .OR. chkfile(cn_fhgr)
+ IF ( ltest ) THEN
+ ! OK
+ ELSE
+ lchk = lchk .OR. chkfile(cf_ufil)
+ lchk = lchk .OR. chkfile(cf_vfil)
+ IF (lheat) THEN
+ lchk = lchk .OR. chkfile(cf_tfil)
+ ENDIF
+ ENDIF
+ IF ( lchk ) STOP ! missing files
+
+ ! adjust the number of output variables according to options
+ IF ( nclass > 1 ) THEN
+ IF ( lheat ) THEN
+ nvarout = 12
+ ELSE
+ nvarout = 8
+ ENDIF
+ IF ( lpm ) nvarout=nvarout+4
+ ELSE
+ IF ( lheat ) THEN
+ nvarout = 9
+ ELSE
+ nvarout = 7
+ ENDIF
+ IF ( lpm ) nvarout=nvarout+2
+ ENDIF
+
+ ALLOCATE ( ilev0(nclass), ilev1(nclass), rclass(nclass) )
+ rclass=(/(jclass, jclass=1,nclass)/)
+
+ npiglo = getdim (cf_ufil,cn_x)
+ npjglo = getdim (cf_ufil,cn_y)
+ npk = getdim (cf_ufil,cn_z)
+ npt = getdim (cf_ufil,cn_t)
+
+ PRINT *, 'npiglo =', npiglo
+ PRINT *, 'npjglo =', npjglo
+ PRINT *, 'npk =', npk
+ PRINT *, 'npt =', npt
+
+ IF ( lobc ) THEN ! if lobc false, l_merid and l_zonal are false (default)
+ IF ( npiglo == 1 ) THEN
+ l_merid=.TRUE.
+ ALLOCATE (zuobc(npjglo,npk), zvobc(npjglo,npk) )
+ PRINT *,' Meridional OBC'
+ ENDIF
+
+ IF ( npjglo == 1 ) THEN
+ l_zonal=.TRUE.
+ ALLOCATE (zuobc(npiglo,npk), zvobc(npiglo,npk) )
+ PRINT *,' Zonal OBC'
+ ENDIF
+ ENDIF
+
+ ALLOCATE ( e31d(npk) )
+
+ ! define new variables for output
+ ALLOCATE ( stypvar(nvarout), ipk(nvarout), id_varout(nvarout) )
+ ALLOCATE ( rdum(1,1) )
+
+ rdum(:,:)=0.e0
+
+ ! Allocate arrays
+ ALLOCATE ( zu(npiglo,npjglo), zv(npiglo,npjglo) )
+ ALLOCATE ( dwku(npiglo,npjglo), dwkv(npiglo,npjglo) )
+ ALLOCATE ( dtrpu(npiglo,npjglo,nclass), dtrpv(npiglo,npjglo,nclass))
+ ALLOCATE ( dvoltrpsum(nclass), dvolallegcl(nclass) )
+
+ IF ( lpm ) THEN
+ ALLOCATE ( dwkup(npiglo,npjglo), dwkvp(npiglo,npjglo) )
+ ALLOCATE ( dwkum(npiglo,npjglo), dwkvm(npiglo,npjglo) )
+ ALLOCATE ( dtrpup(npiglo,npjglo,nclass), dtrpvp(npiglo,npjglo,nclass))
+ ALLOCATE ( dtrpum(npiglo,npjglo,nclass), dtrpvm(npiglo,npjglo,nclass))
+ ALLOCATE ( dvoltrpsump(nclass), dvoltrpsumm(nclass) )
+ ALLOCATE ( dvolallegclp(nclass), dvolallegclm(nclass) )
+ ENDIF
+
+ IF ( lheat ) THEN
+ ALLOCATE ( zut(npiglo,npjglo), zus(npiglo,npjglo) )
+ ALLOCATE ( zvt(npiglo,npjglo), zvs(npiglo,npjglo) )
+ ALLOCATE ( dwkut(npiglo,npjglo), dwkus(npiglo,npjglo) )
+ ALLOCATE ( dwkvt(npiglo,npjglo), dwkvs(npiglo,npjglo) )
+ ALLOCATE ( dtrput(npiglo,npjglo,nclass), dtrpvt(npiglo,npjglo,nclass))
+ ALLOCATE ( dtrpus(npiglo,npjglo,nclass), dtrpvs(npiglo,npjglo,nclass))
+ ALLOCATE ( dheatrpsum(nclass), dsaltrpsum(nclass) )
+ ALLOCATE ( dheatallegcl(nclass), dsaltallegcl(nclass) )
+ ENDIF
+ !
+ ALLOCATE ( e1v(npiglo,npjglo),e3v(npiglo,npjglo) )
+ ALLOCATE ( e2u(npiglo,npjglo),e3u(npiglo,npjglo) )
+ !
+ ALLOCATE ( gphif(npiglo,npjglo) )
+ ALLOCATE ( glamf(npiglo,npjglo) )
+ ALLOCATE ( gdepw(npk) , tim(npt) )
+ !
+ ! read metrics and grid position
+ e1v(:,:) = getvar(cn_fhgr, cn_ve1v, 1, npiglo, npjglo)
+ e2u(:,:) = getvar(cn_fhgr, cn_ve2u, 1, npiglo, npjglo)
+
+ glamf(:,:) = getvar(cn_fhgr, cn_glamf, 1,npiglo, npjglo)
+ gphif(:,:) = getvar(cn_fhgr, cn_gphif, 1,npiglo, npjglo)
+
+ gdepw(:) = getvare3(cn_fzgr, cn_gdepw, npk)
+ e31d(:) = getvare3(cn_fzgr, cn_ve3t, npk) ! used only for full step
+
+ ! look for nearest level to imeter and setup ilev0 and ilev1 (t-index of class limit)
+ ik = 1
+ ilev0(1) = 1 ; ilev1(nclass) = npk-1 ! default value if nclass=1
+
+ IF ( lobc ) THEN
+ ! read u, v on OBC
+ IF ( l_zonal ) THEN ! (jpiglo,jpk)
+ zuobc(:,:)= getvarxz(cf_ufil, cn_vozocrtx, 1, npiglo, npk)
+ zvobc(:,:)= getvarxz(cf_vfil, cn_vomecrty, 1, npiglo, npk)
+ ENDIF
+ IF ( l_merid ) THEN ! (jpjglo,jpk)
+ zuobc(:,:)= getvaryz(cf_ufil, cn_vozocrtx, 1, npjglo, npk)
+ zvobc(:,:)= getvaryz(cf_vfil, cn_vomecrty, 1, npjglo, npk)
+ ENDIF
+ ENDIF
+
+ DO jclass = 1, nclass -1
+ DO WHILE ( gdepw(ik) < imeter(jclass) )
+ ik = ik +1
+ END DO
+
+ rd1 = ABS(gdepw(ik-1) - imeter(jclass) )
+ rd2 = ABS(gdepw(ik ) - imeter(jclass) )
+ IF ( rd2 < rd1 ) THEN
+ ilev1(jclass ) = ik - 1 ! t-levels index
+ ilev0(jclass+1) = ik
+ ELSE
+ ilev1(jclass ) = ik - 2 ! t-levels index
+ ilev0(jclass+1) = ik - 1
+ END IF
+ END DO
+
+ PRINT *, 'Limits : '
+ DO jclass = 1, nclass
+ PRINT *, ilev0(jclass),ilev1(jclass), gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
+ END DO
+
+ ! compute the transports at each grid cell
+ dtrpu (:,:,:)= 0.d0 ; dtrpv (:,:,:)= 0.d0 ! initialization to 0
+
+ IF ( lpm ) THEN
+ dtrpup(:,:,:)= 0.d0 ; dtrpvp(:,:,:)= 0.d0
+ dtrpum(:,:,:)= 0.d0 ; dtrpvm(:,:,:)= 0.d0
+ ENDIF
+ IF ( lheat ) THEN
+ dtrput(:,:,:)= 0.d0 ; dtrpvt(:,:,:)= 0.d0
+ dtrpus(:,:,:)= 0.d0 ; dtrpvs(:,:,:)= 0.d0
+ ENDIF
+
+ DO jclass = 1, nclass
+ DO jk = ilev0(jclass),ilev1(jclass)
+ PRINT *,'level ',jk
+ ! Get velocities, temperature and salinity fluxes at jk
+ IF ( ltest ) THEN
+ zu (:,:) = udum ; zv (:,:) = vdum
+ IF (lheat) THEN
+ zut(:,:) = udum ; zvt(:,:) = vdum
+ zus(:,:) = udum ; zvs(:,:) = vdum
+ ENDIF
+ ELSEIF ( lobc ) THEN
+ IF ( l_zonal ) THEN ; zu(:,1)=zuobc(:,jk) ; zv(:,1)=zvobc(:,jk)
+ ELSE IF ( l_merid ) THEN ; zu(1,:)=zuobc(:,jk) ; zv(1,:)=zvobc(:,jk) ; ENDIF
+ ELSE
+ zu (:,:) = getvar(cf_ufil, cn_vozocrtx, jk, npiglo, npjglo, ktime=itime)
+ zv (:,:) = getvar(cf_vfil, cn_vomecrty, jk, npiglo, npjglo, ktime=itime)
+ IF (lheat) THEN
+ zut(:,:) = getvar(cf_tfil, cn_vozout, jk, npiglo, npjglo, ktime=itime)
+ zvt(:,:) = getvar(cf_tfil, cn_vomevt, jk, npiglo, npjglo, ktime=itime)
+ zus(:,:) = getvar(cf_tfil, cn_vozous, jk, npiglo, npjglo, ktime=itime)
+ zvs(:,:) = getvar(cf_tfil, cn_vomevs, jk, npiglo, npjglo, ktime=itime)
+ ENDIF
+ ENDIF
+
+ ! get e3u, e3v at level jk
+ IF ( lfull ) THEN
+ e3v(:,:) = e31d(jk)
+ e3u(:,:) = e31d(jk)
+ ELSE
+ e3v(:,:) = getvar(cn_fzgr, 'e3v_ps', jk, npiglo, npjglo, ldiom=.TRUE.)
+ e3u(:,:) = getvar(cn_fzgr, 'e3u_ps', jk, npiglo, npjglo, ldiom=.TRUE.)
+ ENDIF
+
+ dwku (:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
+ dwkv (:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
+
+ IF ( lpm ) THEN
+ dwkup = 0.d0 ; dwkum = 0.d0
+ WHERE ( zu >= 0. )
+ dwkup(:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
+ ELSEWHERE
+ dwkum(:,:) = zu (:,:)*e2u(:,:)*e3u(:,:)*1.d0
+ END WHERE
+
+ dwkvp = 0.d0 ; dwkvm = 0.d0
+ WHERE ( zv >= 0. )
+ dwkvp(:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
+ ELSEWHERE
+ dwkvm(:,:) = zv (:,:)*e1v(:,:)*e3v(:,:)*1.d0
+ END WHERE
+ ENDIF
+
+ IF ( lheat ) THEN
+ dwkut(:,:) = zut(:,:)*e2u(:,:)*e3u(:,:)*1.d0
+ dwkvt(:,:) = zvt(:,:)*e1v(:,:)*e3v(:,:)*1.d0
+ dwkus(:,:) = zus(:,:)*e2u(:,:)*e3u(:,:)*1.d0
+ dwkvs(:,:) = zvs(:,:)*e1v(:,:)*e3v(:,:)*1.d0
+ ENDIF
+
+ ! integrates vertically
+ dtrpu (:,:,jclass) = dtrpu (:,:,jclass) + dwku (:,:)
+ dtrpv (:,:,jclass) = dtrpv (:,:,jclass) + dwkv (:,:)
+
+ IF ( lpm ) THEN
+ dtrpup(:,:,jclass) = dtrpup(:,:,jclass) + dwkup(:,:)
+ dtrpvp(:,:,jclass) = dtrpvp(:,:,jclass) + dwkvp(:,:)
+ dtrpum(:,:,jclass) = dtrpum(:,:,jclass) + dwkum(:,:)
+ dtrpvm(:,:,jclass) = dtrpvm(:,:,jclass) + dwkvm(:,:)
+ ENDIF
+
+ IF ( lheat ) THEN
+ dtrput(:,:,jclass) = dtrput(:,:,jclass) + dwkut(:,:) * rau0*rcp
+ dtrpvt(:,:,jclass) = dtrpvt(:,:,jclass) + dwkvt(:,:) * rau0*rcp
+ dtrpus(:,:,jclass) = dtrpus(:,:,jclass) + dwkus(:,:)
+ dtrpvs(:,:,jclass) = dtrpvs(:,:,jclass) + dwkvs(:,:)
+ ENDIF
+
+ END DO ! loop to next level
+ END DO ! next class
+
+ OPEN(numout,FILE=cf_out)
+ ! also dump the results on txt files without any comments, some users like it !
+ OPEN(numvtrp,FILE=cf_vtrp)
+ IF ( lheat ) THEN
+ OPEN(numhtrp,FILE=cf_htrp) ; OPEN(numstrp,FILE=cf_strp)
+ ENDIF
+
+ !################################################################################
+ ! enter interactive part
+ !################################################################################
+ ! initialize all legs arrays and variable to 0
+ dvolalleg = 0.d0 ; dvolallegcl(:) = 0.d0
+ IF ( lpm ) THEN
+ dvolallegp = 0.d0 ; dvolallegclp(:) = 0.d0
+ dvolallegm = 0.d0 ; dvolallegclm(:) = 0.d0
+ ENDIF
+ IF ( lheat ) THEN
+ dheatalleg = 0.d0 ; dheatallegcl(:) = 0.d0
+ dsaltalleg = 0.d0 ; dsaltallegcl(:) = 0.d0
+ ENDIF
+ DO
+ PRINT *, ' Give name of section (EOF to finish)'
+ READ(*,'(a)') cline
+ ii = 0
+ cldumt(:) = 'none'
+ ipos = index(cline,' ')
+ DO WHILE ( ipos > 1 )
+ ii = ii + 1
+ cldumt(ii) = cline(1:ipos - 1 )
+ cline = TRIM ( cline(ipos+1:) )
+ ipos = index( cline,' ' )
+ IF ( ii >= 3 ) EXIT
+ END DO
+ csection = TRIM(cldumt(1) )
+ cvarname = TRIM(cldumt(2) )
+ clongname = TRIM(cldumt(3) )
+
+ IF (TRIM(csection) == 'EOF' ) THEN
+ CLOSE(numout) ; CLOSE(numvtrp)
+ IF ( lheat ) THEN
+ CLOSE(numhtrp) ; CLOSE(numstrp)
+ ENDIF
+ EXIT ! infinite DO loop
+ ENDIF
+
+ ! create output fileset
+ CALL set_typvar( stypvar, csection, cvarname, clongname )
+ cf_outnc = TRIM(csection)//'_transports.nc'
+ ncout = create (cf_outnc, 'none', ikx, iky, nclass, cdep='depth_class')
+ ierr = createvar (ncout, stypvar, nvarout, ipk, id_varout, cdglobal=TRIM(cglobal) )
+ ierr = putheadervar(ncout, cf_ufil, ikx, iky, nclass, pnavlon=rdum, pnavlat=rdum, pdep=rclass )
+ tim = getvar1d (cf_ufil, cn_vtimec, npt )
+ ierr = putvar1d (ncout, tim, npt, 'T')
+
+ PRINT *, ' Give iimin, iimax, ijmin, ijmax '
+ READ(*,*) iimin, iimax, ijmin, ijmax
+ !! Find the broken line between P1 (iimin,ijmin) and P2 (iimax, ijmax)
+ ! ... Initialization
+ ii0 = iimin ; ij0 = ijmin ; ii1 = iimax ; ij1 = ijmax
+ rxi0 = ii0 ; ryj0 = ij0 ; rxi1 = ii1 ; ryj1 = ij1
+
+ ! compute direction of integrations and signs
+ !The transport across the section is the dot product of
+ !integral(line){(Mx,My)*dS}
+ !Mx=integral(u*dz) My=integral(v*dz)) and dS=(dy,-dx)}
+
+ !By defining the direction of the integration as
+ idirx = SIGN(1,ii1-ii0) !positive to the east or if ii1=ii0
+ idiry = SIGN(1,ij1-ij0) !positive to the north or if ij1=ij0
+
+ !Then dS=(e2u*idiry,-e1v*idirx)
+ !This will produce the following sign convention:
+ ! West-to-est line (dx>0, dy=0)=> -My*dx (-ve for a northward flow)
+ ! South-to-north (dy>0, dx=0)=> Mx*dy (+ve for an eastward flow)
+ norm_u = idiry
+ norm_v = -idirx
+
+ ! .. Compute equation: ryj = aj rxi + bj [valid in the (i,j) plane]
+ IF ( (rxi1 -rxi0) /= 0 ) THEN
+ aj = (ryj1 - ryj0 ) / (rxi1 -rxi0)
+ bj = ryj0 - aj * rxi0
+ ELSE
+ aj = 10000. ! flag value
+ bj = 0.
+ END IF
+
+ ! .. Compute equation: rxi = ai ryj + bi [valid in the (i,j) plane]
+ IF ( (ryj1 -ryj0) /= 0 ) THEN
+ ai = (rxi1 - rxi0 ) / ( ryj1 -ryj0 )
+ bi = rxi0 - ai * ryj0
+ ELSE
+ ai = 10000. ! flag value
+ bi = 0.
+ END IF
+
+ ! .. Compute the integer pathway: a succession of F points
+ np=0
+ ! .. Chose the strait line with the smallest slope
+ IF (ABS(aj) <= 1 ) THEN
+ ! ... Here, the best line is y(x)
+ ! ... If ii1 < ii0 swap points [ always describe section from left to right ]
+ IF (ii1 < ii0 ) THEN
+ iitmp = ii0 ; ijtmp = ij0
+ ii0 = ii1 ; ij0 = ij1
+ ii1 = iitmp ; ij1 = ijtmp
+ END IF
+
+ ! iist,ijst is the grid offset to pass from F point to either U/V point
+ IF ( ij1 >= ij0 ) THEN ! line heading NE
+ iist = 1 ; ijst = 1
+ ELSE ! line heading SE
+ iist = 1 ; ijst = 0
+ END IF
+
+ ! ... compute the nearest ji point on the line crossing at ji
+ DO ji=ii0, ii1
+ np=np+1
+ IF (np > jpseg) STOP 'np > jpseg !'
+ ij=NINT(aj*ji + bj )
+ yypt(np) = CMPLX(ji,ij)
+ END DO
+ ELSE
+ ! ... Here, the best line is x(y)
+ ! ... If ij1 < ij0 swap points [ always describe section from bottom to top ]
+ IF (ij1 < ij0 ) THEN
+ iitmp = ii0 ; ijtmp = ij0
+ ii0 = ii1 ; ij0 = ij1
+ ii1 = iitmp ; ij1 = ijtmp
+ END IF
+
+ ! iist,ijst is the grid offset to pass from F point to either U/V point
+ IF ( ii1 >= ii0 ) THEN
+ iist = 1 ; ijst = 1
+ ELSE
+ iist = 0 ; ijst = 1
+ END IF
+
+ ! ... compute the nearest ji point on the line crossing at jj
+ DO jj=ij0,ij1
+ np=np+1
+ IF (np > jpseg) STOP 'np > jpseg !'
+ ii=NINT(ai*jj + bi)
+ yypt(np) = CMPLX(ii,jj)
+ END DO
+ END IF
+
+ !!
+ !! Look for intermediate points to be added.
+ ! .. The final positions are saved in rxx,ryy
+ rxx(1) = REAL(yypt(1))
+ ryy(1) = IMAG(yypt(1))
+ nn = 1
+
+ DO jk=2,np
+ ! .. distance between 2 neighbour points
+ rd=ABS(yypt(jk)-yypt(jk-1))
+ ! .. intermediate points required if rd > 1
+ IF ( rd > 1 ) THEN
+ CALL interm_pt(yypt, jk, ai, bi, aj, bj, yypti)
+ nn=nn+1
+ IF (nn > jpseg) STOP 'nn>jpseg !'
+ rxx(nn) = REAL(yypti)
+ ryy(nn) = IMAG(yypti)
+ END IF
+ nn=nn+1
+ IF (nn > jpseg) STOP 'nn>jpseg !'
+ rxx(nn) = REAL(yypt(jk))
+ ryy(nn) = IMAG(yypt(jk))
+ END DO
+ ! record longitude and latitude of initial en endind point of the section
+ gla (1) = glamf( INT(rxx(1)), INT(ryy(1)) )
+ gphi(1) = gphif( INT(rxx(1)), INT(ryy(1)) )
+ gla (2) = glamf( INT(rxx(nn)), INT(ryy(nn)) )
+ gphi(2) = gphif( INT(rxx(nn)), INT(ryy(nn)) )
+
+ ! Now extract the transport through a section
+ ! ... Check whether we need a u velocity or a v velocity
+ ! Think that the points are f-points and delimit either a U segment
+ ! or a V segment (iist and ijst are set in order to look for the correct
+ ! velocity point on the C-grid
+ PRINT *, TRIM(csection)
+ PRINT *, 'IMIN IMAX JMIN JMAX', iimin, iimax, ijmin, ijmax
+ WRITE(numout,*) '% Transport along a section by levels' ,TRIM(csection)
+ WRITE(numout,*) '% ---- IMIN IMAX JMIN JMAX'
+
+ dvoltrpbrtp = 0.d0
+ dvoltrpbrtpp = 0.d0
+ dvoltrpbrtpm = 0.d0
+ dheatrpbrtp = 0.d0
+ dsaltrpbrtp = 0.d0
+ DO jclass=1,nclass
+ dvoltrpsum(jclass) = 0.d0
+ IF ( lpm ) THEN
+ dvoltrpsump(jclass) = 0.d0
+ dvoltrpsumm(jclass) = 0.d0
+ ENDIF
+ IF ( lheat ) THEN
+ dheatrpsum(jclass) = 0.d0
+ dsaltrpsum(jclass) = 0.d0
+ ENDIF
+
+ ! segment jseg is a line between (rxx(jseg),ryy(jseg)) and (rxx(jseg+1),ryy(jseg+1))
+ DO jseg = 1, nn-1
+ ii0=rxx(jseg)
+ ij0=ryy(jseg)
+ IF ( rxx(jseg) == rxx(jseg+1) ) THEN ! meridional segment, use U velocity
+ dvoltrp(jseg)= dtrpu (ii0,ij0+ijst,jclass)*norm_u
+
+ IF ( lpm ) THEN
+ dvoltrpp(jseg)= dtrpup(ii0,ij0+ijst,jclass)*norm_u
+ dvoltrpm(jseg)= dtrpum(ii0,ij0+ijst,jclass)*norm_u
+ ENDIF
+
+ IF ( lheat ) THEN
+ dheatrp(jseg)= dtrput(ii0,ij0+ijst,jclass)*norm_u
+ dsaltrp(jseg)= dtrpus(ii0,ij0+ijst,jclass)*norm_u
+ ENDIF
+ ELSE IF ( ryy(jseg) == ryy(jseg+1) ) THEN ! zonal segment, use V velocity
+ dvoltrp(jseg)=dtrpv (ii0+iist,ij0,jclass)*norm_v
+
+ IF ( lpm ) THEN
+ dvoltrpp(jseg)=dtrpvp(ii0+iist,ij0,jclass)*norm_v
+ dvoltrpm(jseg)=dtrpvm(ii0+iist,ij0,jclass)*norm_v
+ ENDIF
+
+ IF ( lheat ) THEN
+ dheatrp(jseg)=dtrpvt(ii0+iist,ij0,jclass)*norm_v
+ dsaltrp(jseg)=dtrpvs(ii0+iist,ij0,jclass)*norm_v
+ ENDIF
+ ELSE
+ PRINT *,' ERROR :', rxx(jseg),ryy(jseg),rxx(jseg+1),ryy(jseg+1) ! likely to never happen !
+ END IF
+ dvoltrpsum(jclass) = dvoltrpsum(jclass) + dvoltrp(jseg)
+ IF ( lpm ) THEN
+ dvoltrpsump(jclass) = dvoltrpsump(jclass) + dvoltrpp(jseg)
+ dvoltrpsumm(jclass) = dvoltrpsumm(jclass) + dvoltrpm(jseg)
+ ENDIF
+ IF ( lheat ) THEN
+ dheatrpsum(jclass) = dheatrpsum(jclass) + dheatrp(jseg)
+ dsaltrpsum(jclass) = dsaltrpsum(jclass) + dsaltrp(jseg)
+ ENDIF
+ END DO ! next segment
+
+ ! Ascii outputs :
+ IF (jclass == 1 ) THEN ! print header when it is the first class
+ PRINT '(a,2f8.2,a,2f8.2)', 'FROM (LON LAT): ', gla(1),gphi(1),' TO (LON LAT) ', gla(2), gphi(2)
+ WRITE(numout,*) '% ---- LONmin LATmin LONmax LATmax'
+ WRITE(numout,*) '% Top(m) Bottom(m) MassTrans(Sv) HeatTrans(PW) SaltTrans(kt/s)'
+ WRITE(numout,*) 0 ,iimin, iimax, ijmin, ijmax
+ WRITE(numout,9003) 0. ,gla(1),gphi(1), gla(2), gphi(2)
+ ENDIF
+
+ PRINT *, gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
+ PRINT *, ' Mass transport : ', dvoltrpsum(jclass)/1.e6,' SV'
+ WRITE(numvtrp,'(e12.6)') dvoltrpsum(jclass)
+ IF ( lpm ) THEN
+ PRINT *, ' Positive Mass transport : ', dvoltrpsump(jclass)/1.e6,' SV'
+ PRINT *, ' Negative Mass transport : ', dvoltrpsumm(jclass)/1.e6,' SV'
+ WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), &
+ & dvoltrpsum(jclass)/1.e6, dvoltrpsump(jclass)/1.e6, dvoltrpsumm(jclass)/1.e6
+ WRITE(numvtrp,'(e12.6)') dvoltrpsump(jclass)
+ WRITE(numvtrp,'(e12.6)') dvoltrpsumm(jclass)
+ ENDIF
+
+ IF ( lheat ) THEN
+ PRINT *, ' Heat transport : ', dheatrpsum(jclass)/1.e15,' PW'
+ PRINT *, ' Salt transport : ', dsaltrpsum(jclass)/1.e6,' kT/s'
+ WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), &
+ & dvoltrpsum(jclass)/1.e6, dheatrpsum(jclass)/1.e15, dsaltrpsum(jclass)/1.e6
+ WRITE(numhtrp,'(e12.6)') dheatrpsum(jclass)
+ WRITE(numstrp,'(e12.6)') dsaltrpsum(jclass)
+ ELSE
+ IF ( .NOT. lpm ) WRITE(numout,9002) gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1), dvoltrpsum(jclass)/1.e6
+ ENDIF
+
+ ! netcdf output
+ IF ( nclass > 1 ) THEN
+ rdum(1,1) = REAL(dvoltrpsum(jclass)/1.e6)
+ ierr = putvar(ncout,id_varout(ivtrpcl), rdum, jclass, 1, 1, ktime=itime )
+ IF ( lpm ) THEN
+ rdum(1,1) = REAL(dvoltrpsump(jclass)/1.e6)
+ ierr = putvar(ncout,id_varout(iptrpcl), rdum, jclass, 1, 1, ktime=itime )
+ rdum(1,1) = REAL(dvoltrpsumm(jclass)/1.e6)
+ ierr = putvar(ncout,id_varout(imtrpcl), rdum, jclass, 1, 1, ktime=itime )
+ ENDIF
+ IF ( lheat ) THEN
+ rdum(1,1) = REAL(dheatrpsum(jclass)/1.e15)
+ ierr = putvar(ncout,id_varout(ihtrpcl), rdum, jclass, 1, 1, ktime=itime )
+ rdum(1,1) = REAL(dsaltrpsum(jclass)/1.e6)
+ ierr = putvar(ncout,id_varout(istrpcl), rdum, jclass, 1, 1, ktime=itime )
+ ENDIF
+ ENDIF
+ rdum(1,1) = REAL(gdepw(ilev0(jclass)))
+ ierr = putvar(ncout,id_varout(itop), rdum, jclass, 1, 1, ktime=itime )
+ rdum(1,1) = REAL(gdepw(ilev1(jclass)+1))
+ ierr = putvar(ncout,id_varout(ibot), rdum, jclass, 1, 1, ktime=itime )
+
+ dvoltrpbrtp = dvoltrpbrtp + dvoltrpsum(jclass)
+ IF ( lpm ) THEN
+ dvoltrpbrtpp = dvoltrpbrtpp + dvoltrpsump(jclass)
+ dvoltrpbrtpm = dvoltrpbrtpm + dvoltrpsumm(jclass)
+ ENDIF
+ IF ( lheat) THEN
+ dheatrpbrtp = dheatrpbrtp + dheatrpsum(jclass)
+ dsaltrpbrtp = dsaltrpbrtp + dsaltrpsum(jclass)
+ ENDIF
+ ! save sum over legs
+ dvolallegcl(jclass) = dvolallegcl(jclass) + dvoltrpsum(jclass)
+ IF ( lpm ) THEN
+ dvolallegclp(jclass) = dvolallegclp(jclass) + dvoltrpsump(jclass)
+ dvolallegclm(jclass) = dvolallegclm(jclass) + dvoltrpsumm(jclass)
+ ENDIF
+ IF ( lheat ) THEN
+ dheatallegcl(jclass) = dheatallegcl(jclass) + dheatrpsum(jclass)
+ dsaltallegcl(jclass) = dsaltallegcl(jclass) + dsaltrpsum(jclass)
+ ENDIF
+ END DO ! next class
+ ! save sum over legs
+ dvolalleg = dvolalleg + dvoltrpbrtp
+ IF ( lpm ) THEN
+ dvolallegp = dvolallegp + dvoltrpbrtpp
+ dvolallegm = dvolallegm + dvoltrpbrtpm
+ ENDIF
+ IF ( lheat ) THEN
+ dheatalleg = dheatalleg + dheatrpbrtp
+ dsaltalleg = dsaltalleg + dsaltrpbrtp
+ ENDIF
+
+ IF ( nclass > 1 ) THEN
+ PRINT *, ' ====================================================='
+ PRINT *, ' total Mass transport : ', dvoltrpbrtp/1.e6,' SV'
+ IF ( lpm ) THEN
+ PRINT *, ' total positive transport : ', dvoltrpbrtpp/1.e6,' SV'
+ PRINT *, ' total negative transport : ', dvoltrpbrtpm/1.e6,' SV'
+ ENDIF
+ IF ( lheat ) THEN
+ PRINT *, ' total Heat transport : ', dheatrpbrtp/1.e15,' PW'
+ PRINT *, ' total Salt transport : ', dsaltrpbrtp/1.e6,' kT/s'
+ ENDIF
+ ENDIF
+ ierr = putvar0d(ncout,id_varout(ivtrp), REAL(dvoltrpbrtp/1.e6) )
+ IF ( lpm ) THEN
+ ierr = putvar0d(ncout,id_varout(iptrp), REAL(dvoltrpbrtpp/1.e6) )
+ ierr = putvar0d(ncout,id_varout(imtrp), REAL(dvoltrpbrtpm/1.e6) )
+ ENDIF
+ IF ( lheat ) THEN
+ ierr = putvar0d(ncout,id_varout(ihtrp), REAL(dheatrpbrtp/1.e15) )
+ ierr = putvar0d(ncout,id_varout(istrp), REAL(dsaltrpbrtp/1.e6 ) )
+ ENDIF
+ ierr = putvar0d(ncout,id_varout(ilonmin), REAL(gla(1)) )
+ ierr = putvar0d(ncout,id_varout(ilonmax), REAL(gla(2)) )
+ ierr = putvar0d(ncout,id_varout(ilatmin), REAL(gphi(1)) )
+ ierr = putvar0d(ncout,id_varout(ilatmax), REAL(gphi(2)) )
+ ierr = closeout(ncout)
+ END DO ! infinite loop : gets out when input is EOF
+
+
+ PRINT *,' '
+ PRINT *,' Overall transports (sum of all legs done so far)'
+ DO jclass = 1, nclass
+ PRINT *, gdepw(ilev0(jclass)), gdepw(ilev1(jclass)+1)
+ PRINT *, ' Mass transport : ', dvolallegcl(jclass)/1.e6,' SV'
+ IF ( lpm ) THEN
+ PRINT *, ' Positive Mass transport : ', dvolallegclp(jclass)/1.e6,' SV'
+ PRINT *, ' Negative Mass transport : ', dvolallegclm(jclass)/1.e6,' SV'
+ ENDIF
+
+ IF ( lheat ) THEN
+ PRINT *, ' Heat transport : ', dheatallegcl(jclass)/1.e15,' PW'
+ PRINT *, ' Salt transport : ', dsaltallegcl(jclass)/1.e6,' kT/s'
+ ENDIF
+ ENDDO
+
+ IF ( nclass > 1 ) THEN
+ PRINT *, ' ====================================================='
+ PRINT *, ' Mass transport : ', dvolalleg/1.e6,' SV'
+ IF ( lpm ) THEN
+ PRINT *, ' positive transport : ', dvolallegp/1.e6,' SV'
+ PRINT *, ' negative transport : ', dvolallegm/1.e6,' SV'
+ ENDIF
+ IF ( lheat ) THEN
+ PRINT *, ' heat transport : ', dheatalleg/1.e15,' PW'
+ PRINT *, ' salt transport : ', dsaltalleg/1.e6,' kT/s'
+ ENDIF
+ ENDIF
+
+
+9000 FORMAT(I4,6(f9.3,f8.4))
+9001 FORMAT(I4,6(f9.2,f9.3))
+9002 FORMAT(f9.0,f9.0,f9.2,f9.2,f9.2)
+9003 FORMAT(f9.2,f9.2,f9.2,f9.2,f9.2)
+
+ CONTAINS
+ SUBROUTINE set_typvar ( sd_typvar, cdsection, cdvarname, cdlongname )
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE set_typvar ***
+ !!
+ !! ** Purpose : Initialize typvar structure for netcdfoutput at a given section
+ !!
+ !! ** Method : use varname longname to suffix variable name and attributes
+ !! If varname and/or logname are not given (ie 'none') take
+ !! standard default names
+ !! Netcdf id for variables are passed as global variables
+ !!----------------------------------------------------------------------
+ TYPE(variable), DIMENSION(:), INTENT(out) :: sd_typvar ! structure of output
+ CHARACTER(LEN=*), INTENT(in ) :: cdsection
+ CHARACTER(LEN=*), INTENT(in ) :: cdvarname
+ CHARACTER(LEN=*), INTENT(in ) :: cdlongname
+ !!
+ INTEGER(KIND=4) :: ivar
+ CHARACTER(LEN=255) :: csuffixvarnam
+ CHARACTER(LEN=255) :: cprefixlongnam
+ !!----------------------------------------------------------------------
+ ! set suffixes according to variable/longname
+ IF ( cdvarname /= 'none' ) THEN
+ csuffixvarnam = '_'//TRIM(cdvarname)
+ ELSE
+ csuffixvarnam = ''
+ ENDIF
+
+ IF ( cdlongname /= 'none' ) THEN
+ cprefixlongnam = TRIM(cdlongname)//'_'
+ ELSE
+ cprefixlongnam = ''
+ ENDIF
+
+ ! set common values
+ sd_typvar%rmissing_value=99999.
+ sd_typvar%scale_factor= 1.
+ sd_typvar%add_offset= 0.
+ sd_typvar%savelog10= 0.
+ sd_typvar%conline_operation='N/A'
+ sd_typvar%caxis='T'
+
+ ! set particular values for individual variables
+ ivar = 1 ; ivtrp = ivar
+ ipk(ivar) = 1
+ sd_typvar(ivar)%cname = 'vtrp'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'Sverdrup'
+ sd_typvar(ivar)%valid_min = -500.
+ sd_typvar(ivar)%valid_max = 500.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Volume_Transport'
+ sd_typvar(ivar)%cshort_name = 'vtrp'
+
+ IF ( lpm ) THEN
+ ivar = ivar + 1 ; iptrp = ivar ; imtrp = ivar+1
+ ipk(ivar) = 1 ; ipk(ivar+1) = 1
+ sd_typvar(ivar)%cname = 'ptrp'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'mtrp'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'Sverdrup' ; sd_typvar(ivar+1)%cunits = 'Sverdrup'
+ sd_typvar(ivar)%valid_min = -500. ; sd_typvar(ivar+1)%valid_min = -500.
+ sd_typvar(ivar)%valid_max = 500. ; sd_typvar(ivar+1)%valid_max = 500.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Positive_volume_transport' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'Negative_volume_transport'
+ sd_typvar(ivar)%cshort_name = 'ptrp' ; sd_typvar(ivar+1)%cshort_name = 'mtrp'
+ ivar = ivar + 1
+ ENDIF
+
+ IF ( lheat ) THEN
+ ivar = ivar + 1 ; ihtrp = ivar ; istrp = ivar+1
+ ipk(ivar) = 1 ; ipk(ivar+1) = 1
+ sd_typvar(ivar)%cname = 'htrp'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'strp'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'PW' ; sd_typvar(ivar+1)%cunits = 'kt/s'
+ sd_typvar(ivar)%valid_min = -1000. ; sd_typvar(ivar+1)%valid_min = -1000.
+ sd_typvar(ivar)%valid_max = 1000. ; sd_typvar(ivar+1)%valid_max = 1000.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Heat_Transport' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'Salt_Transport'
+ sd_typvar(ivar)%cshort_name = 'htrp' ; sd_typvar(ivar+1)%cshort_name = 'strp'
+ ivar = ivar + 1
+ ENDIF
+
+ ivar = ivar + 1 ; ilonmin = ivar ; ilonmax = ivar+1
+ ipk(ivar) = 1 ; ipk(ivar+1) = 1
+ sd_typvar(ivar)%cname = 'lonmin'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'lonmax'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'deg' ; sd_typvar(ivar+1)%cunits = 'deg'
+ sd_typvar(ivar)%valid_min = -180. ; sd_typvar(ivar+1)%valid_min = -180.
+ sd_typvar(ivar)%valid_max = 180. ; sd_typvar(ivar+1)%valid_max = 180.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'begin_longitude' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'end_longitude'
+ sd_typvar(ivar)%cshort_name = 'lonmin' ; sd_typvar(ivar+1)%cshort_name = 'lonmax'
+ ivar = ivar + 1
+
+ ivar = ivar + 1 ; ilatmin = ivar ; ilatmax = ivar+1
+ ipk(ivar) = 1 ; ipk(ivar+1) = 1
+ sd_typvar(ivar)%cname = 'latmin'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'latmax'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'deg' ; sd_typvar(ivar+1)%cunits = 'deg'
+ sd_typvar(ivar)%valid_min = -90. ; sd_typvar(ivar+1)%valid_min = -90.
+ sd_typvar(ivar)%valid_max = 90. ; sd_typvar(ivar+1)%valid_max = 90.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'begin_latitude' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'end_latitude'
+ sd_typvar(ivar)%cshort_name = 'latmin' ; sd_typvar(ivar+1)%cshort_name = 'latmax'
+ ivar = ivar + 1
+
+ ivar = ivar + 1 ; itop = ivar ; ibot = ivar+1
+ ipk(ivar) = nclass ; ipk(ivar+1) = nclass
+ sd_typvar(ivar)%cname = 'top' ; sd_typvar(ivar+1)%cname = 'bottom'
+ sd_typvar(ivar)%cunits = 'meters' ; sd_typvar(ivar+1)%cunits = 'meters'
+ sd_typvar(ivar)%valid_min = 0. ; sd_typvar(ivar+1)%valid_min = 0.
+ sd_typvar(ivar)%valid_max = 10000. ; sd_typvar(ivar+1)%valid_max = 10000.
+ sd_typvar(ivar)%clong_name = 'class_min_depth' ; sd_typvar(ivar+1)%clong_name = 'class_max_depth'
+ sd_typvar(ivar)%cshort_name = 'top' ; sd_typvar(ivar+1)%cshort_name = 'bottom'
+ ivar = ivar + 1
+
+ ivtrpcl = -1 ; ihtrpcl = -1 ; istrpcl = -1
+ IF ( nclass > 1 ) THEN ! define additional variable for vertical profile of transport (per class)
+ ivar = ivar + 1 ; ivtrpcl = ivar
+ ipk(ivar) = nclass
+ sd_typvar(ivar)%cname = 'vtrp_dep'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'SV'
+ sd_typvar(ivar)%valid_min = 0.
+ sd_typvar(ivar)%valid_max = 10000.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Volume_Transport_per_class'
+ sd_typvar(ivar)%cshort_name = 'vtrp_dep'
+
+ IF ( lpm ) THEN
+ ivar = ivar + 1 ; iptrpcl = ivar ; imtrpcl = ivar+1
+ ipk(ivar) = nclass ; ipk(ivar+1) = nclass
+ sd_typvar(ivar)%cname = 'ptrp_dep'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'mtrp_dep'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'SV' ; sd_typvar(ivar+1)%cunits = 'SV'
+ sd_typvar(ivar)%valid_min = -500. ; sd_typvar(ivar+1)%valid_min = -500.
+ sd_typvar(ivar)%valid_max = 500. ; sd_typvar(ivar+1)%valid_max = 500.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Positive_trp_per_class' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'Negative_trp_per_class'
+ sd_typvar(ivar)%cshort_name = 'ptrp_dep' ; sd_typvar(ivar+1)%cshort_name = 'mtrp_dep'
+ ivar = ivar + 1
+ ENDIF
+
+ IF ( lheat ) THEN
+ ivar = ivar + 1 ; ihtrpcl = ivar ; istrpcl = ivar+1
+ ipk(ivar) = nclass ; ipk(ivar+1) = nclass
+ sd_typvar(ivar)%cname = 'htrp_dep'//TRIM(csuffixvarnam) ; sd_typvar(ivar+1)%cname = 'strp_dep'//TRIM(csuffixvarnam)
+ sd_typvar(ivar)%cunits = 'PW' ; sd_typvar(ivar+1)%cunits = 'kt/s'
+ sd_typvar(ivar)%valid_min = -1000. ; sd_typvar(ivar+1)%valid_min = -1000.
+ sd_typvar(ivar)%valid_max = 1000. ; sd_typvar(ivar+1)%valid_max = 1000.
+ sd_typvar(ivar)%clong_name = TRIM(cprefixlongnam)//'Heat_Transport_per_class' ; sd_typvar(ivar+1)%clong_name = TRIM(cprefixlongnam)//'Salt_Transport_per_class'
+ sd_typvar(ivar)%cshort_name = 'htrp_dep' ; sd_typvar(ivar+1)%cshort_name = 'strp_dep'
+ ivar = ivar + 1
+ ENDIF
+ ENDIF
+
+ END SUBROUTINE set_typvar
+ SUBROUTINE interm_pt (ydpt, kk, pai, pbi, paj, pbj, ydpti)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE nterm_pt ***
+ !!
+ !! ** Purpose : Find the best intermediate points on a pathway.
+ !!
+ !! ** Method : ydpt : complex vector of the positions of the nearest points
+ !! kk : current working index
+ !! pai, pbi : slope and original ordinate of x(y)
+ !! paj, pbj : slope and original ordinate of y(x)
+ !! ydpti : Complex holding the position of intermediate point
+ !!
+ !! ** Reference : 19/07/1999 : J.M. Molines in Clipper
+ !!----------------------------------------------------------------------
+ COMPLEX, DIMENSION(:), INTENT(in ) :: ydpt
+ COMPLEX, INTENT(out) :: ydpti
+ REAL(KIND=4), INTENT(in ) :: pai, pbi, paj, pbj
+ INTEGER(KIND=4), INTENT(in ) :: kk
+ ! ... local
+ COMPLEX :: ylptmp1, ylptmp2
+ REAL(KIND=4) :: za0, zb0
+ REAL(KIND=4) :: za1, zb1
+ REAL(KIND=4) :: zd1, zd2
+ REAL(KIND=4) :: zxm, zym
+ REAL(KIND=4) :: zxp, zyp
+ !!----------------------------------------------------------------------
+ ! ... Determines whether we use y(x) or x(y):
+ IF (ABS(paj) <= 1) THEN
+ ! ..... use y(x)
+ ! ... possible intermediate points:
+ ylptmp1=ydpt(kk-1)+(1.,0.) ! M1
+ ylptmp2=ydpt(kk-1)+CMPLX(0.,SIGN(1.,paj)) ! M2
+ !
+ ! ... M1 is the candidate point:
+ zxm=REAL(ylptmp1)
+ zym=IMAG(ylptmp1)
+ za0=paj
+ zb0=pbj
+ !
+ za1=-1./za0
+ zb1=zym - za1*zxm
+ ! ... P1 is the projection of M1 on the strait line
+ zxp=-(zb1-zb0)/(za1-za0)
+ zyp=za0*zxp + zb0
+ ! ... zd1 is the distance M1P1
+ zd1=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
+ !
+ ! ... M2 is the candidate point:
+ zxm=REAL(ylptmp2)
+ zym=IMAG(ylptmp2)
+ za1=-1./za0
+ zb1=zym - za1*zxm
+ ! ... P2 is the projection of M2 on the strait line
+ zxp=-(zb1-zb0)/(za1-za0)
+ zyp=za0*zxp + zb0
+ ! ... zd2 is the distance M2P2
+ zd2=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
+ ! ... chose the smallest (zd1,zd2)
+ IF (zd2 <= zd1) THEN
+ ydpti=ylptmp2 ! use M2
+ ELSE
+ ydpti=ylptmp1 ! use M1
+ END IF
+ !
+ ELSE
+ ! ... use x(y)
+ ! ... possible intermediate points:
+ ylptmp1=ydpt(kk-1)+CMPLX(SIGN(1.,pai),0.) ! M1
+ ylptmp2=ydpt(kk-1)+(0.,1.) ! M2
+ !
+ ! ... M1 is the candidate point:
+ zxm=REAL(ylptmp1)
+ zym=IMAG(ylptmp1)
+ za0=pai
+ zb0=pbi
+ !
+ za1=-1./za0
+ zb1=zxm - za1*zym
+ zyp=-(zb1-zb0)/(za1-za0)
+ zxp=za0*zyp + zb0
+ zd1=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
+ !
+ zxm=REAL(ylptmp2)
+ zym=IMAG(ylptmp2)
+ za1=-1./za0
+ zb1=zxm - za1*zym
+ zyp=-(zb1-zb0)/(za1-za0)
+ zxp=za0*zyp + zb0
+ zd2=(zxm-zxp)*(zxm-zxp) + (zym-zyp)*(zym-zyp)
+ IF (zd2 <= zd1) THEN
+ ydpti=ylptmp2
+ ELSE
+ ydpti=ylptmp1
+ END IF
+ END IF
+ END SUBROUTINE interm_pt
+
+ END PROGRAM cdftransport
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/cdftools.git
More information about the debian-science-commits
mailing list