source: trunk/procs/meshes/mesh_lmdz.pro @ 2

PRO mesh_lmdz, res, NO_SHIFT = no_shift, WHOLE_ARRAYS = whole_arrays
@common
@com_eg
;
; init grid, sf, masks for LMDZ grid
;
   CASE res OF
      'low': print, ' LMDZ low resolution inits '
      'high':print, ' LMDZ high resolution inits '
      ELSE: stop
   ENDCASE 

   CASE res OF
      'low': BEGIN & jpi = 96L & jpj = 72L & jpk = 19L & key_shift = 0 & END 
      'high': BEGIN & jpi = 144L & jpj = 97L & jpk = 19L & key_shift = 0 & END 
   ENDCASE 
   IF keyword_set(NO_SHIFT) THEN key_shift = 0

   jpiglo = jpi
   jpjglo = jpj
   jpkglo = jpk

   ixminmesh  =0
   ixmaxmesh  =jpi-1
;
   iyminmesh  =0
   iymaxmesh  =jpj-1
;
   izminmesh  =0
   izmaxmesh  =jpk-1

; 1.  Define longitudes 

   glamt = 360.0*findgen(jpi)/float(jpi)-180
   glamt = glamt#replicate(1, jpj)

; 2 Define latitudes
   CASE res OF
      'low': BEGIN
         gphit = 4*findgen(jpj/2)+2
         gphit = [reverse(-gphit), gphit]
         gphit = reverse(gphit)
      END 
      'high': BEGIN
         gphit = 2.5*findgen(jpj/2)+2.5
         gphit = [reverse(-gphit), 0, gphit]
         gphit = reverse(gphit)
      END 
   ENDCASE 
   gphit = replicate(1, jpi)#gphit

; 3 Define scale factors

   zrad=6371229.0

   ; Exact method: integration on a sphere

   e1t = abs(2*!pi*zrad*sin(gphit*!pi/180.0)/jpi)

   ytop = (shift(gphit, 0, -1)+gphit)*0.5*!pi/180.0
   ybot = (shift(gphit, 0, +1)+gphit)*0.5*!pi/180.0
   
   ytop(*, jpj-1) = !pi/2
   ybot(*,     0) = -!pi/2

   e2t = zrad*(ytop-ybot)

; 4. Define mask

; open mask file
   CASE res OF
      'low': BEGIN
;          model = 'lmdzl'
;          s_file = hom_idl+'grids/grids_'+model+'.nc'
;          tmask = nc_get(s_file, 'sftlf') 
;          tmask = shift(tmask,jpi/2, 0)
;          tmask = reverse(tmask, 2) 
;          tmask =  tmask < 1
;          tmask =  1-tmask
         nummsk = 12
         s_file = hom_idl+'grids/mask_lmdz_'+res
         openr, nummsk, s_file, /get_lun, /f77
         tmask = lonarr(jpi, jpj)
         readu, nummsk, tmask
; convention : 1 on sea, 0 on land
         tmask = 1-tmask
         close, nummsk
         free_lun, nummsk
      END 
      'high': BEGIN 
;         nummsk = 12
;         s_file = hom_idl+'grids/mask_lmdz_'+res
;         openr, nummsk, s_file, /get_lun, /f77
         tmask = lonarr(jpi, jpj)
;         readu, nummsk, tmask
; convention : 1 on sea, 0 on land
         tmask(*, *)= 1
;         close, nummsk
;         free_lun, nummsk
      END
   ENDCASE 

; 5. vertical grid (hPa)
   gdept = [ 100426.5, 98327.56, 95346.5, 90966.84, 84776.86, 76536.46, $ 
    66292.22, 54559.26, 42501.76, 31805.98, 23787.52, 18252.75, 13995.97, $ 
    10320.8, 7191.053, 4661.738, 2732.851, 1345.618, 388.2433 ]

   gdept = gdept/100.
   gdepw = gdept

   e3t = shift(gdept, 1)-gdept
   e3t[0] = e3t[1]

; Shift re-init U,V,F

; key_shit to have first longitude at 20E

   

   glamt = shift(glamt,key_shift, 0)
;   glamt(where(glamt LT 20) ) = glamt(where(glamt LT 20) )+360.
   gphit = shift(gphit,key_shift, 0)
   e1t = shift(e1t,key_shift, 0)
   e2t = shift(e2t,key_shift, 0)
   tmask = reform(shift(tmask,key_shift, 0), jpi*jpj)
   tmask = reform(tmask#replicate(1, jpk), jpi, jpj, jpk)

   glamu = glamt
   gphiu = gphit
   e1u = e1t
   e2u = e2t
   glamv = glamt
   gphiv = gphit
   e1v = e1t
   e2v = e2t

   glamf = 0.5*(shift(glamt, -1, 0)+glamt)
   gphif = 0.5*(shift(gphit, 0, -1)+gphit)
   glamf(jpi-1, *) =  glamf(jpi-2, *) + (glamf(jpi-2, *)-glamf(jpi-3, *))
   gphif(*, jpj-1) =  gphif(*, jpj-2) + (gphif(*, jpj-2)-gphif(*, jpj-3))
   e1f = e1t
   e2f = e2t
   
   e3w = e3t
   
   key_periodique=1

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

  return
END