[cdftools] 175/228: JMM: propset and cosmetics in cdf2levitusgrid2d. propset and change some var to dble prec for cdfgeostrophy.f90
Alastair McKinstry
mckinstry at moszumanska.debian.org
Fri Jun 12 08:21:46 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 00d612eb33335f19e14f0af294be7b08f46c2973
Author: molines <molines at 1055176f-818a-41d9-83e1-73fbe5b947c5>
Date: Wed Oct 31 20:37:17 2012 +0000
JMM: propset and cosmetics in cdf2levitusgrid2d.
propset and change some var to dble prec for cdfgeostrophy.f90
git-svn-id: http://servforge.legi.grenoble-inp.fr/svn/CDFTOOLS/trunk@623 1055176f-818a-41d9-83e1-73fbe5b947c5
---
cdf2levitusgrid2d.f90 | 1015 +++++++++++++++++++++++++------------------------
cdfgeostrophy.f90 | 23 +-
2 files changed, 526 insertions(+), 512 deletions(-)
diff --git a/cdf2levitusgrid2d.f90 b/cdf2levitusgrid2d.f90
index e77b69e..10632a9 100644
--- a/cdf2levitusgrid2d.f90
+++ b/cdf2levitusgrid2d.f90
@@ -1,507 +1,520 @@
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)
+ !!======================================================================
+ !! *** 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 2012
+ !! $Id$
+ !! Copyright (c) 2012, 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 ! IJ coordinates of the closest points
+ INTEGER(KIND=4) :: ijmin, ijmax ! " " "
+ INTEGER(KIND=4) :: imethod=1 ! interpolation method
+ INTEGER(KIND=4) :: iter_shap=3 ! number of Shapiro iteration
+ 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) :: 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'
+ iter_shap = 3 ! number of shapiro iteration
+ CALL shapiro_fill_smooth ( REAL(d_out), npilev, npjlev, iter_shap, 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
+ 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 = DBLE(plonr) ; IF ( dl_r0 < 0.d0 ) dl_r0 = dl_r0 + 360.d0
+ dl_r1 = DBLE(plon ) ; IF ( dl_r1 < 0.d0 ) dl_r1 = dl_r1 + 360.d0
+
+ dl_rx = dl_r1 - dl_r0
+ IF ( dl_rx > 180. ) THEN
+ dl_rx = dl_rx - 360.d0
+ ELSE IF ( dl_rx < -180.d0 ) THEN
+ dl_rx = dl_rx + 360.d0
+ ELSE IF ( ABS(dl_rx) == 180.d0 ) THEN
+ dl_rx = 180.d0
+ ENDIF
+ dl_r2 = DBLE(plat) - DBLE(platr)
+
+ dl_rx = dp_p2 * dl_rx * COS(dp_p1*platr)
+ dl_ry = dp_p2 * dl_r2
+
+ etobd = REAL(SQRT(dl_rx*dl_rx + dl_ry*dl_ry))
+
+ END FUNCTION etobd
+
+
+ SUBROUTINE btoe(plonr, platr, plon, plat, plx, ply)
+ !!---------------------------------------------------------------------
+ !! *** ROUTINE btoe ***
+ !!
+ !! ** Purpose : Return position (lon/lat) of a point located at
+ !! a given distance from the original point
+ !!
+ !! ** Method : Trigo ...
+ !!
+ !!----------------------------------------------------------------------
+ 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), PARAMETER :: dp_p2= 111.1940 ! dp_p1 * earth radius in km (6370.949)
+ REAL(KIND=8) :: dl_rx, dl_ry
+ REAL(KIND=8) :: dl_r0
+ REAL(KIND=8) :: dl_r1, dl_r2, dl_r3
+ !!----------------------------------------------------------------------
+ dl_r0 = DBLE(plonr) ; IF ( dl_r0 < 0.d0 ) dl_r0 = dl_r0 + 360.d0
+ dl_rx = DBLE(plx )
+ dl_r2 = DBLE(platr)
+ dl_r1 = dl_r0 + dl_rx / ( dp_p2 * COS(dp_p1*dl_r2) )
+ dl_r3 = dl_r2 + ply / dp_p2
+ plon = REAL(dl_r1)
+ IF ( plon < 0. ) plon = plon + 360.
+ IF ( plon >= 360. ) plon = plon - 360.
+ plat = REAL(dl_r3)
+ IF ( plat > 90. ) plat = 90.
+ IF ( plat < -90. ) plat = -90.
+
+ END SUBROUTINE btoe
END PROGRAM cdf2levitusgrid2d
diff --git a/cdfgeostrophy.f90 b/cdfgeostrophy.f90
index 9e73825..e255e24 100644
--- a/cdfgeostrophy.f90
+++ b/cdfgeostrophy.f90
@@ -68,7 +68,7 @@ PROGRAM cdfgeostrophy
USE modcdfnames
!!----------------------------------------------------------------------
!! CDFTOOLS_3.0 , MEOM 2011
- !! $Id: cdfgeo-uv.f90 539 2011-07-11 10:33:35Z molines $
+ !! $Id$
!! Copyright (c) 2011, J.-M. Molines
!! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
!!----------------------------------------------------------------------
@@ -94,18 +94,19 @@ PROGRAM cdfgeostrophy
REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: e3 ! vertic metrics
REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamu, gphiu ! longitude latitude u-point
REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: glamv, gphiv ! longitude latitude v-point
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: un, vn ! velocity components
- REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zsshn ! ssh
REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: umask, vmask ! mask at u and v points
REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmask ! mask at t points
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zsigsurf ! density at first level (used for zpsurf)
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zsiglevel ! density at current level (used for zplevel/zphalflevel)
+ REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zt, zsal ! temporary arrays for temperature and salinity
+
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: un, vn ! velocity components
+ REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zsshn ! ssh
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zpupper ! total pressure above current level
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zphalflevel ! pressure at T-point of current level
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zplevel ! pressure at bottom W-point of current level
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zpsurf ! pressure due to SSH
REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: zptot ! total pressure at current level
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zsigsurf ! density at first level (used for zpsurf)
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zsiglevel ! density at current level (used for zplevel/zphalflevel)
- REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE :: zt, zsal ! temporary arrays for temperature and salinity
REAL(KIND=4) :: zhlevel ! thickness of current level
REAL(KIND=4) :: zhhalflevel ! thickness of half the current level
@@ -252,7 +253,7 @@ PROGRAM cdfgeostrophy
! Compute psurf (pressure due to SSH)
zpsurf(:,:) = zsigsurf * grav * zsshn
- zpupper(:,:) = 0.
+ zpupper(:,:) = 0.d0
DO jk=1,npk
@@ -286,8 +287,8 @@ PROGRAM cdfgeostrophy
!! 2. We compute the velocities from geostrophic balance
- un(:,:) = 0.
- vn(:,:) = 0.
+ un(:,:) = 0.d0
+ vn(:,:) = 0.d0
DO jj=2,npjglo-1
DO ji=2,npiglo-1
@@ -314,8 +315,8 @@ PROGRAM cdfgeostrophy
ENDDO
ENDDO
- WHERE ( ABS(ff) < 1.e-5 ) un(:,:) = 0.
- WHERE ( ABS(ff) < 1.e-5 ) vn(:,:) = 0.
+ WHERE ( ABS(ff) < 1.e-5 ) un(:,:) = 0.d0
+ WHERE ( ABS(ff) < 1.e-5 ) vn(:,:) = 0.d0
un(:,:) = un(:,:) * umask(:,:)
vn(:,:) = vn(:,:) * vmask(:,:)
--
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