Add capability to use RAP/HRRR MASSDEN as smoke IC/LBCs

parm/varmap_tables/GSDphys_smoke_var_map.txt

@@ -0,0 +1,28 @@
+dzdt                dzdt      set_to_fill   0         D
+delta_p             delp      skip          0         D
+sphum               sphum     intrp         1E-7      T
+liq_wat             liq_wat   intrp         0         T
+o3mr                o3mr      intrp         1E-7      T
+rainwat             rainwat   intrp         0         T
+ice_wat             ice_wat   intrp         0         T
+snowwat             snowwat   intrp         0         T
+graupel             graupel   intrp         0         T
+ice_nc              ice_nc    intrp         -1.0      T
+rain_nc             rain_nc   intrp         -1.0      T
+water_nc            water_nc  intrp         -1.0      T
+liq_aero            liq_aero  intrp         0         T
+ice_aero            ice_aero  intrp         0         T
+sgs_tke             sgs_tke   intrp         0         T
+vtype               vtype     skip          0         S
+sotype              stype     skip          0         S
+vfrac               vfrac     skip          0         S
+fricv               uustar    skip          0         S
+sfcr                zorl      set_to_fill   0.01      S
+soilw               smc       stop          0         S
+soilt               stc       stop          0         S
+cnwat               cnwat     set_to_fill   0.0       S
+icetk               icetk     set_to_fill   265.0     S
+weasd               weasd     set_to_fill   0.0       S
+snod                snod      set_to_fill   0.0       S
+tprcp               tprcp     set_to_fill   0.0       S
+massden             smoke     set_to_fill   1E-12     T	

parm/varmap_tables/GSDphys_var_map.txt

@@ -24,4 +24,4 @@ cnwat               cnwat     set_to_fill   0.0       S
 icetk               icetk     set_to_fill   265.0     S
 weasd               weasd     set_to_fill   0.0       S
 snod                snod      set_to_fill   0.0       S
-tprcp               tprcp     set_to_fill   0.0       S       	
+tprcp               tprcp     set_to_fill   0.0       S

reg_tests/chgres_cube/13km.conus.rap-smoke.grib2.sh

@@ -0,0 +1,102 @@
+#!/bin/bash
+
+#-----------------------------------------------------------------------------
+# Invoke chgres to create 13-km CONUS coldstart files using RAP GRIB2 data
+# as input.  The coldstart files are then compared to baseline files
+# using the 'nccmp' utility.  This script is run by the machine specific 
+# driver script.
+#-----------------------------------------------------------------------------
+
+set -x
+
+export DATA=$OUTDIR/13km_conus_rap-smoke_grib2
+rm -fr $DATA
+
+export CRES=778
+export KMRES=13km
+export FIXfv3=${HOMEreg}/fix/RRFS_CONUS_${KMRES}
+
+export COMIN=${HOMEreg}/input_data/rap.grib2
+
+export GRIB2_FILE_INPUT=1921221000900
+export VCOORD_FILE=${HOMEufs}/fix/am/global_hyblev.l64.txt
+export VARMAP_FILE=${HOMEufs}/parm/varmap_tables/GSDphys_smoke_var_map.txt
+export INPUT_TYPE='grib2'
+export CONVERT_NST=".false."
+export OROG_FILES_TARGET_GRID="C778_oro_data.tile7.nc"
+export REGIONAL=1
+export HALO_BLEND=0
+export HALO_BNDY=4 
+export CDATE=2019080100
+export EXTERNAL_MODEL="RAP"
+export NSOILL_OUT=9
+export TRACERS_TARGET='"NULL"'
+export TRACERS_INPUT='"NULL"'
+export SOTYP_FROM_CLIMO=.false.
+export VGTYP_FROM_CLIMO=.false.
+export VGFRC_FROM_CLIMO=.true.
+export MINMAX_VGFRC_FROM_CLIMO=.true.
+export TG3_FROM_SOIL=.true.
+export LAI_FROM_CLIMO=.true.
+export GEOGRID_FILE_INPUT=${HOMEufs}/fix/am/geo_em.d01.nc_RAPX
+
+export OMP_NUM_THREADS_CH=${OMP_NUM_THREADS:-1}
+
+NCCMP=${NCCMP:-$(which nccmp)}
+
+#-----------------------------------------------------------------------------
+# Invoke chgres program.
+#-----------------------------------------------------------------------------
+
+echo "Starting at: " `date`
+
+${HOMEufs}/ush/chgres_cube.sh
+
+iret=$?
+if [ $iret -ne 0 ]; then
+  set +x
+  echo "<<< 13-km CONUS RAP-SMOKE W/ GSD PHYSICS AND SFC FROM FILE GRIB2 TEST FAILED. <<<"
+  exit $iret
+fi
+
+echo "Ending at: " `date`
+
+#-----------------------------------------------------------------------------
+# Compare output from chgres to baseline set of data.
+#
+# orion's nccmp utility does not work with the netcdf
+# required to run ufs_utils.  So swap it.
+#-----------------------------------------------------------------------------
+
+machine=${machine:-NULL}
+if [ $machine == 'orion' ]; then
+  module unload netcdfp/4.7.4.release
+  module load netcdf/4.7.2
+fi
+
+cd $DATA
+
+test_failed=0
+for files in *.nc
+do
+  if [ -f $files ]; then
+    echo CHECK $files
+    $NCCMP -dmfqS $files $HOMEreg/baseline_data/13km_conus_rap_grib2/$files
+    iret=$?
+    if [ $iret -ne 0 ]; then
+      test_failed=1
+    fi
+  fi
+done
+
+set +x
+if [ $test_failed -ne 0 ]; then
+  echo "<<< 13-km CONUS RAP-SMOKE W/ GSD PHYSICS AND SFC FROM FILE GRIB2 TEST FAILED. >>>"
+  if [ "$UPDATE_BASELINE" = "TRUE" ]; then
+    $HOMEufs/reg_tests/update_baseline.sh $HOMEreg "13km_conus_rap-smoke_grib2" $commit_num
+  fi
+else
+  echo "<<< 13-km CONUS RAP-SMOKE W/ GSD PHYSICS AND SFC FROM FILE GRIB2 TEST PASSED. >>>"
+fi
+
+exit 0

reg_tests/chgres_cube/driver.jet.sh

@@ -183,14 +183,23 @@ export OMP_NUM_THREADS=1
 TEST13=$(sbatch --parsable --partition=xjet --nodes=1 --ntasks-per-node=6 -t 0:05:00 -A $PROJECT_CODE -q $QUEUE -J c96.gefs.grib2 \
       --exclusive -o $LOG_FILE -e $LOG_FILE ./c96.gefs.grib2.sh)
 
+#-----------------------------------------------------------------------------
+# Initialize CONUS 13-KM USING RAP-SMOKE GRIB2 file WITH GSD PHYSICS .
+#-----------------------------------------------------------------------------
+
+LOG_FILE=consistency.log14
+export OMP_NUM_THREADS=1   # should match cpus-per-task
+TEST14=$(sbatch --parsable --partition=xjet --ntasks-per-node=6 --nodes=1 -t 0:05:00 -A $PROJECT_CODE -q $QUEUE -J 13km.conus.rap-smoke.grib2.conus \
+      --exclusive -o $LOG_FILE -e $LOG_FILE ./13km.conus.rap-smoke.grib2.sh)
+
 #-----------------------------------------------------------------------------
 # Create summary log.
 #-----------------------------------------------------------------------------
 
 LOG_FILE=consistency.log
 sbatch --partition=xjet --nodes=1  -t 0:01:00 -A $PROJECT_CODE -J chgres_summary -o $LOG_FILE -e $LOG_FILE \
        --open-mode=append -q $QUEUE -d\
-       afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6:$TEST7:$TEST8:$TEST9:$TEST10:$TEST11:$TEST12:$TEST13 << EOF
+       afterok:$TEST1:$TEST2:$TEST3:$TEST4:$TEST5:$TEST6:$TEST7:$TEST8:$TEST9:$TEST10:$TEST11:$TEST12:$TEST13:$TEST14 << EOF
 #!/bin/bash
 grep -a '<<<' $LOG_FILE*  > $SUM_FILE
 EOF

sorc/chgres_cube.fd/atm_input_data.F90

@@ -1990,7 +1990,7 @@ subroutine read_input_atm_grib2_file(localpet)
 
  integer, intent(in)                   :: localpet
 
- integer, parameter                    :: ntrac_max=14
+ integer, parameter                    :: ntrac_max=15
  integer, parameter                    :: max_levs=1000
 
  character(len=300)                    :: the_file
@@ -2027,7 +2027,8 @@ subroutine read_input_atm_grib2_file(localpet)
  real(esmf_kind_r8), allocatable       :: rlevs(:)
  real(esmf_kind_r4), allocatable       :: dummy2d(:,:)
  real(esmf_kind_r8), allocatable       :: dummy3d(:,:,:), dummy2d_8(:,:),&
-                                          u_tmp_3d(:,:,:), v_tmp_3d(:,:,:)
+                                          u_tmp_3d(:,:,:), v_tmp_3d(:,:,:),&
+                                          dummy3d_pres(:,:,:)
  real(esmf_kind_r8), pointer           :: presptr(:,:,:), psptr(:,:),tptr(:,:,:), &
                                           qptr(:,:,:), wptr(:,:,:),  &
                                           uptr(:,:,:), vptr(:,:,:)
@@ -2042,18 +2043,19 @@ subroutine read_input_atm_grib2_file(localpet)
 
  tracers(:) = "NULL"
 
- trac_names_oct10 = (/1,  1,  14,  1,  1,  1,  1, 6,  6,   1,  6,  13,  13, 2 /)
- trac_names_oct11 = (/0, 22, 192, 23, 24, 25, 32, 1, 29, 100, 28, 193, 192, 2 /)
+ trac_names_oct10 = (/1,  1,  14,  1,  1,  1,  1, 6,  6,   1,  6,  13,  13, 2, 20 /)
+ trac_names_oct11 = (/0, 22, 192, 23, 24, 25, 32, 1, 29, 100, 28, 193, 192, 2,  0 /)
+
 
  trac_names_vmap = (/"sphum   ", "liq_wat ", "o3mr    ", "ice_wat ", &
                      "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc  ", &
                      "rain_nc ", "water_nc", "liq_aero", "ice_aero", &
-                     "sgs_tke "/)
+                     "sgs_tke ", "massden "/)
 
  tracers_default = (/"sphum   ", "liq_wat ", "o3mr    ", "ice_wat ", &
                      "rainwat ", "snowwat ", "graupel ", "cld_amt ", "ice_nc  ", &
                      "rain_nc ", "water_nc", "liq_aero", "ice_aero", &
-                     "sgs_tke "/)
+                     "sgs_tke ", "smoke   "/)
 
  the_file = trim(data_dir_input_grid) // "/" // trim(grib2_file_input_grid)
 
@@ -2405,11 +2407,19 @@ subroutine read_input_atm_grib2_file(localpet)
    allocate(dummy2d(i_input,j_input))
    allocate(dummy2d_8(i_input,j_input))
    allocate(dummy3d(i_input,j_input,lev_input))
+   if (trim(external_model) .eq. 'RAP' .or. &  ! for smoke conversion
+       trim(external_model) .eq. 'HRRR' ) then
+      allocate(dummy3d_pres(i_input,j_input,lev_input))
+    endif
    allocate(dum2d_1(i_input,j_input))
  else
    allocate(dummy2d(0,0))
    allocate(dummy2d_8(0,0))
    allocate(dummy3d(0,0,0))
+   if (trim(external_model) .eq. 'RAP' .or. &  ! for smoke conversion
+       trim(external_model) .eq. 'HRRR' ) then
+      allocate(dummy3d_pres(0,0,0))
+   endif
    allocate(dum2d_1(0,0))
  endif
 
@@ -2469,12 +2479,61 @@ subroutine read_input_atm_grib2_file(localpet)
    call get_var_cond(vname,this_miss_var_method=method, this_miss_var_value=value, &
                        this_field_var_name=tmpstr,loc=varnum)
 
-   if (n==1 .and. .not. hasspfh) then 
+   if (n==1 .and. .not. hasspfh .or. &
+     ( (trim(external_model) .eq. 'RAP' .or. &  ! for smoke conversion
+       trim(external_model) .eq. 'HRRR' ) .and. &
+       tracers_input_vmap(n) == trac_names_vmap(15) )) then 
      print*,"- CALL FieldGather TEMPERATURE." 
      call ESMF_FieldGather(temp_input_grid,dummy3d,rootPet=0, tile=1, rc=rc)
      if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
      call error_handler("IN FieldGet", rc) 
    endif
+
+   if ( (trim(external_model) .eq. 'RAP' .or. &  ! for smoke conversion
+       trim(external_model) .eq. 'HRRR' ) .and. &
+       tracers_input_vmap(n) == trac_names_vmap(15)) then
+
+    if (localpet == 0) then
+
+    print*,"- READ PRESSURE FOR SMOKE CONVERSION."
+
+    jdisc   = 0     ! search for discipline - meteorological products
+    j = 0           ! search at beginning of file.
+    jpdt    = -9999  ! array of values in product definition template, set towildcard
+    jids    = -9999  ! array of values in identification section, set towildcard
+    jgdt    = -9999  ! array of values in grid definition template, set towildcard
+    jgdtn   = -1     ! search for any grid definition number.
+    jpdtn   =  pdt_num ! Search for the product definition template number.
+    jpdt(1) = 3      ! Sect4/oct 10 - parameter category - mass
+    jpdt(2) = 0      ! Sect4/oct 11 - parameter number - pressure
+    jpdt(10) = octet_23 ! Sect4/oct 23 - type of level.
+    unpack=.true.
+
+    do vlev = 1, lev_input
+
+      jpdt(12) = nint(rlevs(vlev))
+      call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
+             unpack, k, gfld, iret)
+      if (iret /= 0) then
+        call error_handler("READING IN PRESSURE AT LEVEL"//trim(slevs(vlev)),iret)
+      endif
+
+      dum2d_1 = reshape(gfld%fld, (/i_input,j_input/) )
+
+      dummy3d_pres(:,:,vlev) = dum2d_1
+
+    enddo
+
+    endif  ! localpet == 0
+
+   endif   ! read pressure for smoke conversion
+
+   if (tracers_input_vmap(n) == trac_names_vmap(15) .and. &
+       (trim(external_model) .ne. 'RAP' .and. &  ! for smoke conversion
+       trim(external_model) .ne. 'HRRR' ) ) then
+       cycle ! Do not process smoke for non RAP/HRRR
+   endif
+
 
    if (localpet == 0) then
 
@@ -2569,6 +2628,16 @@ subroutine read_input_atm_grib2_file(localpet)
         end if
       endif
 
+      ! Convert smoke from mass density (RAP/HRRR = kg/m^3) to mixing ratio (ug/kg)
+      if ( tracers_input_vmap(n) == trac_names_vmap(15) ) then
+         do i = 1, i_input
+            do j = 1, j_input
+               dummy2d(i,j) = dummy2d(i,j) * 1.0d9 * &
+                              (287.05 * dummy3d(i,j,vlev) / dummy3d_pres(i,j,vlev))
+            enddo
+         enddo
+      endif
+
        dummy3d(:,:,vlev) = real(dummy2d,esmf_kind_r8)
 
      enddo !vlev
@@ -2881,6 +2950,7 @@ subroutine read_input_atm_grib2_file(localpet)
  endif
 
  deallocate(dummy3d, dum2d_1) 
+ if (allocated(dummy3d_pres)) deallocate(dummy3d_pres)
 
 !---------------------------------------------------------------------------
 ! Convert from 2-d to 3-d component winds.