source: trunk/procs/meshes/mesh_regular.pro @ 45

PRO mesh_regular, deltax, deltay, INILON = inilon,  EQUATOR = equator,  NOMASK = nomask, MASK_FILE = mask_file, I_INDEX = i_index, DELTA_I = delta_i, J_INDEX = j_index, DELTA_J = delta_j, NO_SHIFT = no_shift, WHOLE_ARRAYS = whole_arrays, REVERSE_Y = reverse_y, GLAMBOUNDARY = glamboundary
@common
@com_eg
;
; init grid, sf, masks for regular grid (includes equator)
;
   print,' Regular grid inits.',deltax, deltay

   jpi = floor(360./deltax)
   jpj = floor(180./deltay)
   IF keyword_set(EQUATOR) THEN jpj = jpj + 1
   jpk = 1

 ; initialisation of character variables used in the execution of computegrid
   shift_txt =  ''
   periodic_txt =  ''
   idx_txt = ''

; 1.  Define longitudes 

   IF NOT keyword_set(inilon) THEN  inilon = 0.

   glamt = 360.0*findgen(jpi)/float(jpi)+inilon
   glamt = glamt#replicate(1, jpj)

; 2.  Define latitudes

   IF keyword_set(EQUATOR) THEN BEGIN 

      gphit = 90.0*findgen((jpj-1)/2)/float((jpj-1)/2)+deltay
      gphitn = reverse(-gphit)
      gphit = [gphitn, 0.0, gphit]

   END ELSE BEGIN

      gphit = 90.0*findgen(jpj/2)/float(jpj/2)+deltay/2.
      gphitn = reverse(gphit)
      gphit = [gphitn, -gphit]

   ENDELSE 

; 3. define mask

   tmask = lonarr(jpi, jpj)
   IF NOT keyword_set(NOMASK) THEN BEGIN  
      s_file = hom_idl+'grids/mask_'+mask_file+'.dat'
      print, '     Read mask from ',s_file 
      restore, s_file
;      nummsk = 12
;      s_file = hom_idl+'grids/mask_regular_'+ $
;       strtrim(string(jpi), 2)+'x'+strtrim(string(jpj), 2)
;      openr, nummsk, s_file, /get_lun, /f77
;      readu, nummsk, tmask
;   tmask = reverse(1-tmask, 2)
;      close, nummsk
;      free_lun, nummsk

   ENDIF ELSE BEGIN 
      tmask[*, *]= 1
      print, '     Warning, No mask read'
   ENDELSE 

; 4. Reduce grid in longitude

   IF keyword_set(i_index) THEN BEGIN 
      idx_txt = ',XMINMESH='+string(i_index)+',XMAXMESH='+string(i_index+delta_i-1)
      idx_txt = idx_txt+',XMINDTA='+string(i_index)+',XMAXDTA='+string(i_index+delta_i-1)
   ENDIF 
   IF keyword_set(j_index) THEN BEGIN 
      idx_txt = idx_txt+',YMINMESH='+string(j_index)+',YMAXMESH='+string(j_index+delta_j-1)
      idx_txt = idx_txt+',YMINDTA='+string(j_index)+',YMAXDTA='+string(j_index+delta_j-1)
   ENDIF 

   IF keyword_set(REVERSE_Y) THEN gphit = reverse(gphit)
   gphit = replicate(1, jpi)#gphit

   masked_data = 0
   mesh_type =  'atm'
;
; definition of key_shift
;
   IF keyword_set(NO_SHIFT) THEN BEGIN
      shift_txt = ', SHIFT = 0 '
   ENDIF
   IF keyword_set(WHOLE_ARRAYS) THEN BEGIN
      shift_txt = ', SHIFT = 0 '
      periodic_txt =  ', PERIODIC = 0'
   ENDIF

; Use the computegrid routine
   cmd_grid =  'computegrid, XAXIS = glamt, YAXIS = gphit, MASK = tmask, GLAMBOUNDARY = glamboundary, /FULLCGRID'+shift_txt+periodic_txt+idx_txt
   IF debug_w THEN print, cmd_grid
   res =  execute( cmd_grid )

   print, '   key_shift =', key_shift


   key_offset = [0, 0, 0]
;
; indice i pour grille j moyenne zonale
;
   diaznl_idx = 1

  return
END