[cdftools] 174/228: JMM +NF : add cdf2levitusgrid2d.f90 and changes in the Makefile
Alastair McKinstry
mckinstry at moszumanska.debian.org
Fri Jun 12 08:21:45 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 9734e4de5786d93ae426dcb1fc9f8490a585fbed
Author: molines <molines at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date: Wed Oct 31 20:15:11 2012 +0000
JMM +NF : add cdf2levitusgrid2d.f90 and changes in the Makefile
git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@622 1055176f-818a-41d9-83e1-73fbe5b947c5
---
Makefile | 3 +
cdf2levitusgrid2d.f90 | 507 ++++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 510 insertions(+)
diff --git a/Makefile b/Makefile
index 8daa49d..9118b02 100644
--- a/Makefile
+++ b/Makefile
@@ -306,6 +306,9 @@ cdfcoloc2: cdfio.o cdfcoloc2.f90
cdfcoloc3: cdfio.o cdfcoloc3.f90
$(F90) cdfcoloc3.f90 -o $(BINDIR)/cdfcoloc3 cdfio.o $(FFLAGS)
+cdf2levitusgrid2d: cdfio.o cdftools.o modutils.o cdf2levitusgrid2d.f90
+ $(F90) cdf2levitusgrid2d.f90 -o $(BINDIR)/cdf2levitusgrid2d cdfio.o modcdfnames.o cdftools.o modutils.o $(FFLAGS)
+
cdfstatcoord: cdfio.o cdfstatcoord.f90
$(F90) cdfstatcoord.f90 -o $(BINDIR)/cdfstatcoord cdfio.o modcdfnames.o $(FFLAGS)
diff --git a/cdf2levitusgrid2d.f90 b/cdf2levitusgrid2d.f90
new file mode 100644
index 0000000..e77b69e
--- /dev/null
+++ b/cdf2levitusgrid2d.f90
@@ -0,0 +1,507 @@
+PROGRAM cdf2levitusgrid2d
+ !!======================================================================
+ !! *** PROGRAM cdf2levitusgrid2d ***
+ !!=====================================================================
+ !! ** Purpose : remaps (bin) 2D high resolution (finer than 1x1 deg)
+ !! fields on Levitus 2D 1x1 deg grid
+ !!
+ !! ** Method : data surface averaging
+ !! It assumes that Levitus grid SW grid cell center
+ !! is 0.5W,89.5S
+ !!
+ !! History : 3.0 : 06/2012 : N. Ferry : Original code
+ !!----------------------------------------------------------------------
+ USE cdfio
+ USE modcdfnames
+ USE cdftools
+ USE modutils
+ !!----------------------------------------------------------------------
+ !! CDFTOOLS_3.0 , MEOM 2011
+ !! $Id: cdf2levitusgrid2d.f90 553 2011-10-17 13:46:56Z molines $
+ !! Copyright (c) 2011, J.-M. Molines
+ !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
+ !!----------------------------------------------------------------------
+ IMPLICIT NONE
+
+ INTEGER(KIND=4) :: ji, jj, jk, jt ! dummy loop index
+ INTEGER(KIND=4) :: jilev,jjlev ! dummy loop index
+ INTEGER(KIND=4) :: jvar, numvar0 ! dummy loop index
+ INTEGER(KIND=4) :: ii, ij ! array index (not loop)
+ INTEGER(KIND=4) :: iilev, ijlev ! array index (not loop)
+ INTEGER(KIND=4) :: npiglo, npjglo ! size of the domain
+ INTEGER(KIND=4) :: npk, npt ! size of the domain
+ INTEGER(KIND=4) :: npilev, npjlev ! size of the Levitus domain
+ INTEGER(KIND=4) :: narg, iargc, ijarg ! browse line
+ INTEGER(KIND=4) :: ncout ! ncid of output file
+ INTEGER(KIND=4) :: ierr ! error status
+ INTEGER(KIND=4) :: nvars ! number of variables in the input file
+ INTEGER(KIND=4) :: nvarsout ! number of variables in the output file
+ INTEGER(KIND=4) :: iimin, iimax, ijmin, ijmax
+ INTEGER(KIND=4) :: imethod=1 ! interpolation method
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipk, id_var ! levels and varid's of input vars
+ INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE :: ipkout, id_varout ! levels and varid's of output vars
+
+ REAL(KIND=4) :: zradius=120. ! Distance (km) for the search bubble (FHZ)
+! REAL(KIND=4) :: etobd ! Beta-plane distance (FHZ)
+ REAL(KIND=4) :: rlon1, rlon2, rlat1, rlat2, rpos
+ REAL(KIND=4) :: gphitmin
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e1t, e2t ! horizontal T metrics
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: z_in ! input field
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: z_fill ! output
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamt, gphit ! T longitude latitude
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: rlonlev, rlatlev ! Levitus grid longitude latitude
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zbt
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmask ! input mask
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmasklev ! output mask
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: tim ! time counter
+ REAL(KIND=4), DIMENSION(:), ALLOCATABLE :: gdept ! depth axis
+
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: d_out, d_n ! output field and weighting field
+
+ CHARACTER(LEN=256) :: cf_in ! input file name
+ CHARACTER(LEN=256) :: cf_out ! output file name ( output)
+ CHARACTER(LEN=256) :: cf_levitus_mask='levitus_mask.nc' ! Levitus mask filename
+ CHARACTER(LEN=256) :: cv_nam ! variable name
+ CHARACTER(LEN=256) :: cldum ! dummy string
+ CHARACTER(LEN=256) :: ctcalendar ! time attributes
+ CHARACTER(LEN=256) :: cttitle ! time attributes
+ CHARACTER(LEN=256) :: ctlong_name ! time attributes
+ CHARACTER(LEN=256) :: ctaxis ! time attributes
+ CHARACTER(LEN=256) :: ctunits ! time attributes
+ CHARACTER(LEN=256) :: cttime_origin ! time attributes
+
+ CHARACTER(LEN=256), DIMENSION(:), ALLOCATABLE :: cv_names ! array of var name
+ CHARACTER(LEN=6) :: ctyp ! 'fill' or 'smooth' for shapiro
+
+ TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvar ! input attributes
+ TYPE(variable), DIMENSION(:), ALLOCATABLE :: stypvarout ! output attributes
+
+ LOGICAL :: lchk ! missing files flag
+ LOGICAL :: ltest
+ !!----------------------------------------------------------------------
+ CALL ReadCdfNames()
+
+ narg = iargc()
+ IF ( narg < 3 ) THEN
+ PRINT *,' usage : cdf2levitusgrid2d IN-file OUT-file VAR-name2D'
+ PRINT *,' '
+ PRINT *,' PURPOSE :'
+ PRINT *,' remaps (bin) 2D high resolution (finer than 1x1 deg) '
+ PRINT *,' fields on Levitus 2D 1x1 deg grid '
+ PRINT *,' (does not work for vector fields) '
+ PRINT *,' It assumes that Levitus grid SW grid cell center '
+ PRINT *,' is (0.5W,89.5S) '
+ PRINT *,' '
+ PRINT *,' ARGUMENTS :'
+ PRINT *,' IN-file : netcdf input file '
+ PRINT *,' OUT-file : netcdf output file '
+ PRINT *,' VAR-name2D : input variable name for interpolation '
+ PRINT *,' '
+ PRINT *,' OPTIONS :'
+ PRINT *,' '
+ PRINT *,' REQUIRED FILES :'
+ PRINT *,' ',TRIM(cn_fhgr)
+ PRINT *,' ',TRIM(cf_levitus_mask)
+ PRINT *,' '
+ PRINT *,' OUTPUT : '
+ PRINT *,' netcdf file : name given as second argument'
+ PRINT *,' variables : 2d_var_name'
+ STOP
+ ENDIF
+
+ ijarg = 1
+ CALL getarg(ijarg, cf_in) ; ijarg = ijarg + 1
+ CALL getarg(ijarg, cf_out) ; ijarg = ijarg + 1
+ CALL getarg(ijarg, cv_nam) ; ijarg = ijarg + 1
+
+ lchk = chkfile (cn_fhgr)
+ lchk = chkfile (cf_levitus_mask) .OR. lchk
+ lchk = chkfile (cf_in) .OR. lchk
+ IF ( lchk ) STOP ! missing files
+
+ npiglo = getdim(cf_in,cn_x)
+ npjglo = getdim(cf_in,cn_y)
+ npk = getdim(cf_in,cn_z)
+ npt = getdim(cf_in,cn_t)
+ npilev = getdim(cf_levitus_mask,cn_x)
+ npjlev = getdim(cf_levitus_mask,cn_y)
+
+ nvars = getnvar(cf_in)
+ ALLOCATE (cv_names(nvars) )
+ ALLOCATE (stypvar(nvars) , stypvarout(1))
+ ALLOCATE (id_var(nvars), ipk(nvars), id_varout(1), ipkout(1))
+ ALLOCATE ( zbt(npiglo,npjglo) , z_in(npiglo,npjglo) )
+
+ ! get list of variable names and collect attributes in stypvar (optional)
+ cv_names(:) = getvarname(cf_in, nvars, stypvar)
+
+ id_var(:) = (/(jvar, jvar=1,nvars)/)
+
+ ! ipk gives the number of level or 0 if not a T[Z]YX variable
+ ipk(:) = getipk(cf_in, nvars)
+ WHERE( ipk == 0 ) cv_names='none'
+ stypvar(:)%cname = cv_names
+
+ ! select variables to output:
+ ii=1
+ DO jk=1,nvars
+ IF ( TRIM(cv_names(jk)) == TRIM(cv_nam) ) THEN
+ ipkout(ii) = ipk(jk)
+ stypvarout(ii) = stypvar(jk)
+ stypvarout(ji)%rmissing_value=getspval ( cf_in, TRIM(cv_nam) )
+ PRINT*, 'rmissing_value = ', stypvarout(ji)%rmissing_value
+ nvarsout = ii
+ numvar0 = jk
+ ENDIF
+ ENDDO
+ z_in(:,:) = getvar(cf_in, cv_names(numvar0), 1, npiglo, npjglo)
+
+ ! Allocate the memory
+ ALLOCATE ( e1t(npiglo,npjglo), e2t(npiglo,npjglo) )
+ ALLOCATE ( glamt(npiglo,npjglo), gphit(npiglo,npjglo) )
+ ALLOCATE ( d_out(npilev,npjlev) , d_n(npilev,npjlev) )
+ ALLOCATE ( tmask(npiglo,npjglo) , tmasklev(npilev,npjlev))
+ ALLOCATE ( rlonlev(npilev,npjlev), rlatlev(npilev,npjlev) )
+ ALLOCATE ( gdept(1), tim(npt) )
+
+ ! Read the metrics from the mesh_hgr file
+ e1t = getvar(cn_fhgr, cn_ve1t, 1, npiglo, npjglo)
+ e2t = getvar(cn_fhgr, cn_ve2t, 1, npiglo, npjglo)
+
+ ! and the coordinates from the mesh_hgr file
+ glamt = getvar(cn_fhgr, cn_glamt, 1, npiglo, npjglo)
+ gphit = getvar(cn_fhgr, cn_gphit, 1, npiglo, npjglo)
+ gphitmin = MINVAL(gphit(:,1))
+ WHERE ( glamt < 0. )
+ glamt = glamt + 360.
+ END WHERE
+
+ ! get the tmask from the byte_mask file
+ tmask(:,:) = getvar(cn_fmsk, 'tmask', 1, npiglo, npjglo)
+
+ ! get the longitude,latitude,mask from the input Levitus mask file
+ rlonlev(:,:) = getvar(cf_levitus_mask, 'nav_lon', 1, npilev, npjlev)
+ rlatlev(:,:) = getvar(cf_levitus_mask, 'nav_lat' , 1, npilev, npjlev)
+ tmasklev(:,:) = getvar(cf_levitus_mask, 'mask', 1, npilev, npjlev)
+
+ ! create output fileset
+ ncout = create (cf_out, cf_levitus_mask, npilev, npjlev, 0 )
+ ierr = createvar (ncout , stypvarout, 1, ipkout, id_varout )
+ ierr = putheadervar(ncout , 'dummy', npilev, npjlev, 0 , pnavlon=rlonlev, pnavlat=rlatlev )
+
+ tim = getvar1d(cf_in, cn_vtimec, npt )
+ ierr = putvar1d(ncout, tim, npt, 'T')
+ ierr = gettimeatt(cf_in, cn_vtimec, ctcalendar, cttitle, ctlong_name, ctaxis, ctunits, cttime_origin )
+ ierr = puttimeatt(ncout, cn_vtimec, ctcalendar, cttitle, ctlong_name, ctaxis, ctunits, cttime_origin )
+
+ zbt(:,:) = e1t(:,:) * e2t(:,:) ! for surface weighting
+
+ DO jt = 1, npt
+ PRINT *,'jt = ', jt
+ ! Compute spatial mean by bin
+ !-----------------------------
+ ! Perform bining of the input file on the Levitus grid.
+ ! Input area weighted values are summed up into a Levitus 1x1 bin
+ d_out(:,:) = 0.d0
+ d_n (:,:) = 0.d0
+ DO jj=1,npjglo
+ DO ji=1,npiglo
+ iilev = MIN( 360, INT( glamt(ji,jj) ) + 1)
+ ijlev = MIN (180 , INT( gphit(ji,jj) + 90. ) + 1)
+ !IF ( iilev < 1 .OR. iilev .GT. 360 ) print*, 'iilev, glamt = ',iilev,glamt(ji,jj)
+ !IF ( ijlev < 1 .OR. ijlev .GT. 180 ) print*, 'ijlev, gphit = ',ijlev,gphit(ji,jj)
+ IF ( z_in(ji,jj) /= stypvarout(1)%rmissing_value ) THEN
+ d_out(iilev,ijlev) = d_out(iilev,ijlev) + (z_in(ji,jj)*tmask(ji,jj))*zbt(ji,jj)*tmasklev(iilev,ijlev)*1.d0
+ d_n (iilev,ijlev) = d_n (iilev,ijlev) + tmask(ji,jj) *zbt(ji,jj)*tmasklev(iilev,ijlev)*1.d0
+ ENDIF
+ ENDDO
+ ENDDO
+
+ WHERE ( d_n > 0. )
+ d_out = d_out / d_n
+ ELSEWHERE
+ d_out = stypvarout(1)%rmissing_value
+ END WHERE
+
+ ! Check if there are points with missing values on Levitus grid
+ IF ( COUNT( d_out == stypvarout(1)%rmissing_value .AND. tmasklev == 1. ) /= 0. ) THEN
+ ALLOCATE ( z_fill(npilev,npjlev) )
+ z_fill(:,:) = 0.
+ !
+ imethod = 1
+ SELECT CASE (imethod)
+ CASE ( 1 ) ! Method 1: fill missing data with shapiro
+ ctyp='fill'
+ CALL shapiro_fill_smooth ( REAL(d_out), npilev, npjlev, 3, ctyp, stypvarout(1)%rmissing_value, INT(tmasklev), z_fill )
+
+ DO jjlev = 1 , npjlev
+ DO jilev = 1 , npilev
+ IF ( z_fill(jilev,jjlev) .NE. stypvarout(1)%rmissing_value &
+ & .AND. tmasklev(jilev,jjlev) == 1 &
+ & .AND. d_out(jilev,jjlev) == stypvarout(1)%rmissing_value ) &
+ & d_out(jilev,jjlev) = z_fill(jilev,jjlev)
+ ENDDO
+ ENDDO
+
+ CASE ( 2 ) ! Method 2: compute with influence bubble
+ ! For each point of Levitus grid, a data screening is performed
+ ! in a influence bubble of radius zradius, centered on Levitus point
+ ! and the weighted average of the data in the bubble is computed
+ z_fill(:,:) = 0.
+ DO jjlev = 1 , npjlev-1
+ DO jilev = 1 , npilev
+ ierr = 0
+ IF ( tmasklev(jilev,jjlev) == 1 .AND. d_out(jilev,jjlev) == stypvarout(1)%rmissing_value ) THEN
+ ! for the South pole, no treatment performed if data too far from southern most orca points
+ CALL btoe(rlonlev(jilev,jjlev),rlatlev(jilev,jjlev),rlon1,rlat1,-1.2*zradius,-1.2*zradius)
+ IF ( rlat1 > gphitmin ) THEN
+ ! Search the closest point of ORCA grid for this Levitus point
+ CALL cdf_findij (rlonlev(jilev,jjlev), rlonlev(jilev,jjlev), rlatlev(jilev,jjlev), rlatlev(jilev,jjlev), &
+ & iimin, iimax, ijmin, ijmax,cd_coord=cn_fhgr,cd_point='T' ) !, cd_verbose='N')
+
+ ! Next valid grid point going northward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin
+ DO WHILE ( ij <= npjglo .AND. ltest .AND. &
+ & etobd(rlonlev(jilev,jjlev),rlatlev(jilev,jjlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,jjlev) = z_fill(jilev,jjlev) + &
+ &(z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,jjlev)
+ d_n (jilev,jjlev) = d_n (jilev,jjlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,jjlev)
+ ierr = ierr + 1
+ ENDIF
+ ij = ij + 1
+ END DO
+ ! Next valid grid point going southward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin-1 ; ii = iimin
+ DO WHILE ( ij >= 1 .AND. ltest .AND. &
+ & etobd(rlonlev(jilev,jjlev),rlatlev(jilev,jjlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,jjlev) = z_fill(jilev,jjlev) + &
+ & (z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,jjlev)
+ d_n (jilev,jjlev) = d_n (jilev,jjlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,jjlev)
+ ierr = ierr + 1
+ ENDIF
+ ij = ij - 1
+ END DO
+ ! Next valid grid point going westward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin+1 ; IF ( ii > npiglo ) ii = ii - npiglo
+ DO WHILE ( ltest .AND. &
+ & etobd(rlonlev(jilev,jjlev),rlatlev(jilev,jjlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,jjlev) = z_fill(jilev,jjlev) + &
+ & (z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,jjlev)
+ d_n (jilev,jjlev) = d_n (jilev,jjlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,jjlev)
+ ierr = ierr + 1
+ ENDIF
+ ii = ii + 1
+ IF ( ii > npiglo ) ii = ii - npiglo
+ END DO
+ ! Next valid grid point going eastward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin-1 ; IF ( ii < 1) ii = ii + npiglo
+ DO WHILE ( ltest .AND. &
+ & etobd(rlonlev(jilev,jjlev),rlatlev(jilev,jjlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,jjlev) = z_fill(jilev,jjlev) + &
+ & (z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,jjlev)
+ d_n (jilev,jjlev) = d_n (jilev,jjlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,jjlev)
+ ierr = ierr + 1
+ ENDIF
+ ii = ii - 1
+ IF ( ii < 1 ) ii = ii + npiglo
+ END DO
+
+ ! computing d_out value
+ IF ( z_fill(jilev,jjlev) .NE. stypvarout(1)%rmissing_value &
+ & .AND. d_n(jilev,jjlev) > 0 &
+ & .AND. d_out(jilev,jjlev) == stypvarout(1)%rmissing_value ) &
+ & d_out(jilev,jjlev) = z_fill(jilev,jjlev) / d_n(jilev,jjlev)
+
+ ENDIF ! rlat1 > gphitmin
+
+ ENDIF ! tmasklev(jilev,jjlev) == 1 .AND. d_out(jilev,jjlev) == stypvarout(1)%rmissing_value
+
+ ENDDO
+ ENDDO
+
+ ! Case of the North Pole
+ ijlev = npjlev
+ DO jilev = 1 , npilev
+ ierr = 0
+
+ IF ( tmasklev(jilev,ijlev) == 1 .AND. d_out(jilev,ijlev) == stypvarout(1)%rmissing_value ) THEN
+
+ ! Search the closest point of ORCA grid for this Levitus point
+ CALL cdf_findij (rlonlev(jilev,ijlev), rlonlev(jilev,ijlev), rlatlev(jilev,ijlev), rlatlev(jilev,ijlev), &
+ & iimin, iimax, ijmin, ijmax,cd_coord=cn_fhgr,cd_point='T') !, cd_verbose='N')
+
+ ! Next valid grid point going southward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin
+ DO WHILE ( ij >= 1 .AND. ltest .AND. &
+ & etobd(rlonlev(jilev,ijlev),rlatlev(jilev,ijlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,ijlev) = z_fill(jilev,ijlev) + &
+ & (z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,ijlev)
+ d_n (jilev,ijlev) = d_n (jilev,ijlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,ijlev)
+ ierr = ierr + 1
+ ENDIF
+ ij = ij - 1
+ END DO
+ ! Next valid grid point going westward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin+1 ; IF ( ii > npiglo) ii = ii - npiglo
+ DO WHILE ( ltest .AND. &
+ & etobd(rlonlev(jilev,ijlev),rlatlev(jilev,ijlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,ijlev) = z_fill(jilev,ijlev) + &
+ & (z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,ijlev)
+ d_n (jilev,ijlev) = d_n (jilev,ijlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,ijlev)
+ ierr = ierr + 1
+ ENDIF
+ ii = ii + 1
+ IF ( ii > npiglo ) ii = ii - npiglo
+ END DO
+ ! Next valid grid point going eastward on ORCA grid
+ ltest = .TRUE. ; ij = ijmin ; ii = iimin-1 ; IF ( ii < 1 ) ii = ii + npiglo
+ DO WHILE ( ltest .AND. &
+ & etobd(rlonlev(jilev,ijlev),rlatlev(jilev,ijlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,ijlev) = z_fill(jilev,ijlev) + &
+ & (z_in(ii,ij)*tmask(ji,ij))*zbt(ii,ij)*tmasklev(jilev,ijlev)
+ d_n (jilev,ijlev) = d_n (jilev,ijlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,ijlev)
+ ierr = ierr + 1
+ ENDIF
+ ii = ii - 1 ; IF ( ii < 1 ) ii = ii + npiglo
+ END DO
+ ! Going "northward" and crossing the pole
+ ltest = .TRUE.
+ rpos = 1.
+ ij = ijmin + 1 * rpos ; ii = iimin
+ IF ( ij > npjglo ) THEN
+ ii = ii + npiglo/2 ; IF ( ii > npiglo ) ii = ii - npiglo
+ rpos = -1.
+ ij = ij + 1 * rpos
+ ENDIF
+ DO WHILE ( ij >= 1 .AND. ltest .AND. &
+ & etobd(rlonlev(jilev,ijlev),rlatlev(jilev,ijlev),glamt(ii,ij), gphit(ii,ij)) <= zradius )
+ IF ( tmask(ii,ij) == 1 ) THEN
+ ltest = .FALSE.
+ z_fill(jilev,ijlev) = z_fill(jilev,ijlev) + &
+ &(z_in(ii,ij)*tmask(ii,ij))*zbt(ii,ij)*tmasklev(jilev,ijlev)
+ d_n (jilev,ijlev) = d_n (jilev,ijlev) + &
+ & tmask(ii,ij) *zbt(ii,ij)*tmasklev(jilev,ijlev)
+ ierr = ierr + 1
+ ENDIF
+ ij = ij + 1 * rpos
+ IF ( ij > npjglo ) THEN
+ ii = ii + npiglo/2 ; IF ( ii > npiglo ) ii = ii - npiglo
+ rpos = -1.
+ ij = ij + 1 * rpos
+ ENDIF
+ END DO
+
+ ! computing d_out value
+ IF ( z_fill(jilev,ijlev) .NE. stypvarout(1)%rmissing_value &
+ & .AND. d_n(jilev,ijlev) > 0 ) &
+ & d_out(jilev,ijlev) = z_fill(jilev,ijlev) / d_n(jilev,ijlev)
+
+ ENDIF
+
+ ENDDO
+
+
+ CASE DEFAULT
+ PRINT *, ' METHOD ', imethod ,'is not recognized in this program'
+ STOP
+
+ END SELECT ! imethod
+ IF ( ALLOCATED(z_fill) ) DEALLOCATE( z_fill )
+
+ ENDIF ! filling points
+ ! ----------------------------------------------------------------------------------------
+ ! write
+ ierr = putvar(ncout, id_varout(1), REAL(d_out(:,:)), 1, npilev, npjlev, ktime=jt)
+
+ END DO ! loop on time
+
+ ierr = closeout(ncout)
+
+CONTAINS
+
+REAL(KIND=4) FUNCTION etobd(plonr, platr, plon, plat)
+ !!---------------------------------------------------------------------
+ !! *** FUNCTION etobd ***
+ !!
+ !! ** Purpose :
+ !! Compute the beta plane distance between two lon/lat values
+ !!
+ !! ** Method :
+ !! Return the distance (km) between 2 points given in the
+ !! arguments with their latitudes and longitudes
+ !!----------------------------------------------------------------------
+ REAL(KIND=4), INTENT(in) :: plonr, platr, plon, plat ! lon/lat of the 2 input points
+
+ REAL(KIND=8), PARAMETER :: dp_p1 = 1.745329D-02 ! radians per degree
+ REAL(KIND=8), PARAMETER :: dp_p2 = 111.1940D0 ! dp_p1 * earth radius in km (6370.949)
+
+ REAL(KIND=8) :: dl_rx, dl_ry ! working double prec variables
+ REAL(KIND=8) :: dl_r0
+ REAL(KIND=8) :: dl_r1, dl_r2
+ !!----------------------------------------------------------------------
+
+ dl_r0 = plonr ; IF ( dl_r0 < 0. ) dl_r0 = dl_r0 + 360.
+ dl_r1 = plon ; IF ( dl_r1 < 0. ) dl_r1 = dl_r1 + 360.
+
+ dl_rx = dl_r1 - dl_r0
+ IF ( dl_rx > 180. ) THEN
+ dl_rx = dl_rx - 360.
+ ELSE IF ( dl_rx < -180. ) THEN
+ dl_rx = dl_rx + 360.
+ ELSE IF ( ABS(dl_rx) == 180. ) THEN
+ dl_rx = 180.
+ ENDIF
+ dl_r2 = plat - platr
+
+ dl_rx = dp_p2*dl_rx*COS(dp_p1*platr)
+ dl_ry = dp_p2*dl_r2
+
+ etobd = SQRT(dl_rx*dl_rx + dl_ry*dl_ry)
+
+END FUNCTION etobd
+
+SUBROUTINE btoe(plonr, platr, plon, plat, plx, ply)
+
+ REAL(KIND=4), INTENT(in) :: plonr,platr,plx,ply
+ REAL(KIND=4), INTENT(out) :: plon,plat
+
+ REAL(KIND=8),PARAMETER :: dp_p1= 1.745329D-02 ! radians per degree
+ REAL(KIND=8) :: dp_p2= 111.1940 ! dp_p1 * earth radius in km (6370.949)
+ REAL(KIND=8) :: dl_rx,dl_ry,dl_r0,dl_r1,dl_r2,dl_r3
+
+ dl_r0 = plonr ; IF ( dl_r0 < 0. ) dl_r0 = dl_r0 + 360.
+ dl_rx = plx
+ dl_r2 = platr
+ dl_r1 = dl_r0 + dl_rx / ( dp_p2*COS(dp_p1*dl_r2) )
+ dl_r3 = dl_r2 + ply / dp_p2
+ plon = dl_r1
+ IF ( plon < 0. ) plon = plon + 360.
+ IF ( plon >= 360. ) plon = plon - 360.
+ plat = dl_r3
+ IF ( plat > 90. ) plat = 90.
+ IF ( plat < -90. ) plat = -90.
+
+END SUBROUTINE btoe
+
+END PROGRAM cdf2levitusgrid2d
--
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