source: trunk/SRC/Grid/micromeshmask.pro @ 297

;+
;
; @file_comments
; Reduce the size of the NetCDF meshmask created by OPA by
; using bit (and not byte) format for the masks and the float format
; for the other fields.
;
; @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
;
;+
;
; @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='micromeshmask.nc'}
; the name of the uniq reduced meshmask file.
;
; @keyword IODIR
; to define the files path.
;
; @examples
; IDL> meshdir='/d1fes2-raid2/smasson/DATA/ORCA05/'
; IDL> micromeshmask, 'meshmask_ORCA_R05.nc',iodir=meshdir
;
;-
;
PRO micromeshmask, ncfilein, ncfileout, IODIR = iodir
;
  compile_opt idl2, strictarrsubs
;
  filein = isafile(FILE = ncfilein, IODIR = iodir, /NEW)
  test = (findfile(filein))[0]
  IF test EQ '' THEN BEGIN
    filein_hgr = (findfile(filein+'.mesh_hgr.nc'))[0]
    filein_zgr = (findfile(filein+'.mesh_zgr.nc'))[0]
    filein_msk = (findfile(filein+'.mask.nc'))[0]
    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
; for the mask, we use "its byte" representation -> its y dimension
; will be extended to be a multiple of 8, then it will be divided by
; 8. -> if (jpj mod 8) eq 0 the jpj_m=jpi/8 else jpj_m=jpi/8 + 1
  jpj_m = (jpj+7)/8
; 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
  varlist_ps = ncdf_listvars(cdfid)
  varlist_ps = strtrim(strlowcase(varlist_ps), 2)
;------------------------------------------------------
;------------------------------------------------------
;
;------------------------------------------------------
; define the output file
;------------------------------------------------------
  IF n_elements(ncfileout) EQ 0  THEN ncfileout = 'micromeshmask.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)
  dimidy_m = ncdf_dimdef(cdfidout, 'y_m', jpj_m)
  dimidz = ncdf_dimdef(cdfidout, 'z', jpk)
; global attributs
  ncdf_attput, cdfidout, 'IDL_Program_Name', 'micromeshmask.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
  IF (where(varlist_ps EQ 'hdept'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'hdept', [dimidx, dimidy], /float)]
  IF (where(varlist_ps EQ 'hdepw'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'hdepw', [dimidx, dimidy], /float)]
; old variable name. keep for compatibility with old run. Change e3tp to e3t_ps
  IF (where(varlist_ps EQ 'e3tp'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'e3t_ps', [dimidx, dimidy], /float)]
; old variable name. keep for compatibility with old run. Change e3wp to e3w_ps
  IF (where(varlist_ps EQ 'e3wp'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'e3w_ps', [dimidx, dimidy], /float)]
;
  IF (where(varlist_ps EQ 'e3t_ps'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'e3t_ps', [dimidx, dimidy], /float)]
  IF (where(varlist_ps EQ 'e3w_ps'))[0] NE -1 THEN $
    varid = [varid, ncdf_vardef(cdfidout, 'e3w_ps', [dimidx, dimidy], /float)]
;   IF (where(varlist_ps EQ 'e3u_ps'))[0] NE -1 THEN $
;     varid = [varid, ncdf_vardef(cdfidout, 'e3u_ps', [dimidx, dimidy], /float)]
;   IF (where(varlist_ps EQ 'e3v_ps'))[0] NE -1 THEN $
;     varid = [varid, ncdf_vardef(cdfidout, 'e3v_ps', [dimidx, dimidy], /float)]
; mask variable
  msklist = ['tmask', 'umask', 'vmask', 'fmask']
  FOR m = 0, n_elements(msklist)-1 DO BEGIN
    varid = [varid, ncdf_vardef(cdfidout, msklist[m] $
                                , [dimidx, dimidy_m, dimidz], /byte)]
;     ncdf_attput, cdfidout, varid[n_elements(varid)-1] $
;       , 'Comment', 'the mask is stored as bit. You must use ' $
;       +'the binary representation of the byte to get back the data.'
  ENDFOR
;------------------------------------------------------
;------------------------------------------------------
  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
;
  FOR z = 0, n_elements(zgrlist)-1 DO $
    ncdf_transfer, cdfid, cdfidout, zgrlist[z]
; partial step variables
  IF (where(varlist_ps EQ 'hdept'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'hdept'
  IF (where(varlist_ps EQ 'hdepw'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'hdepw'
  IF (where(varlist_ps EQ 'e3tp'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'e3tp', 'e3t_ps'
  IF (where(varlist_ps EQ 'e3wp'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'e3wp', 'e3w_ps'
  IF (where(varlist_ps EQ 'e3t_ps'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'e3t_ps'
  IF (where(varlist_ps EQ 'e3w_ps'))[0] NE -1 THEN $
    ncdf_transfer, cdfid, cdfidout, 'e3w_ps'
;   IF (where(varlist_ps EQ 'e3u_ps'))[0] NE -1 THEN $
;     ncdf_transfer, cdfid, cdfidout, 'e3u_ps'
;   IF (where(varlist_ps EQ 'e3v_ps'))[0] NE -1 THEN $
;     ncdf_transfer, cdfid, cdfidout, 'e3v_ps'
;
; mask
;
  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, 3 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
      zzz = byte(temporary(zzz))
; zzz must contain only 0 or 1
      zzz = temporary(zzz) MOD 2
; we transpose zzz because we need to work with the y dimension as the
; first dimension
      zzz = transpose(temporary(zzz))
; extend jpj to be a multiple of 8
      jpjadd = jpj_m*8-jpj
      IF jpjadd NE 0 THEN zzz = [temporary(zzz), bytarr(jpjadd, jpi)]
; reform zzz, to look like output of binary.pro
      zzz = reform(zzz, 8, 1, jpj_m, jpi, /over)
; convert into "its byte form"
      zzz = inverse_binary(temporary(zzz))
      ncdf_varput, cdfidout, msklist[m], transpose(temporary(zzz)) $
        , offset = [0, 0, k], count = [jpi, jpj_m, 1]
    ENDFOR
  ENDFOR
;------------------------------------------------------
;------------------------------------------------------
  ncdf_close, cdfid
  ncdf_close, cdfidout

  RETURN
END