Subversion Repositories lagranto.ecmwf

Rev

Blame | Last modification | View Log | Download | RSS feed

; ecmwf_grib_to_netcdf_lagranto.ncl
; 
; NCL script to read in ECMWF grib data and output in netCDF format 
; for lagranto.
;
; Written by Andrew Penny - November 26, 2013
;=========================================================================

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

  ml_grb_DIR            = "/Volumes/Quadra4/model_data_4/ECMWF_YOTC_GLOBAL/model_level/"
  sfc_grb_DIR   = "/Volumes/Quadra4/model_data_4/ECMWF_YOTC_GLOBAL/surface/"
  out_DIR               = "/Volumes/Quadra4/model_data_4/ECMWF_YOTC_GLOBAL/netcdf/"
  basetime              = "hours since 1979-01-01 00:00" 
  write_const_file      = True
  
; START TIME
  s_yr  = 2008
  s_mo  = 8
  s_day = 25
  s_hr          = 0
  s_mn          = 0
  s_sec = 0

; END TIME  
  e_yr  = 2008
  e_mo  = 8
  e_day = 30
  e_hr          = 12
  e_mn          = 0
  e_sec = 0
  
; PARSE STARTING AND ENDING DATES
  start_date    = ut_inv_calendar(s_yr, s_mo, s_day, s_hr, s_mn, s_sec, basetime, 0)
  start_date_utc = ut_calendar(start_date,0) 
  start_date_out = start_date_utc(:,0) + "-" + where(start_date_utc(:,1).lt.10,"0"+start_date_utc(:,1),start_date_utc(:,1)) + "-" + where(start_date_utc(:,2).lt.10,"0"+start_date_utc(:,2),start_date_utc(:,2))+ "-" + where(start_date_utc(:,3).lt.10,"0"+start_date_utc(:,3),start_date_utc(:,3)) + where(start_date_utc(:,4).lt.10,"0"+start_date_utc(:,4),start_date_utc(:,4))

  end_date      = ut_inv_calendar(e_yr, e_mo, e_day, e_hr, e_mn, e_sec, basetime, 0)
  end_date_utc = ut_calendar(end_date,0) 
  end_date_out = end_date_utc(:,0) + "-" + where(end_date_utc(:,1).lt.10,"0"+end_date_utc(:,1),end_date_utc(:,1)) + "-" + where(end_date_utc(:,2).lt.10,"0"+end_date_utc(:,2),end_date_utc(:,2))+ "-" + where(end_date_utc(:,3).lt.10,"0"+end_date_utc(:,3),end_date_utc(:,3)) + where(end_date_utc(:,4).lt.10,"0"+end_date_utc(:,4),end_date_utc(:,4))      
  
  ml_FILES_ALL  = systemfunc ("csh -c 'ls "+ ml_grb_DIR + "*.grib'")    
  ml_FILEnames  = str_sub_str(systemfunc ("csh -c 'ls "+ ml_grb_DIR + "*.grib'"),ml_grb_DIR,"") 
    
  sfc_FILES_ALL         = systemfunc ("csh -c 'ls "+ sfc_grb_DIR + "*.grib'")    
  sfc_FILEnames         = str_sub_str(systemfunc ("csh -c 'ls "+ sfc_grb_DIR + "*.grib'"),sfc_grb_DIR,"")     
  
; PARSE FILE DATES
  ml_files_yr   = stringtoint(str_get_cols(ml_FILEnames, 0,3))
  ml_files_mo   = stringtoint(str_get_cols(ml_FILEnames, 4,5))
  ml_files_day          = stringtoint(str_get_cols(ml_FILEnames, 6,7))
  ml_files_hr   = stringtoint(str_get_cols(ml_FILEnames, 8,9))      
  ml_files_mn   = ispan(0,dimsizes(ml_files_yr)-1,1)*0
  ml_files_sec          = ispan(0,dimsizes(ml_files_yr)-1,1)*0
  ml_files_date         = ut_inv_calendar(ml_files_yr, ml_files_mo, ml_files_day, ml_files_hr, ml_files_mn, ml_files_sec, basetime, 0)
  ml_files_date_utc  = ut_calendar(ml_files_date, 0)
  ml_files_date_out     = ml_files_date_utc(:,0) + "-" + where( ml_files_date_utc(:,1).lt.10, "0" + ml_files_date_utc(:,1), ml_files_date_utc(:,1)) + "-" + where( ml_files_date_utc(:,2).lt.10, "0" + ml_files_date_utc(:,2), ml_files_date_utc(:,2)) + " " + where( ml_files_date_utc(:,3).lt.10, "0" + ml_files_date_utc(:,3), ml_files_date_utc(:,3)) + where( ml_files_date_utc(:,4).lt.10, "0" + ml_files_date_utc(:,4), ml_files_date_utc(:,4))         

  sfc_files_yr          = stringtoint(str_get_cols(sfc_FILEnames, 0,3))
  sfc_files_mo          = stringtoint(str_get_cols(sfc_FILEnames, 4,5))
  sfc_files_day         = stringtoint(str_get_cols(sfc_FILEnames, 6,7))
  sfc_files_hr          = stringtoint(str_get_cols(sfc_FILEnames, 8,9))      
  sfc_files_mn          = ispan(0,dimsizes(sfc_files_yr)-1,1)*0
  sfc_files_sec         = ispan(0,dimsizes(sfc_files_yr)-1,1)*0
  sfc_files_date        = ut_inv_calendar(sfc_files_yr, sfc_files_mo, sfc_files_day, sfc_files_hr, sfc_files_mn, sfc_files_sec, basetime, 0)
  sfc_files_date_utc    = ut_calendar(sfc_files_date, 0)
  sfc_files_date_out    = sfc_files_date_utc(:,0) + "-" + where( sfc_files_date_utc(:,1).lt.10, "0" + sfc_files_date_utc(:,1), sfc_files_date_utc(:,1)) + "-" + where( sfc_files_date_utc(:,2).lt.10, "0" + sfc_files_date_utc(:,2), sfc_files_date_utc(:,2)) + " " + where( sfc_files_date_utc(:,3).lt.10, "0" + sfc_files_date_utc(:,3), sfc_files_date_utc(:,3)) + where( sfc_files_date_utc(:,4).lt.10, "0" + sfc_files_date_utc(:,4), sfc_files_date_utc(:,4))         
  
; TAG FILES WITHINN START AND END DATES  
  ml_ind_files  = ind(ml_files_date.ge.start_date.and.ml_files_date.le.end_date)
  ml_FILES              = ml_FILES_ALL(ml_ind_files)
  ml_numFILES   = dimsizes(ml_FILES)    

  sfc_ind_files = ind(sfc_files_date.ge.start_date.and.sfc_files_date.le.end_date)
  sfc_FILES             = sfc_FILES_ALL(sfc_ind_files)
  sfc_numFILES  = dimsizes(sfc_FILES)   
  
  if(ml_numFILES.ne.sfc_numFILES)then
    print("model level / surface file number mismatch")
    exit
  end if
  
; DEFINE SOME CONSTANTS
  R                     = 287.05
  g                     = 9.80665  
  
; OPEN ECMWF YOTC COEFFICIENTS FILE AND READ IN VALUES
  strs = asciiread("/Users/abpenny/code/ncl/ecmwf_coeffs",-1,"string")

  delim = ","                                              
  N  = stringtoint(str_get_field(strs, 1, delim))       
  a  = stringtofloat(str_get_field(strs, 2, delim))     
  b  = stringtofloat(str_get_field(strs, 3, delim))       
  
  dims_N = dimsizes(N)  
  
; FLIP SO BOTTOM --> TOP
  aklev = a(dims_N-2:0)*0.01    ; need to multiply "a" coefficients by 0.01 so that they match surface pressure in hPa
  bklev = b(dims_N-2:0)

  aklay = new((/(dims_N-1)/),float,"No_FillValue") 
  bklay = new((/(dims_N-1)/),float,"No_FillValue")   

  do k = 0, dims_N - 3  
    aklay(k+1) = 0.5*(aklev(k)+aklev(k+1))
    bklay(k+1) = 0.5*(bklev(k)+bklev(k+1))
  end do
    
  aklay(0)=0.5*(0. + aklev(0))
  bklay(0)=0.5*(1. + bklev(0))    
  
  aklev!0               = "nz"    
  copy_VarCoords(aklev,bklev)        
  copy_VarCoords(aklev,aklay)        
  copy_VarCoords(aklev,bklay)            
      
;<><><><><><><><><><><><><><>
; --- LOOP THROUGH FILES ---

do ff = 0, ml_numFILES - 1

  print(" ")
  print("PROCESSING DATE: " + ml_files_date_out(ml_ind_files(ff)) ) 
     
  ml_grb_file   = addfile(ml_FILES(ff),"r")
  sfc_grb_file  = addfile(sfc_FILES(ff),"r")  
;  print(getfilevarnames(ml_grb_file)) 
;  print(getfilevarnames(sfc_grb_file)) 

; LATITUDE/LONGITUDE
  lat_ml                = (/ml_grb_file->g0_lat_0(::-1)/)                       ; latitude [degrees_north] (Cylindrical Equidistant Projection Grid), lat
  lon_subset            = (/ml_grb_file->g0_lon_1/)                             ; longitude [degrees_east] (Cylindrical Equidistant Projection Grid), lon
  dims_y                = dimsizes(lat_ml)
  dims_x                = dimsizes(lon_subset)+1
  dims_t                = 1
  dims_sfc              = 1
   
; MAKE LONGITUDE WRAP AROUND
  lon_ml                = new((/dims_x/),float,-999.99)   
  lon_ml(0:dims_x-2)    = lon_subset
  lon_ml(dims_x-1)      = 180.0  
  
; VERTICAL LEVELS
  lev_ml                = (/ml_grb_file->lv_HYBL2/)                             ; Hybrid level (atmosphere hybrid sigma pressure coordinate)
  dims_z                = dimsizes(lev_ml)
  
; 3D VARIABLES (flip latitude so monotonically increasing)
  t_ml                  = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   
  u_ml                  = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   
  v_ml                  = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   
  q_ml                  = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   
  w_ml                  = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   

  t_ml(0,:,:,0:dims_x-2)        = (/ml_grb_file->T_GDS0_HYBL(:,::-1,:)/)                ; temperature [K], lv_HYBL2/lat/lon
  u_ml(0,:,:,0:dims_x-2)        = (/ml_grb_file->U_GDS0_HYBL(:,::-1,:)/)                ; u wind [m s**-1], lv_HYBL2/lat/lon
  v_ml(0,:,:,0:dims_x-2)        = (/ml_grb_file->V_GDS0_HYBL(:,::-1,:)/)                ; v wind [m s**-1], lv_HYBL2/lat/lon
  q_ml(0,:,:,0:dims_x-2)        = (/ml_grb_file->Q_GDS0_HYBL(:,::-1,:)/)                ; specific humidity [kg kg**-1], lv_HYBL2/lat/lon
  w_ml(0,:,:,0:dims_x-2)        = (/ml_grb_file->W_GDS0_HYBL(:,::-1,:)/)                ; vertical velocity (omega) [Pa s**-1], lv_HYBL2/lat/lon

; WRAP AROUND
  t_ml(0,:,:,dims_x-1)  = (/ml_grb_file->T_GDS0_HYBL(:,::-1,0)/)                ; temperature [K], lv_HYBL2/lat/lon
  u_ml(0,:,:,dims_x-1)  = (/ml_grb_file->U_GDS0_HYBL(:,::-1,0)/)                ; u wind [m s**-1], lv_HYBL2/lat/lon
  v_ml(0,:,:,dims_x-1)  = (/ml_grb_file->V_GDS0_HYBL(:,::-1,0)/)                ; v wind [m s**-1], lv_HYBL2/lat/lon
  q_ml(0,:,:,dims_x-1)  = (/ml_grb_file->Q_GDS0_HYBL(:,::-1,0)/)                ; specific humidity [kg kg**-1], lv_HYBL2/lat/lon
  w_ml(0,:,:,dims_x-1)  = (/ml_grb_file->W_GDS0_HYBL(:,::-1,0)/)                ; vertical velocity (omega) [Pa s**-1], lv_HYBL2/lat/lon

; SURFACE VARIABLES
  sst_sfc               = new((/dims_t,dims_sfc,dims_y,dims_x/),float,-999.99) 
  z_sfc                 = new((/dims_t,dims_sfc,dims_y,dims_x/),float,-999.99)   
  ps_sfc                = new((/dims_t,dims_sfc,dims_y,dims_x/),float,-999.99)   
  mslp_sfc              = new((/dims_t,dims_sfc,dims_y,dims_x/),float,-999.99)   
  sst_sfc(0,0,:,0:dims_x-2)             = (/sfc_grb_file->SSTK_GDS0_SFC(::-1,:)/)               ; sea surface temperature [K], lat/lon
  z_sfc(0,0,:,0:dims_x-2)               = (/sfc_grb_file->Z_GDS0_SFC(::-1,:)/)          ; geopotential [m**2 s**-2], lat/lon  
  ps_sfc(0,0,:,0:dims_x-2)              = (/sfc_grb_file->SP_GDS0_SFC(::-1,:)/)         ; surface pressure [Pa], lat/lon - keep in Pa for GHT calculation below
  mslp_sfc(0,0,:,0:dims_x-2)    = (/sfc_grb_file->MSL_GDS0_SFC(::-1,:)/)                ; mean sea-level pressure [Pa], lat/lon    

; WRAP AROUND
  sst_sfc(0,0,:,dims_x-1)       = (/sfc_grb_file->SSTK_GDS0_SFC(::-1,0)/)               ; sea surface temperature [K], lat/lon
  z_sfc(0,0,:,dims_x-1)         = (/sfc_grb_file->Z_GDS0_SFC(::-1,0)/)          ; geopotential [m**2 s**-2], lat/lon  
  ps_sfc(0,0,:,dims_x-1)        = (/sfc_grb_file->SP_GDS0_SFC(::-1,0)/)         ; surface pressure [Pa], lat/lon - keep in Pa for GHT calculation below
  mslp_sfc(0,0,:,dims_x-1)      = (/sfc_grb_file->MSL_GDS0_SFC(::-1,0)/)                ; mean sea-level pressure [Pa], lat/lon    
  
;<><><><><><><><><><><><><><><>  
; COMPUTE GEOPOTENTIAL HEIGHT
;<><><><><><><><><><><><><><><>  
  z_ml          = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)   
  hgtprev       = z_sfc(0,0,:,:)

  do j = 1, dims_N-1     
  k = (dims_N-1) - j + 1  
  pend   = 0.5 * (a(k-1) + a(k)) + ps_sfc(0,0,:,:) * 0.5 * (b(k-1) + b(k))   
   
; IF USING vinth2p_ecmwf where we need to work with geopotential, not geopotential height
  if (k .eq. (dims_N-1) ) then               
    pstart = ps_sfc(0,0,:,:)
    z_ml(0,k-1,:,:) = hgtprev + R * t_ml(0,k-1,:,:) * (1. + 0.61 * q_ml(0,k-1,:,:)) * log(pstart/pend)     
  else
    pstart = 0.5 * (a(k) + a(k+1)) + ps_sfc(0,0,:,:) * 0.5 * (b(k) + b(k+1))     
    z_ml(0,k-1,:,:) = hgtprev + R * (t_ml(0,k,:,:) + t_ml(0,k-1,:,:))/2.*(1. + 0.61 * (q_ml(0,k,:,:)+q_ml(0,k-1,:,:))/2.) * log(pstart/pend) 
  end if      
    hgtprev = (/z_ml(0,k-1,:,:)/)    
  end do    
  
  z_ml = (/z_ml/g/)
  
; CONVERT SURFACE PRESSURE AND SLP TO hPA FROM Pa
  ps_sfc   = ps_sfc*0.01
  mslp_sfc = mslp_sfc*0.01
  
; COMPUTE PRESSURE
  prs_half      = new((/dims_N,dims_y,dims_x/),float,"No_FillValue") 
  p_ml          = new((/(dims_N-1),dims_y,dims_x/),float,"No_FillValue") 
  
  do k = 0, dims_N - 1
    prs_half(k,:,:) = a(k) + b(k)*ps_sfc(0,0,:,:)*100.0    
  end do
  
  do k = 0, dims_N - 2
    p_ml(k,:,:) = 0.5*(prs_half(k,:,:) + prs_half(k+1,:,:))     
  end do

; COMPUTE MIXING RATIO
  mixr_ml               = q_ml/(1.0-q_ml)
  
; COMPUTE RELATIVE HUMIDITY  
  rh_ml                 = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)     
  rh_ml(0,:,:,:)        = (/relhum(t_ml(0,:,:,:), mixr_ml(0,:,:,:), p_ml)/)
  rh_ml         = where(rh_ml.gt.100.0,100.0,rh_ml)

; COMPUTE THETA-E
  eth_ml                = new((/dims_t,dims_z,dims_y,dims_x/),float,-999.99)     
  eth_ml(0,:,:,:)       = (/wrf_eth(q_ml(0,:,:,:), t_ml(0,:,:,:), p_ml)/)
  
; FLIP SO BOTTOM --> TOP
  eth_ml(0,:,:,:)       = eth_ml(0,::-1,:,:)
  q_ml(0,:,:,:)         = q_ml(0,::-1,:,:)
  rh_ml(0,:,:,:)        = rh_ml(0,::-1,:,:)
  t_ml(0,:,:,:)         = t_ml(0,::-1,:,:)  
  z_ml(0,:,:,:)         = z_ml(0,::-1,:,:)    
  u_ml(0,:,:,:)         = u_ml(0,::-1,:,:)    
  v_ml(0,:,:,:)         = v_ml(0,::-1,:,:)    
  w_ml(0,:,:,:)         = w_ml(0,::-1,:,:)          
  
; PARSE TIME AGAIN JUST FOR KICKS
  time = ml_grb_file->T_GDS0_HYBL@initial_time  
  mo    = str_get_cols(time,0,1)
  day   = str_get_cols(time,3,4)  
  year = str_get_cols(time,6,9) 
  hr    = str_get_cols(time,12,13)
  mn    = str_get_cols(time,15,16)
  sec   = 0
  
  date  = ut_inv_calendar(stringtoint(year), stringtoint(mo), stringtoint(day), stringtoint(hr), stringtoint(mn), sec, basetime, 0)  
  date_utc      = ut_calendar(date,0)   
  
; SETUP COORDINATES FOR 2D VARIABLES 
  ps_sfc!0                      = "time"  
  ps_sfc&time                   = tofloat(date)
  ps_sfc&time@units             = basetime
  ps_sfc&time@calendar          = "proleptic_gregorian"
  ps_sfc&time@standard_name     = "time"

  ps_sfc!1                      = "dimz_PS"  

  ps_sfc!2                      = "dimy_Q"  
  ps_sfc&dimy_Q                         = lat_ml
  ps_sfc&dimy_Q@units           = "degrees_north"
  ps_sfc&dimy_Q@standard_name   = "latitude"
  ps_sfc&dimy_Q@_FillValue      = -999.99  

  ps_sfc!3                      = "dimx_Q"  
  ps_sfc&dimx_Q                         = lon_ml
  ps_sfc&dimx_Q@units           = "degrees_east"
  ps_sfc&dimx_Q@standard_name   = "longitude"
  ps_sfc&dimx_Q@_FillValue      = -999.99  

  ps_sfc@xmin                   = min(lon_ml)
  ps_sfc@xmax                   = max(lon_ml)  
  ps_sfc@xstag                  = 0.0

  ps_sfc@ymin                   = min(lat_ml)
  ps_sfc@ymax                   = max(lat_ml)  
  ps_sfc@ystag                  = 0.0

  ps_sfc@zmin                   = 1050.0
  ps_sfc@zmax                   = 0.0
  ps_sfc@zstag                  = -0.5

  ps_sfc@missing_data           = ps_sfc@_FillValue

  copy_VarMeta(ps_sfc,z_sfc)
  copy_VarMeta(ps_sfc,sst_sfc)  
  copy_VarMeta(ps_sfc,mslp_sfc)    
  
; SETUP COORDINATES FOR 3D VARIABLES
  t_ml!0                        = "time"  
  t_ml&time                     = tofloat(date)
  t_ml&time@units               = basetime
  t_ml&time@calendar    = "proleptic_gregorian"
  t_ml&time@standard_name       = "time"
  
  t_ml!1                        = "dimz_Q" 
  
  t_ml!2                        = "dimy_Q"  
  t_ml&dimy_Q           = lat_ml
  t_ml&dimy_Q@units     = "degrees_north"
  t_ml&dimy_Q@standard_name = "latitude"
  t_ml&dimy_Q@_FillValue        = -999.99  

  t_ml!3                        = "dimx_Q"  
  t_ml&dimx_Q           = lon_ml
  t_ml&dimx_Q@units     = "degrees_east"
  t_ml&dimx_Q@standard_name = "longitude"
  t_ml&dimx_Q@_FillValue        = -999.99  
  
  t_ml@xmin                     = min(lon_ml)
  t_ml@xmax                     = max(lon_ml)  
  t_ml@xstag                    = 0.0

  t_ml@ymin                     = min(lat_ml)
  t_ml@ymax                     = max(lat_ml)  
  t_ml@ystag                    = 0.0

  t_ml@zmin                     = 1050.0
  t_ml@zmax                     = 0.0
  t_ml@zstag                    = -0.5

  t_ml@missing_data             = t_ml@_FillValue

  copy_VarMeta(t_ml,q_ml)
  copy_VarMeta(t_ml,u_ml)  
  copy_VarMeta(t_ml,v_ml)    
  copy_VarMeta(t_ml,w_ml)  
  copy_VarMeta(t_ml,z_ml)    
  copy_VarMeta(t_ml,rh_ml)      
  copy_VarMeta(t_ml,eth_ml)      
          
; PRIMARY "P" FILE     
  print("WRITING " + "P"+year+mo+day+"_"+hr)
  outfile_P = addfile(out_DIR + "P"+year+mo+day+"_"+hr+".nc","c")  
  filedimdef(outfile_P, (/"time","dimz_PS","dimz_Q","dimy_Q","dimx_Q"/), (/dims_t,dims_sfc,dims_z,dims_y,dims_x/), (/True,False,False,False,False/)) 
        
; define global attributes:
                outfile_P@domxmin = min(lon_ml)
                outfile_P@domxmax = max(lon_ml)
                outfile_P@domymin = min(lat_ml)
                outfile_P@domymax = max(lat_ml)
                outfile_P@domzmin = 1050.0
                outfile_P@domzmax = 0.0
                outfile_P@constants_file_name = "yotc_cst"
  outfile_P->$("U")$            = u_ml      
  outfile_P->$("V")$            = v_ml        
  outfile_P->$("OMEGA")$        = w_ml          
  outfile_P->$("PS")$   = ps_sfc   
  
  system("mv -f " + out_DIR + "P"+year+mo+day+"_"+hr+".nc" + " " + out_DIR + "P"+year+mo+day+"_"+hr)                 

; SECONDARY "S" FILE     
  print("WRITING " + "S"+year+mo+day+"_"+hr)
  outfile_S = addfile(out_DIR + "S"+year+mo+day+"_"+hr+".nc","c")  
  filedimdef(outfile_S, (/"time","dimz_PS","dimz_Q","dimy_Q","dimx_Q"/), (/dims_t,dims_sfc,dims_z,dims_y,dims_x/), (/True,False,False,False,False/)) 
        
; define global attributes:
                outfile_S@domxmin = min(lon_ml)
                outfile_S@domxmax = max(lon_ml)
                outfile_S@domymin = min(lat_ml)
                outfile_S@domymax = max(lat_ml)
                outfile_S@domzmin = 1050.0
                outfile_S@domzmax = 0.0
                outfile_S@constants_file_name = "yotc_cst"                 
  outfile_S->$("T")$            = t_ml
  outfile_S->$("Q")$            = q_ml
  outfile_S->$("Z")$            = z_ml
  outfile_S->$("RH")$   = rh_ml
  outfile_S->$("THE")$  = eth_ml
  outfile_S->$("ZB")$   = z_sfc
  outfile_S->$("SST")$  = sst_sfc
  outfile_S->$("SLP")$  = mslp_sfc
  outfile_S->$("PS")$   = ps_sfc
  
  system("mv -f " + out_DIR + "S"+year+mo+day+"_"+hr+".nc" + " " + out_DIR + "S"+year+mo+day+"_"+hr)
  
  end do
  
; <><><><><><><><><><><><><><>  
; - GENERATE CONSTANTS FILE -
; <><><><><><><><><><><><><><>  
  if(write_const_file.eq.True)then
  constants_file = "yotc_cst"
  print("CREATING CONSTANTS FILE: " + constants_file)
  outfile_cnst                  = addfile(out_DIR + constants_file + ".nc","c")  
  filedimdef(outfile_cnst, (/"nx","ny","nz"/), (/dims_x,dims_y,dims_z/), (/False,False,False/))   

  lonmin = new(1,"float","No_FillValue")
  lonmax = new(1,"float","No_FillValue")
  latmin = new(1,"float","No_FillValue")
  latmax = new(1,"float","No_FillValue")
  dellon = new(1,"float","No_FillValue")
  dellat = new(1,"float","No_FillValue")  
  pollon = new(1,"float","No_FillValue")  
  pollat = new(1,"float","No_FillValue")      
  starty = new(1,"integer","No_FillValue")    
  startm = new(1,"integer","No_FillValue")      
  startd = new(1,"integer","No_FillValue")      
  starth = new(1,"integer","No_FillValue")      
  starts = new(1,"integer","No_FillValue")        
  dattyp = new(1,"integer","No_FillValue")          
  datver = new(1,"integer","No_FillValue")            
  cstver = new(1,"integer","No_FillValue")  
  
  lonmin = (/min(lon_ml)/)
  lonmax = (/max(lon_ml)/)
  latmin = (/min(lat_ml)/)
  latmax = (/max(lat_ml)/)
  dellon = 0.25
  dellat = 0.25
  pollon = 0.0
  pollat = 90.0
  starty = (/toint(sfc_files_date_utc(ml_ind_files(0),0))/)
  startm = (/toint(sfc_files_date_utc(ml_ind_files(0),1))/)
  startd = (/toint(sfc_files_date_utc(ml_ind_files(0),2))/)
  starth = (/toint(sfc_files_date_utc(ml_ind_files(0),3))/)
  starts = (/toint(sfc_files_date_utc(ml_ind_files(0),4))/)
  dattyp = 2
  datver = 0
  cstver = 0
   
  lonmin!0 = "ncl_scalar"              
  copy_VarMeta(lonmin,lonmax)       
  copy_VarMeta(lonmin,latmin)       
  copy_VarMeta(lonmin,latmax)       
  copy_VarMeta(lonmin,dellon)       
  copy_VarMeta(lonmin,dellat)       
  copy_VarMeta(lonmin,pollon)       
  copy_VarMeta(lonmin,pollat)
         
  starty!0 = "ncl_scalar"                       
  copy_VarMeta(starty,startm)       
  copy_VarMeta(starty,startd)       
  copy_VarMeta(starty,starth)       
  copy_VarMeta(starty,starts)       
  copy_VarMeta(starty,dattyp)       
  copy_VarMeta(starty,datver)       
  copy_VarMeta(starty,cstver)       
    
  outfile_cnst->$("lonmin")$    = lonmin
  outfile_cnst->$("lonmax")$    = lonmax  
  outfile_cnst->$("latmin")$    = latmin
  outfile_cnst->$("latmax")$    = latmax
  outfile_cnst->$("dellon")$    = dellon
  outfile_cnst->$("dellat")$    = dellat
  outfile_cnst->$("pollon")$    = pollon
  outfile_cnst->$("pollat")$    = pollat
  outfile_cnst->$("starty")$    = starty
  outfile_cnst->$("startm")$    = startm
  outfile_cnst->$("startd")$    = startd
  outfile_cnst->$("starth")$    = starth
  outfile_cnst->$("starts")$    = starts
  outfile_cnst->$("dattyp")$    = dattyp
  outfile_cnst->$("datver")$    = datver
  outfile_cnst->$("cstver")$    = cstver
  outfile_cnst->$("aklev")$             = aklev
  outfile_cnst->$("bklev")$             = bklev
  outfile_cnst->$("aklay")$             = aklay
  outfile_cnst->$("bklay")$             = bklay
  
  system("mv -f " + out_DIR + constants_file + ".nc" + " " + out_DIR + constants_file)  
  end if  
  
  end