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

;+
;
; @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
  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 = '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
  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
  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
;
  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'
;
  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
    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, '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 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