updates to real-time for new GFS format, SST

atm_to_cam/atm_to_cam.ncl

@@ -74,7 +74,7 @@ print("Loading variables")
 if (datasource .eq. "GFS") then
   ps = grb_file->PRES_P0_L1_GLL0(:,:)
   p0   = 100000.
-  pblh = grb_file->HPBL_P0_L1_GLL0(:,:)
+  ;pblh = grb_file->HPBL_P0_L1_GLL0(:,:)
   t_gfs = grb_file->TMP_P0_L100_GLL0(:,:,:)
   u_gfs = grb_file->UGRD_P0_L100_GLL0(:,:,:)
   v_gfs = grb_file->VGRD_P0_L100_GLL0(:,:,:)
@@ -102,7 +102,7 @@ end if
 if (datasource .eq. "CFSR") then
   ps = grb_file->PRES_P0_L1_GLL0(:,:)
   p0   = 100000.
-  pblh = grb_file->PRES_P0_L1_GLL0(:,:)
+  ;pblh = grb_file->PRES_P0_L1_GLL0(:,:)
   t_gfs = grb_file->TMP_P0_L100_GLL0(:,:,:)
   u_gfs = grb_file->UGRD_P0_L100_GLL0(:,:,:)
   v_gfs = grb_file->VGRD_P0_L100_GLL0(:,:,:)
@@ -130,7 +130,7 @@ if (datasource .eq. "ERAI") then
   ; We need to use short2flt because ERA-Interim data is stored
   ; in short format and requires scale/offset
   ps = short2flt(grb_file->sp(0,:,:))
-  pblh       = short2flt(grb_file->sp(0,:,:))
+  ;pblh       = short2flt(grb_file->sp(0,:,:))
   t_gfs      = short2flt(grb_file->t(0,:,:,:))
   u_gfs      = short2flt(grb_file->u(0,:,:,:))
   v_gfs      = short2flt(grb_file->v(0,:,:,:))
@@ -147,7 +147,7 @@ print("Max Q: "+max(q_gfs)+"   min Q: "+min(q_gfs))
 print("Max PS: "+max(ps)+"   min PS: "+min(ps))
 print("Max CLDICE: "+max(cldice_gfs)+"   min CLDICE: "+min(cldice_gfs))
 print("Max CLDLIQ: "+max(cldliq_gfs)+"   min CLDLIQ: "+min(cldliq_gfs))
-print("Max PBLH: "+max(pblh)+"   min PBLH: "+min(pblh))
+;print("Max PBLH: "+max(pblh)+"   min PBLH: "+min(pblh))
 print("=================================================================")
 
 
@@ -197,8 +197,8 @@ t_cam!1    = "lat"
 t_cam!2    = "lon"
 ps!0 = "lat"
 ps!1 = "lon"
-pblh!0 = "lat"
-pblh!1 = "lon"
+;pblh!0 = "lat"
+;pblh!1 = "lon"
 
 copy_VarCoords(t_cam,q_cam)
 copy_VarCoords(t_cam,cldice_cam)
@@ -220,7 +220,7 @@ latCell     = dfile->lat
 Opt         = True
 Opt@CopyVarCoords = False
 ps_fv  = ESMF_regrid_with_weights(ps,wgt_filename,Opt)
-pblh_fv  = ESMF_regrid_with_weights(pblh,wgt_filename,Opt)
+;pblh_fv  = ESMF_regrid_with_weights(pblh,wgt_filename,Opt)
 t_fv  = ESMF_regrid_with_weights(t_cam,wgt_filename,Opt)
 u_fv  = ESMF_regrid_with_weights(u_cam,wgt_filename,Opt)
 v_fv  = ESMF_regrid_with_weights(v_cam,wgt_filename,Opt)
@@ -235,7 +235,7 @@ ncol=dim_sePS(1)
  print("Converting floats to doubles")
  ; Convert floats to doubles
  ps_fv_dbl = todouble(ps_fv)
- pblh_fv_dbl = todouble(pblh_fv)
+ ;pblh_fv_dbl = todouble(pblh_fv)
  t_fv_dbl = todouble(t_fv)
  u_fv_dbl = todouble(u_fv)
  v_fv_dbl = todouble(v_fv)
@@ -244,7 +244,7 @@ ncol=dim_sePS(1)
  cldliq_fv_dbl = todouble(cldliq_fv)
 
  copy_VarMeta(ps_fv,ps_fv_dbl)
- copy_VarMeta(pblh_fv,pblh_fv_dbl)
+ ;copy_VarMeta(pblh_fv,pblh_fv_dbl)
  copy_VarMeta(t_fv,t_fv_dbl)
  copy_VarMeta(u_fv,u_fv_dbl)
  copy_VarMeta(v_fv,v_fv_dbl)
@@ -254,7 +254,7 @@ ncol=dim_sePS(1)
 ; 
 ; ; clean up some stuff
  delete(ps_fv)
- delete(pblh_fv)
+ ;delete(pblh_fv)
  delete(t_fv)
  delete(u_fv)
  delete(v_fv)
@@ -267,8 +267,8 @@ ncol=dim_sePS(1)
  ps_fv_dbl_time = new((/1,ncol/),double)
  ps_fv_dbl_time(0,:) = ps_fv_dbl
 
- pblh_fv_dbl_time = new((/1,ncol/),double)
- pblh_fv_dbl_time(0,:) = pblh_fv_dbl
+ ;pblh_fv_dbl_time = new((/1,ncol/),double)
+ ;pblh_fv_dbl_time(0,:) = pblh_fv_dbl
 
  u_fv_dbl_time = new((/1,numlevels,ncol/),double)
  u_fv_dbl_time(0,:,:) = u_fv_dbl

dates.txt

@@ -1,3 +1,4 @@
+
 2013090600
 2013090600
 2013090600

notliveshell.csh

@@ -8,7 +8,7 @@
 #BSUB -e out.%J                  # error filename
 #BSUB -q geyser                 # queue
 #BSUB -J sewx_driver
-#BSUB -W 23:59                    # wall clock limit
+#BSUB -W 10:00                    # wall clock limit
 #BSUB -P P54048000               # account number
 ##=======================================================================
 
@@ -23,7 +23,11 @@ islive=0 ; echo "islive set to $islive" # 0 no, using historical data - 1 yes, r
 sendplots=0 ; echo "sendplots set to $sendplots" # 0 send plots to external server, no, 1 yes
 machineid=1 ; echo "machineid set to $machineid" # 1 = Yellowstone, 2 = Flux, 3 = Agri
 
-atmDataType=1            # 1 = GFS analysis, 2 = ERA-interim, 3 = CFSR
+### islive=0 (historical)
+### GFS analysis available from 8/2004 onward
+### CFSR analysis available from 1979 to Mar 2011
+atmDataType=1              # 1 = GFS analysis, 2 = ERA-interim, 3 = CFSR
+### Use NOAAOI unless running real-time
 sstDataType=3              # 1 = GDAS, 2 = ERA, 3 = NOAAOI
 numLevels=30               # 30 -> CAM5 physics, 26 -> CAM4 physics
 
@@ -227,6 +231,10 @@ else
   sstcyclestr=$cyclestr
 fi
 
+yestmonthstr=`date --date="yesterday" -u +%m`
+yestdaystr=`date --date="yesterday" -u +%d`
+yestyearstr=`date --date="yesterday" -u +%Y`
+
 echo "The current time is $currtime"
 echo "We are using $yearstr $monthstr $daystr $cyclestr Z ($cyclestrsec seconds) for GFS data"
 echo "We are using $sstyearstr $sstmonthstr $sstdaystr $sstcyclestr Z for SST data"
@@ -276,7 +284,8 @@ then
         then
           rm -f gfs.t*pgrb2f00*
           gfsFTPPath=ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.$yearstr$monthstr$daystr$cyclestr/
-          gfsFTPFile='gfs.t'$cyclestr'z.pgrb2f00'
+         #gfsFTPFile='gfs.t'$cyclestr'z.pgrb2f00'
+          gfsFTPFile='gfs.t'$cyclestr'z.pgrb2.0p50.anl'
           echo "Attempting to download ${gfsFTPPath}${gfsFTPFile}"
         else
           rm -f gfs.t*pgrb2f00*
@@ -349,8 +358,36 @@ then
 ############################### GET SST / NCL ############################### 
 
   if [ ${sstDataType} -eq 1 ] ; then
+      SSTTYPE=GDAS
       ## Pull sea surface temps, need to rename and delete (if necess)
-      echo "Live GFS SST hasn't been updated in a while..." ; exit
+      #echo "Live GFS SST hasn't been updated in a while..." ; exit
+      echo "Getting SST data"
+      cd ${sst_files_path}
+      if [ $islive -ne 0 ] ; then
+        # Here is where we get the "live" GDAS SST file
+        rm -f gdas1*sstgrb*
+        sstFTPPath=ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/sst.${yestyearstr}${yestmonthstr}${yestdaystr}/
+        sstFTPFile='rtgssthr_grb_0.5.grib2'
+        echo "Attempting to download ${sstFTPPath}${sstFTPFile}"
+      else
+        echo "NCEP broke support for historical GDAS, use NOAAOI instead."
+        exit
+      fi
+      ## Scrape for files
+      error=1
+      while [ $error != 0 ]
+      do
+        wget -nv $sstFTPPath$sstFTPFile
+        error=`echo $?`
+        if [ $error -ne 0 ]
+        then
+          echo "Cannot get file, will wait 2 min and scrape again"
+          sleep 120
+        fi
+      done
+      sstFile='gfs_sst_'$yearstr$monthstr$daystr$cyclestr'.grib2'
+      mv ${sstFTPFile} ${sstFile}
+
 #       echo "Getting SST data"
 #       cd ${sst_files_path}
 #       if [ $islive -ne 0 ] ; then
@@ -383,7 +420,8 @@ then
 #       mv $sstFTPFile 'gfs_sst_'$yearstr$monthstr$daystr$cyclestr'.grib2'
   elif [ ${sstDataType} -eq 2 ] ; then  
     echo "ERA SST not quite supported yet..." ; exit 1
-  elif [ ${sstDataType} -eq 3 ] ; then  
+  elif [ ${sstDataType} -eq 3 ] ; then
+    SSTTYPE=NOAAOI
     echo "Using NOAAOI SSTs"
     cd ${sst_files_path}
     sstFile=sst.day.mean.${yearstr}.v2.nc
@@ -419,10 +457,15 @@ then
       done        
     fi
 
+
+  else
+      echo "Incorrect SST data type entered" ; exit 1
+  fi
+
     set +e
     cd ${sst_to_cam_path}
     ncl sst_interp.ncl 'initdate="'${yearstr}${monthstr}${daystr}${cyclestr}'"' \
-      'datasource="NOAAOI"' \
+      'datasource="'${SSTTYPE}'"' \
       'sstDataFile = "'${sst_files_path}/${sstFile}'"' \
       'iceDataFile = "'${sst_files_path}/${iceFile}'"' \
       'SST_write_file = "'${sstFileIC}'"'
@@ -433,9 +476,6 @@ then
     fi
     echo "SST NCL completed successfully"
     set -e # Turn error checking back on
-  else
-      echo "Incorrect SST data type entered" ; exit 1
-  fi
 
 ############################### ATM NCL ############################### 
 
@@ -501,8 +541,6 @@ then
 
 fi #End debug if statement
 
-exit
-
 ############################### #### ############################### 
 
 cd $path_to_case

sst_to_cam/sst_interp.ncl

@@ -9,6 +9,7 @@ begin
 ;datasource = "GDAS"
 
 print(initdate)
+KtoC=273.15
 
 if (datasource .eq. "GDAS") then
   do_ice = False
@@ -39,10 +40,11 @@ end if
 end if
 
 if (datasource .eq. "GDAS") then
-  sst_file = addfile(sst_file_full,"r")
+  sst_file = addfile(sstDataFile,"r")
   sstlat = sst_file->lat_0
   sstlon = sst_file->lon_0
   sst_gfs = sst_file->TMP_P0_L1_GLL0(:,:)
+  sst_gfs = sst_gfs - KtoC
 else if (datasource .eq. "NOAAOI") then
   sst_file = addfile(sstDataFile,"r")
   ice_file = addfile(iceDataFile,"r")