Implement SNICAR v3 snow albedo model into Noah-MP

drivers/hrldas/ConfigVarInTransferMod.F90

@@ -95,6 +95,23 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
     noahmp%config%domain%RunoffSlopeType             = NoahmpIO%SLOPETYP
     noahmp%config%domain%DepthSoilTempBottom         = NoahmpIO%ZBOT_TABLE
 
+    !SNICAR
+    if (noahmp%config%nmlist%OptSnowAlbedo == 3 )then
+       noahmp%config%nmlist%OptSnicarSnowShape          = NoahmpIO%SNICAR_SNOWSHAPE_OPT
+       noahmp%config%nmlist%OptSnicarRTSolver           = NoahmpIO%SNICAR_RTSOLVER_OPT
+       noahmp%config%nmlist%FlagSnicarSnowBCIntmix      = NoahmpIO%SNICAR_SNOWBC_INTMIX
+       noahmp%config%nmlist%FlagSnicarSnowDustIntmix    = NoahmpIO%SNICAR_SNOWDUST_INTMIX
+       noahmp%config%nmlist%FlagSnicarUseAerosol        = NoahmpIO%SNICAR_USE_AEROSOL
+       noahmp%config%nmlist%FlagSnicarUseOC             = NoahmpIO%SNICAR_USE_OC                
+       noahmp%config%nmlist%FlagSnicarAerosolReadTable  = NoahmpIO%SNICAR_AEROSOL_READTABLE     
+
+       noahmp%config%domain%idx_T_max                   = NoahmpIO%idx_T_max
+       noahmp%config%domain%idx_Tgrd_max                = NoahmpIO%idx_Tgrd_max
+       noahmp%config%domain%idx_rhos_max                = NoahmpIO%idx_rhos_max
+       noahmp%config%domain%NumSnicarRadBand            = NoahmpIO%snicar_numrad_snw
+       noahmp%config%domain%idx_Mie_snw_mx              = NoahmpIO%idx_Mie_snw_mx
+    endif
+
     ! the following initialization cannot be done in ConfigVarInitMod
     ! because the NumSoilLayer and NumSnowLayerMax are initialized with input values in this module
     if ( .not. allocated(noahmp%config%domain%DepthSoilLayer) )          &

drivers/hrldas/EnergyVarInTransferMod.F90

@@ -40,7 +40,9 @@ subroutine EnergyVarInTransfer(noahmp, NoahmpIO)
               FlagUrban       => noahmp%config%domain%FlagUrban       ,&
               NumSnowLayerMax => noahmp%config%domain%NumSnowLayerMax ,&
               NumSoilLayer    => noahmp%config%domain%NumSoilLayer    ,&
-              NumSwRadBand    => noahmp%config%domain%NumSwRadBand     &
+              NumSwRadBand    => noahmp%config%domain%NumSwRadBand    ,&
+              NumSnicarRadBand=> noahmp%config%domain%NumSnicarRadBand,&
+              idx_Mie_snw_mx  => noahmp%config%domain%idx_Mie_snw_mx   &
              )
 ! -------------------------------------------------------------------------
 
@@ -128,6 +130,46 @@ subroutine EnergyVarInTransfer(noahmp, NoahmpIO)
     noahmp%energy%param%AlbedoSoilDry    (1:NumSwRadBand)         = NoahmpIO%ALBDRY_TABLE(SoilColor,1:NumSwRadBand)
     noahmp%energy%param%AlbedoLakeFrz    (1:NumSwRadBand)         = NoahmpIO%ALBLAK_TABLE(1:NumSwRadBand)
     noahmp%energy%param%ScatterCoeffSnow (1:NumSwRadBand)         = NoahmpIO%OMEGAS_TABLE(1:NumSwRadBand)
+    !SNICAR
+    if (noahmp%config%nmlist%OptSnowAlbedo == 3 )then
+       noahmp%energy%state%AlbedoSoilDir       (1:NumSwRadBand)                      = NoahmpIO%ALBSOILDIRXY(I,1:NumSwRadBand,J)
+       noahmp%energy%state%AlbedoSoilDif       (1:NumSwRadBand)                      = NoahmpIO%ALBSOILDIFXY(I,1:NumSwRadBand,J)
+       noahmp%energy%param%flx_wgt_dif         (1:NumSnicarRadBand)                  = NoahmpIO%flx_wgt_dif(1:NumSnicarRadBand)
+       noahmp%energy%param%flx_wgt_dir         (1:NumSnicarRadBand)                  = NoahmpIO%flx_wgt_dir(1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_snw_drc      (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%ss_alb_snw_drc      (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_snw_drc     (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%asm_prm_snw_drc     (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_snw_drc (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%ext_cff_mss_snw_drc (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_snw_dfs      (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%ss_alb_snw_dfs      (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_snw_dfs     (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%asm_prm_snw_dfs     (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_snw_dfs (1:idx_Mie_snw_mx,1:NumSnicarRadBand) = NoahmpIO%ext_cff_mss_snw_dfs (1:idx_Mie_snw_mx,1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_bc1          (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_bc1       (1:NumSnicarRadBand) 
+       noahmp%energy%param%asm_prm_bc1         (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_bc1      (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_bc1     (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_bc1  (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_bc2          (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_bc2       (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_bc2         (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_bc2      (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_bc2     (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_bc2  (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_oc1          (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_oc1       (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_oc1         (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_oc1      (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_oc1     (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_oc1  (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_oc2          (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_oc2       (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_oc2         (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_oc2      (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_oc2     (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_oc2  (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_dst1         (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_dst1      (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_dst1        (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_dst1     (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_dst1    (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_dst1 (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_dst2         (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_dst2      (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_dst2        (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_dst2     (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_dst2    (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_dst2 (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_dst3         (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_dst3      (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_dst3        (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_dst3     (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_dst3    (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_dst3 (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_dst4         (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_dst4      (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_dst4        (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_dst4     (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_dst4    (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_dst4 (1:NumSnicarRadBand)
+       noahmp%energy%param%ss_alb_dst5         (1:NumSnicarRadBand)                  = NoahmpIO%ss_alb_dst5      (1:NumSnicarRadBand)
+       noahmp%energy%param%asm_prm_dst5        (1:NumSnicarRadBand)                  = NoahmpIO%asm_prm_dst5     (1:NumSnicarRadBand)
+       noahmp%energy%param%ext_cff_mss_dst5    (1:NumSnicarRadBand)                  = NoahmpIO%ext_cff_mss_dst5 (1:NumSnicarRadBand)
+    endif
 
     do SoilLayerIndex = 1, size(SoilType)
        noahmp%energy%param%SoilQuartzFrac(SoilLayerIndex)         = NoahmpIO%QUARTZ_TABLE(SoilType(SoilLayerIndex))

drivers/hrldas/EnergyVarOutTransferMod.F90

@@ -38,6 +38,7 @@ subroutine EnergyVarOutTransfer(noahmp, NoahmpIO)
               NumSoilLayer    => noahmp%config%domain%NumSoilLayer    ,&
               NumSnowLayerMax => noahmp%config%domain%NumSnowLayerMax ,&
               NumSnowLayerNeg => noahmp%config%domain%NumSnowLayerNeg ,&
+              NumSwRadBand    => noahmp%config%domain%NumSwRadBand    ,&
               IndicatorIceSfc => noahmp%config%domain%IndicatorIceSfc  &
              )
 !-----------------------------------------------------------------------
@@ -142,6 +143,17 @@ subroutine EnergyVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%TSLB    (I,1:NumSoilLayer,J)       = noahmp%energy%state%TemperatureSoilSnow(1:NumSoilLayer)
     NoahmpIO%TSNOXY  (I,-NumSnowLayerMax+1:0,J) = noahmp%energy%state%TemperatureSoilSnow(-NumSnowLayerMax+1:0)
 
+    !SNICAR
+    if (noahmp%config%nmlist%OptSnowAlbedo == 3 )then
+       NoahmpIO%ALBSOILDIRXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSoilDir(1:NumSwRadBand)
+       NoahmpIO%ALBSOILDIFXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSoilDif(1:NumSwRadBand)
+    endif
+
+    NoahmpIO%ALBSFCDIRXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSfcDir(1:NumSwRadBand)
+    NoahmpIO%ALBSFCDIFXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSfcDif(1:NumSwRadBand)
+    NoahmpIO%ALBSNOWDIRXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSnowDir(1:NumSwRadBand)
+    NoahmpIO%ALBSNOWDIFXY(I,1:NumSwRadBand,J)=noahmp%energy%state%AlbedoSnowDif(1:NumSwRadBand)
+
     ! New Calculation of total Canopy/Stomatal Conductance Based on Bonan et al. (2011), Inverse of Canopy Resistance (below)
     LeafAreaIndSunlit      = max(noahmp%energy%state%LeafAreaIndSunlit, 0.0)
     LeafAreaIndShade       = max(noahmp%energy%state%LeafAreaIndShade, 0.0)

drivers/hrldas/ForcingVarInTransferMod.F90

@@ -62,6 +62,34 @@ subroutine ForcingVarInTransfer(noahmp, NoahmpIO)
     noahmp%forcing%PrecipNonConvRefHeight  = noahmp%forcing%PrecipNonConvRefHeight + PrecipOtherRefHeight
     noahmp%forcing%PrecipSnowRefHeight     = noahmp%forcing%PrecipSnowRefHeight + PrecipOtherRefHeight * NoahmpIO%SR(I,J)
 
+    noahmp%forcing%DirFrac                 = NoahmpIO%DirFrac(I,J)
+    noahmp%forcing%VisFrac                 = NoahmpIO%VisFrac(I,J)
+
+    !SNICAR
+    if (noahmp%config%nmlist%OptSnowAlbedo == 3 )then
+       if (noahmp%config%nmlist%FlagSnicarAerosolReadTable == .true.) then 
+          noahmp%forcing%DepBChydropho           = NoahmpIO%DepBChydropho_TABLE
+          noahmp%forcing%DepBChydrophi           = NoahmpIO%DepBChydrophi_TABLE
+          noahmp%forcing%DepOChydropho           = NoahmpIO%DepOChydropho_TABLE
+          noahmp%forcing%DepOChydrophi           = NoahmpIO%DepOChydrophi_TABLE
+          noahmp%forcing%DepDust1                = NoahmpIO%DepDust1_TABLE
+          noahmp%forcing%DepDust2                = NoahmpIO%DepDust2_TABLE
+          noahmp%forcing%DepDust3                = NoahmpIO%DepDust3_TABLE
+          noahmp%forcing%DepDust4                = NoahmpIO%DepDust4_TABLE
+          noahmp%forcing%DepDust5                = NoahmpIO%DepDust5_TABLE
+       else
+          noahmp%forcing%DepBChydropho           = NoahmpIO%DepBChydrophoXY(I,J)
+          noahmp%forcing%DepBChydrophi           = NoahmpIO%DepBChydrophiXY(I,J)
+          noahmp%forcing%DepOChydropho           = NoahmpIO%DepOChydrophoXY(I,J)
+          noahmp%forcing%DepOChydrophi           = NoahmpIO%DepOChydrophiXY(I,J)
+          noahmp%forcing%DepDust1                = NoahmpIO%DepDust1XY(I,J)
+          noahmp%forcing%DepDust2                = NoahmpIO%DepDust2XY(I,J)
+          noahmp%forcing%DepDust3                = NoahmpIO%DepDust3XY(I,J)
+          noahmp%forcing%DepDust4                = NoahmpIO%DepDust4XY(I,J)
+          noahmp%forcing%DepDust5                = NoahmpIO%DepDust5XY(I,J)
+       endif
+    endif
+
     end associate
 
   end subroutine ForcingVarInTransfer

drivers/hrldas/ForcingVarOutTransferMod.F90

@@ -37,6 +37,10 @@ subroutine ForcingVarOutTransfer(noahmp, NoahmpIO)
     NoahmpIO%FORCWLSM  (I,J) = sqrt(noahmp%forcing%WindEastwardRefHeight**2 + &
                                     noahmp%forcing%WindNorthwardRefHeight**2)
 
+    NoahmpIO%DirFrac   (I,J) = noahmp%forcing%DirFrac
+    NoahmpIO%VisFrac   (I,J) = noahmp%forcing%VisFrac
+
+
     end associate
 
   end subroutine ForcingVarOutTransfer

drivers/hrldas/NoahmpDriverMainMod.F90

@@ -199,7 +199,6 @@ subroutine NoahmpDriverMain(NoahmpIO)
                  noahmp%config%domain%IndicatorIceSfc = 0   ! land soil point.
                  call NoahmpMain(noahmp)
              endif ! glacial split ends 
-
              !---------------------------------------------------------------------
              !  Transfer 1-D Noah-MP column variables to 2-D output variables
              !---------------------------------------------------------------------