source: trunk/SRC/Grid/smallmeshmask.pro

;+
;
; @categories
; For OPA
;
; @param inid {in}{required}
;
; @param outid {in}{required}
;
; @param inname {in}{required}
;
; @param outname {in}{optional}
;
; @history
;      July 2004 Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO ncdf_transfer, inid, outid, inname, outname
;
  compile_opt idl2, strictarrsubs
;
  IF n_elements(outname) EQ 0 THEN outname = inname
  ncdf_varget, inid, inname, zzz
  ncdf_varput, outid, outname, float(reform(zzz, /over))
  RETURN
END
;
;+
;
; @file_comments
; Reduce the size of the NetCDF meshmask created by OPA by
; using byte format for the masks and the float format for the other
; fields.
;
; @keyword IODIR
; to define the files path.
;
; @param ncfilein {in}{required}
; 1) the name of the meshmask file to be reduced. In that case,
; there is only one meshmask file
;
; OR
;
; 2) the xxx part in the names: xxx.mesh_hgr.nc xxx.mesh_zgr.nc
; xxx.mask.nc. In that case, the meshmask is split into 3 files.
;
; @param ncfileout {in}{optional}{default=smallmeshmask.nc}
; the name of the reduced meshmask file.
;
; @examples
;
;   IDL> meshdir='/d1fes2-raid2/smasson/DATA/ORCA05/'
;   IDL> smallmeshmask, 'meshmask_ORCA_R05.nc',iodir=meshdir
;
; @categories
; For OPA
;
; @history
;      July 2004 Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO smallmeshmask, ncfilein, ncfileout, IODIR = iodir, HGRFILE = hgrfile, ZGRFILE = zgrfile, MSKFILE = mskfile, USEBATHY = usebathy
;
  compile_opt idl2, strictarrsubs
;
  filein = isafile(FILE = ncfilein, IODIR = iodir, /NEW)
  test = (file_search(filein))[0]
  IF test EQ '' THEN BEGIN
    IF keyword_set(hgrfile) THEN filein_hgr = isafile(FILE = hgrfile, IODIR = iodir) ELSE filein_hgr = (file_search(filein+'[._]mesh_hgr.nc'))[0]
    IF keyword_set(zgrfile) THEN filein_zgr = isafile(FILE = zgrfile, IODIR = iodir) ELSE filein_zgr = (file_search(filein+'[._]mesh_zgr.nc'))[0]
    IF keyword_set(usebathy) THEN filein_msk = 'not used' ELSE BEGIN
      IF keyword_set(mskfile) THEN filein_msk = isafile(FILE = mskfile, IODIR = iodir) ELSE filein_msk = (file_search(filein+'[._]mask.nc'))[0]
    ENDELSE
    IF filein_hgr EQ '' OR filein_zgr EQ '' OR filein_msk EQ ''  THEN BEGIN
      ras = report(['meshmask file(s) not found...', $
      filein+' does not exist', $
      filein+'.mesh_hgr.nc does not exist', $
      filein+'.mesh_zgr.nc does not exist', $
      filein+'.mask.nc does not exist'])
      return
    ENDIF
  ENDIF ELSE filein = test
;------------------------------------------------------
;------------------------------------------------------
; get the horizontal dimensions
  IF n_elements(filein_hgr) NE 0  THEN cdfid = ncdf_open(filein_hgr) $
  ELSE cdfid = ncdf_open(filein)
  ncdf_diminq, cdfid, 'x', name, jpi
  ncdf_diminq, cdfid, 'y', name, jpj
; get the vertical dimensions
  IF n_elements(filein_zgr) NE 0  THEN BEGIN
    ncdf_close, cdfid
    cdfid = ncdf_open(filein_zgr)
  ENDIF
  listdims = strlowcase(ncdf_listdims(cdfid))
  IF (where(listdims EQ 'z'))[0] NE -1 THEN ncdf_diminq, cdfid, 'z', name, jpk ELSE BEGIN
    dimid = (where(strmid(listdims, 0, 5) EQ 'depth'))[0]
    IF dimid NE -1 THEN ncdf_diminq, cdfid, dimid, name, jpk ELSE BEGIN
      dummy = report('We could not find the vertical dimension..., its name must be z or start with depth')
      return
    ENDELSE
  ENDELSE
; get the variables list related to the partial steps
  zgr_varlist = ncdf_listvars(cdfid)
  zgr_varlist = strtrim(strlowcase(zgr_varlist), 2)
; is e3t a real 3D array?
  IF (where(zgr_varlist EQ 'e3t'))[0] NE -1 THEN BEGIN
    varinq = ncdf_varinq(cdfid, 'e3t')
    if varinq.ndims GE 3 THEN BEGIN
      ncdf_diminq, cdfid, varinq.dim[0], name, iii
      ncdf_diminq, cdfid, varinq.dim[1], name, jjj
      ncdf_diminq, cdfid, varinq.dim[2], name, kkk
      IF iii EQ jpi AND jjj EQ jpj AND kkk EQ jpk THEN key_e3_3d = 1
    ENDIF
  ENDIF
; is gdept a real 3D array?
  IF (where(zgr_varlist EQ 'gdept'))[0] NE -1 THEN BEGIN
    varinq = ncdf_varinq(cdfid, 'gdept')
    if varinq.ndims GE 3 THEN BEGIN
      ncdf_diminq, cdfid, varinq.dim[0], name, iii
      ncdf_diminq, cdfid, varinq.dim[1], name, jjj
      ncdf_diminq, cdfid, varinq.dim[2], name, kkk
      IF iii EQ jpi AND jjj EQ jpj AND kkk EQ jpk THEN key_gdep_3d = 1
    ENDIF
  ENDIF
;------------------------------------------------------
;
;------------------------------------------------------
; define the output file
;------------------------------------------------------
  IF n_elements(ncfileout) EQ 0  THEN ncfileout = 'smallmeshmask.nc'

  cdfidout = ncdf_create(isafile(FILE = ncfileout, IODIR = iodir, /NEW), /clobber)
  ncdf_control, cdfidout, /nofill
; dimension
  dimidx = ncdf_dimdef(cdfidout, 'x', jpi)
  dimidy = ncdf_dimdef(cdfidout, 'y', jpj)
  dimidz = ncdf_dimdef(cdfidout, 'z', jpk)
; global attributes
  ncdf_attput, cdfidout, 'IDL_Program_Name', 'smallmeshmask.pro', /GLOBAL
  ncdf_attput, cdfidout, 'Creation_Date', systime(), /GLOBAL
; declaration des variables
  varid = lonarr(20)
; horizontal variables
  hgrlist = ['glamt', 'glamu', 'glamv', 'glamf' $
             , 'gphit', 'gphiu', 'gphiv', 'gphif' $
             , 'e1t', 'e1u', 'e1v', 'e1f' $
             , 'e2t', 'e2u', 'e2v', 'e2f']
  FOR h = 0, n_elements(hgrlist)-1 DO $
    varid[h] = ncdf_vardef(cdfidout, hgrlist[h], [dimidx, dimidy], /float)
; vertical variables
  zgrlist = ['e3t', 'e3w', 'gdept', 'gdepw']
  FOR z = 0, n_elements(zgrlist)-1 DO $
    varid[16+z] = ncdf_vardef(cdfidout, zgrlist[z], [dimidz], /float)
; variables related to the partial steps
  with_partial = keyword_set(key_e3_3d) OR (where(zgr_varlist EQ 'hdept'))[0] NE -1 $
                 OR (where(zgr_varlist EQ 'e3tp'))[0] NE -1 $
                 OR (where(zgr_varlist EQ 'e3t_ps'))[0] NE -1
  IF with_partial THEN BEGIN
    varid = [varid, ncdf_vardef(cdfidout, 'hdept', [dimidx, dimidy], /float)]
    varid = [varid, ncdf_vardef(cdfidout, 'hdepw', [dimidx, dimidy], /float)]
    varid = [varid, ncdf_vardef(cdfidout, 'e3t_ps', [dimidx, dimidy], /float)]
    varid = [varid, ncdf_vardef(cdfidout, 'e3w_ps', [dimidx, dimidy], /float)]
  ENDIF
;
  IF (where(zgr_varlist EQ 'mbathy'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'mbathy', [dimidx, dimidy], /short)]
; mask variable
  IF keyword_set(usebathy) THEN msklist = ['tmask'] ELSE msklist = ['tmask', 'umask', 'vmask', 'fmask']
  FOR m = 0, n_elements(msklist)-1 DO $
    varid = [varid, ncdf_vardef(cdfidout, msklist[m], [dimidx, dimidy, dimidz], /byte)]
;------------------------------------------------------
;------------------------------------------------------
  ncdf_control, cdfidout, /endef
;------------------------------------------------------
;
; get the horizontal variables
;
  IF n_elements(filein_hgr) NE 0  THEN BEGIN
    ncdf_close, cdfid
    cdfid = ncdf_open(filein_hgr)
  ENDIF
;
  FOR h = 0, n_elements(hgrlist)-1 DO $
    ncdf_transfer, cdfid, cdfidout, hgrlist[h]
;
; get the vertical variables
;
  IF n_elements(filein_zgr) NE 0  THEN BEGIN
    ncdf_close, cdfid
    cdfid = ncdf_open(filein_zgr)
  ENDIF
;
  inzgrlist =  zgrlist
  IF (where(zgr_varlist EQ 'gdept_0'))[0] NE -1 THEN inzgrlist =  inzgrlist+'_0'
  FOR z = 0, n_elements(zgrlist)-1 DO ncdf_transfer, cdfid, cdfidout, inzgrlist[z], zgrlist[z]
; partial step variables
  IF (where(zgr_varlist EQ 'hdept' ))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'hdept'
  IF (where(zgr_varlist EQ 'hdepw' ))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'hdepw'
  IF (where(zgr_varlist EQ 'e3tp'  ))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'e3tp', 'e3t_ps'
  IF (where(zgr_varlist EQ 'e3wp'  ))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'e3wp', 'e3w_ps'
  IF (where(zgr_varlist EQ 'e3t_ps'))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'e3t_ps'
  IF (where(zgr_varlist EQ 'e3w_ps'))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'e3w_ps'
  IF (where(zgr_varlist EQ 'mbathy'))[0] NE -1 THEN ncdf_transfer, cdfid, cdfidout, 'mbathy'
;
  bat = -1
  IF keyword_set(key_e3_3d)  THEN BEGIN
; get the bottom
    CASE (ncdf_varinq(cdfid, 'mbathy')).ndims OF
      2:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj]
      3:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1]
      4:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1, 1]
    ENDCASE
    bottom = 0L > ( long(bat) - 1L )
    bottom = lindgen(jpi, jpj) + jpi*jpj*temporary(bottom)
;
    CASE (ncdf_varinq(cdfid, 'e3t')).ndims OF
      3:ncdf_varget, cdfid, 'e3t', res, count = [jpi, jpj, jpk]
      4:ncdf_varget, cdfid, 'e3t', res, count = [jpi, jpj, jpk, 1]
    ENDCASE
    ncdf_varput, cdfidout, 'e3t_ps', (temporary(res))[bottom] * (bat NE 0L)
;
    CASE (ncdf_varinq(cdfid, 'e3w')).ndims OF
      3:ncdf_varget, cdfid, 'e3w', res, count = [jpi, jpj, jpk]
      4:ncdf_varget, cdfid, 'e3w', res, count = [jpi, jpj, jpk, 1]
    ENDCASE
    ncdf_varput, cdfidout, 'e3w_ps', (temporary(res))[bottom] * (bat NE 0L)
  ENDIF
;
  IF keyword_set(key_gdep_3d)  THEN BEGIN
; get the bottom
    IF bat[0] EQ -1 THEN BEGIN
      CASE (ncdf_varinq(cdfid, 'mbathy')).ndims OF
        2:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj]
        3:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1]
        4:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1, 1]
      ENDCASE
    ENDIF
    bottom = 0L > ( long(bat) - 1L )
    bottom = lindgen(jpi, jpj) + jpi*jpj*temporary(bottom)
;
    CASE (ncdf_varinq(cdfid, 'gdept')).ndims OF
      3:ncdf_varget, cdfid, 'gdept', res, count = [jpi, jpj, jpk]
      4:ncdf_varget, cdfid, 'gdept', res, count = [jpi, jpj, jpk, 1]
    ENDCASE
    ncdf_varput, cdfidout, 'hdept', (temporary(res))[bottom] * (bat NE 0L)
;
    bottom = jpi*jpj + temporary(bottom)
    CASE (ncdf_varinq(cdfid, 'gdepw')).ndims OF
      3:ncdf_varget, cdfid, 'gdepw', res, count = [jpi, jpj, jpk]
      4:ncdf_varget, cdfid, 'gdepw', res, count = [jpi, jpj, jpk, 1]
    ENDCASE
    ncdf_varput, cdfidout, 'hdepw', (temporary(res))[bottom] * (bat NE 0L)
  ENDIF
;
; mask
;
  IF keyword_set(usebathy) THEN BEGIN
    IF bat[0] EQ -1 THEN BEGIN
      CASE (ncdf_varinq(cdfid, 'mbathy')).ndims OF
        2:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj]
        3:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1]
        4:ncdf_varget, cdfid, 'mbathy', bat, count = [jpi, jpj, 1, 1]
      ENDCASE
    ENDIF
    FOR k = 0, jpk-1 DO BEGIN
      zz = bat GT k
      ncdf_varput, cdfidout, 'tmask', temporary(zz), offset = [0, 0, k], count = [jpi, jpj, 1]
    ENDFOR  
  ENDIF ELSE BEGIN
    IF n_elements(filein_msk) NE 0  THEN BEGIN
      ncdf_close, cdfid
      cdfid = ncdf_open(filein_msk)
    ENDIF
; loop on the vertical levels to limit the memory use
    FOR k = 0, jpk-1 DO BEGIN
      FOR m = 0, n_elements(msklist)-1 DO BEGIN
        CASE (ncdf_varinq(cdfid, msklist[m])).ndims OF
          3:ncdf_varget, cdfid, msklist[m], zzz, offset = [0, 0, k   ], count = [jpi, jpj, 1]
          4:ncdf_varget, cdfid, msklist[m], zzz, offset = [0, 0, k, 0], count = [jpi, jpj, 1, 1]
        ENDCASE
        ncdf_varput, cdfidout, msklist[m], byte(temporary(zzz)), offset = [0, 0, k], count = [jpi, jpj, 1]
      ENDFOR
    ENDFOR
  ENDELSE
;------------------------------------------------------
;------------------------------------------------------
  ncdf_close, cdfid
  ncdf_close, cdfidout

  RETURN
END