phys/module_sf_noah_seaice_drv.F

   1 module module_sf_noah_seaice_drv
   2 #if defined(mpas)
   3 use mpas_atmphys_utilities, only: physics_message,physics_error_fatal
   4 #define FATAL_ERROR(M) call physics_error_fatal( M )
   5 #define WRITE_MESSAGE(M) call physics_message( M )
   6 #else
   7 use module_wrf_error
   8 #define FATAL_ERROR(M) call wrf_error_fatal( M )
   9 #define WRITE_MESSAGE(M) call wrf_message( M )
  10 #endif
  11   use module_sf_noah_seaice
  12   implicit none
  13 contains
  14   subroutine seaice_noah( SEAICE_ALBEDO_OPT, SEAICE_ALBEDO_DEFAULT, SEAICE_THICKNESS_OPT, &
  15        &                  SEAICE_THICKNESS_DEFAULT, SEAICE_SNOWDEPTH_OPT,               &
  16        &                  SEAICE_SNOWDEPTH_MAX, SEAICE_SNOWDEPTH_MIN,                   &
  17        &                  T3D, QV3D, P8W3D, DZ8W, NUM_SOIL_LAYERS, DT, FRPCPN, SR,      &
  18        &                  GLW, SWDOWN, RAINBL, SNOALB2D, QGH, XICE, XICE_THRESHOLD,     &
  19        &                  ALBSI, ICEDEPTH, SNOWSI,                                      &
  20        &                  TSLB, EMISS, ALBEDO, Z02D, TSK, SNOW, SNOWC, SNOWH2D, &
  21        &                  CHS, CHS2, CQS2,                                              &
  22        &                  RIB, ZNT, LH, HFX, QFX, POTEVP, GRDFLX, QSFC, ACSNOW,         &
  23        &                  ACSNOM, SNOPCX, SFCRUNOFF, NOAHRES,                           &
  24        &                  SF_URBAN_PHYSICS, B_T_BEP, B_Q_BEP, RHO,                      &
  25        &                  IDS, IDE, JDS, JDE, KDS, KDE,                                 &
  26        &                  IMS, IME, JMS, JME, KMS, KME,                                 &
  27        &                  ITS, ITE, JTS, JTE, KTS, KTE  )
  28 #if defined(wrfmodel)
  29 #if (NMM_CORE != 1)
  30     USE module_state_description, ONLY : NOAHUCMSCHEME
  31     USE module_state_description, ONLY : BEPSCHEME
  32     USE module_state_description, ONLY : BEP_BEMSCHEME
  33 #endif
  34 #endif
  35     implicit none
  36
  37     INTEGER, INTENT(IN)       ::               SEAICE_ALBEDO_OPT
  38     REAL   , INTENT(IN)       ::               SEAICE_ALBEDO_DEFAULT
  39     INTEGER, INTENT(IN)       ::               SEAICE_THICKNESS_OPT
  40     REAL,    INTENT(IN)       ::               SEAICE_THICKNESS_DEFAULT
  41     INTEGER, INTENT(IN)       ::               SEAICE_SNOWDEPTH_OPT
  42     REAL,    INTENT(IN)       ::               SEAICE_SNOWDEPTH_MAX
  43     REAL,    INTENT(IN)       ::               SEAICE_SNOWDEPTH_MIN
  44
  45     INTEGER, INTENT(IN)       ::                            IDS, &
  46          &                                                  IDE, &
  47          &                                                  JDS, &
  48          &                                                  JDE, &
  49          &                                                  KDS, &
  50          &                                                  KDE
  51
  52     INTEGER, INTENT(IN)       ::                            IMS, &
  53          &                                                  IME, &
  54          &                                                  JMS, &
  55          &                                                  JME, &
  56          &                                                  KMS, &
  57          &                                                  KME
  58
  59     INTEGER, INTENT(IN)       ::                            ITS, &
  60          &                                                  ITE, &
  61          &                                                  JTS, &
  62          &                                                  JTE, &
  63          &                                                  KTS, &
  64          &                                                  KTE
  65
  66     REAL,    DIMENSION( ims:ime, kms:kme, jms:jme )            , &
  67          &   INTENT (IN)      ::                            T3D, &
  68          &                                                 QV3D, &
  69          &                                                P8W3D, &
  70          &                                                 DZ8W
  71
  72     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
  73          &   INTENT (IN)      ::                             SR, &
  74          &                                                  GLW, &
  75          &                                                  QGH, &
  76          &                                               SWDOWN, &
  77          &                                               RAINBL, &
  78          &                                             SNOALB2D, &
  79          &                                                 XICE, &
  80          &                                                  RIB
  81
  82     LOGICAL, INTENT (IN)      :: FRPCPN
  83     REAL   , INTENT (IN)      :: DT
  84     INTEGER, INTENT (IN)      :: NUM_SOIL_LAYERS
  85     REAL   , INTENT (IN)      :: XICE_THRESHOLD
  86
  87     REAL,     DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
  88          INTENT(INOUT)   ::                            TSLB
  89
  90     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
  91          &   INTENT (INOUT)   ::                          EMISS, &
  92          &                                               ALBEDO, &
  93          &                                                ALBSI, &
  94          &                                                 Z02D, &
  95          &                                                 SNOW, &
  96          &                                                  TSK, &
  97          &                                                SNOWC, &
  98          &                                              SNOWH2D, &
  99          &                                                  CHS, &
 100          &                                                 CQS2
 101
 102     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 103          &   INTENT (OUT)     ::                            HFX, &
 104          &                                                   LH, &
 105          &                                                  QFX, &
 106          &                                                  ZNT, &
 107          &                                               POTEVP, &
 108          &                                               GRDFLX, &
 109          &                                                 QSFC, &
 110          &                                               ACSNOW, &
 111          &                                               ACSNOM, &
 112          &                                               SNOPCX, &
 113          &                                            SFCRUNOFF, &
 114          &                                              NOAHRES, &
 115          &                                                 CHS2
 116
 117     REAL,    DIMENSION( ims:ime, jms:jme )                      ,&
 118          &   INTENT(INOUT)    ::                         SNOWSI
 119
 120     REAL,    DIMENSION( ims:ime, jms:jme )                     , &
 121          &   INTENT (INOUT)   ::                        ICEDEPTH
 122
 123     INTEGER, INTENT (IN)      ::               SF_URBAN_PHYSICS
 124     REAL,    OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme )  , &
 125          &   INTENT (INOUT)   ::                        B_Q_BEP, &
 126          &                                              B_T_BEP
 127     REAL,    OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme )  , &
 128          &   INTENT (IN)      ::                            RHO
 129
 130     INTEGER :: I
 131     INTEGER :: J
 132     REAL    :: FFROZP
 133     REAL    :: ZLVL
 134     INTEGER :: NSOIL
 135     REAL    :: LWDN
 136     REAL    :: SOLNET
 137     REAL    :: SFCPRS
 138     REAL    :: PRCP
 139     REAL    :: SFCTMP
 140     REAL    :: Q2
 141     REAL    :: TH2
 142     REAL    :: Q2SAT
 143     REAL    :: DQSDT2
 144     REAL    :: SNOALB
 145     REAL    :: TBOT
 146     REAL    :: SITHICK
 147
 148     REAL    :: ALBEDOK
 149     REAL    :: ALBBRD
 150     REAL    :: Z0BRD
 151     REAL    :: EMISSI
 152     REAL    :: T1
 153     REAL, DIMENSION(1:NUM_SOIL_LAYERS)::  STC
 154     REAL    :: SNOWH
 155     REAL    :: SNEQV
 156     REAL    :: CH
 157     REAL    :: SNCOVR
 158     REAL    :: RIBB
 159
 160     REAL    :: Z0
 161     REAL    :: ETA
 162     REAL    :: SHEAT
 163     REAL    :: ETA_KINEMATIC
 164     REAL    :: FDOWN
 165     REAL    :: ESNOW
 166     REAL    :: DEW
 167     REAL    :: ETP
 168     REAL    :: SSOIL
 169     REAL    :: FLX1
 170     REAL    :: FLX2
 171     REAL    :: FLX3
 172     REAL    :: SNOMLT
 173     REAL    :: RUNOFF1
 174     REAL    :: Q1
 175
 176     REAL    :: APES
 177     REAL    :: APELM
 178     REAL    :: PSFC
 179     REAL    :: SFCTSNO
 180     REAL    :: E2SAT
 181     REAL    :: Q2SATI
 182     INTEGER :: NS
 183     REAL    :: FDTW
 184     REAL    :: FDTLIW
 185     REAL    :: ALBEDOSI
 186     REAL    :: SNOWONSI
 187     REAL, PARAMETER  :: CAPA   = R_D / CP
 188     REAL, PARAMETER  :: A2     = 17.67
 189     REAL, PARAMETER  :: A3     = 273.15
 190     REAL, PARAMETER  :: A4     = 29.65
 191     REAL, PARAMETER  :: A23M4  = A2 * ( A3 - A4 )
 192     REAL, PARAMETER  :: ROW    = 1.E3
 193     REAL, PARAMETER  :: ELIW   = XLF
 194     REAL, PARAMETER  :: ROWLIW = ROW * ELIW
 195
 196     CHARACTER(len=80) :: message
 197
 198     FDTLIW = DT / ROWLIW
 199     FDTW   = DT / ( XLV * RHOWATER )
 200
 201     NSOIL  = NUM_SOIL_LAYERS
 202
 203     SEAICE_JLOOP : do J = JTS, JTE
 204        SEAICE_ILOOP : do I = ITS, ITE
 205
 206           ! Skip the points that are not sea-ice points.
 207           IF ( XICE(I,J) < XICE_THRESHOLD ) THEN
 208               IF ( SEAICE_THICKNESS_OPT == 1 ) THEN
 209                   ICEDEPTH(I,J) = 0.0
 210               ENDIF
 211               IF ( SEAICE_SNOWDEPTH_OPT == 1 ) THEN
 212                   SNOWSI(I,J) = 0.0
 213               ENDIF
 214               CYCLE SEAICE_ILOOP
 215           ENDIF
 216
 217           SELECT CASE ( SEAICE_THICKNESS_OPT )
 218           CASE DEFAULT
 219               WRITE(message,'("Namelist value for SEAICE_THICKNESS_OPT not recognized: ",I6)') SEAICE_THICKNESS_OPT
 220               FATAL_ERROR(message)
 221           CASE (0)
 222               ! Use uniform sea-ice thickness.
 223               SITHICK = SEAICE_THICKNESS_DEFAULT
 224           CASE (1)
 225               ! Use the sea-ice as read in from the input files.
 226               ! Limit the to between 0.10 and 10.0 m.
 227               IF ( ICEDEPTH(I,J) < -1.E6 ) THEN
 228                   WRITE_MESSAGE("Field ICEDEPTH not found in input files.")
 229                   WRITE_MESSAGE(".... Namelist SEAICE_THICKNESS_OPT=1 requires ICEDEPTH field.")
 230                   WRITE_MESSAGE(".... Try namelist option SEAICE_THICKNESS_OPT=0.")
 231                   FATAL_ERROR("SEAICE_THICKNESS_OPT")
 232               ENDIF
 233               SITHICK = MIN ( MAX ( 0.10 , ICEDEPTH(I,J) ) , 10.0 )
 234               ICEDEPTH(I,J) = SITHICK
 235           END SELECT
 236
 237           SFCTMP = T3D(I,1,J)
 238           T1     = TSK(I,J)
 239           IF ( SEAICE_ALBEDO_OPT == 2 ) THEN
 240               IF ( ALBSI(I,J) < -1.E6 ) THEN
 241                   FATAL_ERROR("Field ALBSI not found in input.  Field ALBSI is required if SEAICE_ALBEDO_OPT=2")
 242               ENDIF
 243               SNOALB = ALBSI(I,J)
 244               ALBEDO(I,J) = ALBSI(I,J)
 245               ALBEDOK = ALBSI(I,J)
 246               ALBBRD = ALBSI(I,J)
 247               ALBEDOSI = ALBSI(I,J)
 248           ELSE
 249               SNOALB = SNOALB2D(I,J)
 250           ENDIF
 251           ZLVL   = 0.5 * DZ8W(I,1,J)
 252           EMISSI = EMISS(I,J)               ! But EMISSI might change in SFLX_SEAICE
 253           LWDN   = GLW(I,J) * EMISSI        ! But EMISSI might change in SFLX_SEAICE
 254
 255           ! convert snow water equivalent from mm to meter
 256           SNEQV = SNOW(I,J) * 0.001
 257
 258           ! snow depth in meters
 259           SNOWH = SNOWH2D(I,J)
 260           SNCOVR = SNOWC(I,J)
 261
 262           ! Use mid-day albedo to determine net downward solar (no solar zenith angle correction)
 263           SOLNET = SWDOWN(I,J) * (1.-ALBEDO(I,J))   ! But ALBEDO might change after SFLX_SEAICE
 264
 265           ! Pressure in middle of lowest layer.  Why don't we use the true surface pressure?
 266           ! Are there places where we would need to use the true surface pressure?
 267           SFCPRS = ( P8W3D(I,KTS+1,j) + P8W3D(I,KTS,J) ) * 0.5
 268
 269           ! surface pressure
 270           PSFC   = P8W3D(I,1,J)
 271
 272           ! Convert lowest model level humidity from mixing ratio to specific humidity
 273           Q2     = QV3D(I,1,J) / ( 1.0 + QV3D(I,1,J) )
 274
 275           ! Calculate TH2 via Exner function
 276           APES   = ( 1.E5 / PSFC )   ** CAPA
 277           APELM  = ( 1.E5 / SFCPRS ) ** CAPA
 278           TH2    = ( SFCTMP * APELM ) / APES
 279
 280           ! Q2SAT is specific humidity
 281           Q2SAT  = QGH(I,J) / ( 1.0 + QGH(I,J) )
 282           DQSDT2 = Q2SAT * A23M4 / ( SFCTMP - A4 ) ** 2
 283
 284           SELECT CASE ( SEAICE_SNOWDEPTH_OPT )
 285           CASE DEFAULT
 286
 287               WRITE(message,'("Namelist value for SEAICE_SNOWDEPTH_OPT not recognized: ",I6)') SEAICE_SNOWDEPTH_OPT
 288               FATAL_ERROR(message)
 289
 290           CASE ( 0 )
 291
 292               ! Minimum and maximum bounds on snow depth are enforced in SFLX_SEAICE
 293
 294           CASE ( 1 )
 295
 296               ! Snow depth on sea ice comes from a 2D array, SNOWSI, bounded by user-specified
 297               ! minimum and maximum values.  No matter what anybody else says about snow
 298               ! accumulation and melt, we want the snow depth on sea ice to be specified
 299               ! as SNOWSI (bounded by SEAICE_SNOWDEPTH_MIN and SEAICE_SNOWDEPTH_MAX).
 300               SNOWONSI = MAX ( SEAICE_SNOWDEPTH_MIN , MIN ( SNOWSI(I,J) , SEAICE_SNOWDEPTH_MAX ) )
 301               SNEQV = SNOWONSI * 0.3
 302               SNOWH2D(I,J) = SNOWONSI
 303
 304           END SELECT
 305
 306           IF ( SNOW(I,J) .GT. 0.0 ) THEN
 307              ! If snow on surface, use ice saturation properties
 308              SFCTSNO = SFCTMP ! Lowest model Air temperature
 309              E2SAT = 611.2 * EXP ( 6174. * ( 1./273.15 - 1./SFCTSNO ) )
 310              Q2SATI = 0.622 * E2SAT / ( SFCPRS - E2SAT )
 311              Q2SATI = Q2SATI / ( 1.0 + Q2SATI )    ! Convert to specific humidity
 312              ! T1 is skin temperature
 313              IF (T1 .GT. 273.14) THEN
 314                 ! Warm ground temps, weight the saturation between ice and water according to SNOWC
 315                 Q2SAT = Q2SAT * (1.-SNOWC(I,J)) + Q2SATI * SNOWC(I,J)
 316                 DQSDT2 = DQSDT2 * (1.-SNOWC(I,J)) + Q2SATI * 6174. / (SFCTSNO**2) * SNOWC(I,J)
 317              ELSE
 318                 ! Cold ground temps, use ice saturation only
 319                 Q2SAT = Q2SATI
 320                 DQSDT2 = Q2SATI * 6174. / (SFCTSNO**2)
 321              ENDIF
 322              IF ( ( T1 .GT. 273. ) .AND. ( SNOWC(I,J) .GT. 0.0 ) ) THEN
 323                 ! If (SNOW > 0) can we have (SNOWC <= 0) ?  Perhaps not, so the check on
 324                 ! SNOWC here might be superfluous.
 325                 DQSDT2 = DQSDT2 * ( 1. - SNOWC(I,J) )
 326              ENDIF
 327           ENDIF
 328
 329           PRCP = RAINBL(I,J) / DT
 330
 331           ! If "SR" is present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
 332           ! SR from e.g. Ferrier microphysics
 333           ! otherwise define from 1st atmos level temperature
 334
 335           IF (FRPCPN) THEN
 336              FFROZP = SR(I,J)
 337           ELSE
 338              IF (SFCTMP <=  273.15) THEN
 339                 FFROZP = 1.0
 340              ELSE
 341                 FFROZP = 0.0
 342              ENDIF
 343           ENDIF
 344
 345           ! Sea-ice point has deep-level temperature of about -1.8 C
 346           TBOT = 271.36
 347           ! TBOT=273.15  ! appropriate value for lake ice.
 348
 349           ! INTENT(IN) for SFLX_SEAICE, values unchanged by SFLX_SEAICE
 350           !       I           --
 351           !       J           --
 352           !       FFROZP      --
 353           !       DT          --
 354           !       ZLVL        --
 355           !       NSOIL       --
 356           !       LWDN        --
 357           !       SOLNET      --
 358           !       SFCPRS      --
 359           !       PRCP        --
 360           !       SFCTMP      --
 361           !       Q2          --
 362           !       TH2         --
 363           !       Q2SAT       --
 364           !       DQSDT2      --
 365           !       SNOALB      --
 366           !       TBOT        --
 367
 368           Z0BRD  = Z02D(I,J)
 369
 370           DO NS = 1, NSOIL
 371              STC(NS) = TSLB(I,NS,J)
 372           ENDDO
 373
 374           CH = CHS(I,J)
 375           RIBB = RIB(I,J)
 376
 377           ! INTENT(INOUT) for SFLX_SEAICE, values updated by SFLX_SEAICE
 378           !       Z0BRD       --
 379           !       EMISSI      --
 380           !       T1          --
 381           !       STC         --
 382           !       SNOWH       --
 383           !       SNEQV       --
 384           !       SNCOVR      --
 385           !       CH          -- but the result isn't used for anything.
 386           !                      Might as well be intent in to SFLX_SEAICE and changed locally in
 387           !                      that routine?
 388           !       RIBB        -- but the result isn't used for anything.
 389           !                      Might as well be intent in to SFLX_SEAICE and changed locally in
 390           !                      that routine?
 391
 392           ! INTENT(OUT) for SFLX_SEAICE.  Input value should not matter.
 393           Z0               = -1.E36
 394           ETA              = -1.E36
 395           SHEAT            = -1.E36
 396           ETA_KINEMATIC    = -1.E36
 397           FDOWN            = -1.E36  ! Returned value unused.  Might as well be local to SFLX_SEAICE ?
 398           ESNOW            = -1.E36  ! Returned value unused.  Might as well be local to SFLX_SEAICE ?
 399           DEW              = -1.E36  ! Returned value unused.  Might as well be local to SFLX_SEAICE ?
 400           ETP              = -1.E36
 401           SSOIL            = -1.E36
 402           FLX1             = -1.E36
 403           FLX2             = -1.E36
 404           FLX3             = -1.E36
 405           SNOMLT           = -1.E36
 406           RUNOFF1          = -1.E36
 407           Q1               = -1.E36
 408
 409           call sflx_seaice(I, J, SEAICE_ALBEDO_OPT, SEAICE_ALBEDO_DEFAULT,  &    !C
 410                &           SEAICE_SNOWDEPTH_OPT, SEAICE_SNOWDEPTH_MAX,      &    !C
 411                &           SEAICE_SNOWDEPTH_MIN,                            &    !C
 412                &           FFROZP, DT, ZLVL, NSOIL,                         &    !C
 413                &           SITHICK,                                         &
 414                &           LWDN, SOLNET, SFCPRS, PRCP, SFCTMP, Q2,          &    !F
 415                &           TH2, Q2SAT, DQSDT2,                              &    !I
 416                &           SNOALB, TBOT, Z0BRD, Z0, EMISSI,                 &    !S
 417                &           T1, STC, SNOWH, SNEQV, ALBEDOK, CH,              &    !H
 418                &           ALBEDOSI, SNOWONSI,                              &
 419                &           ETA, SHEAT, ETA_KINEMATIC, FDOWN,                &    !O
 420                &           ESNOW, DEW, ETP, SSOIL, FLX1, FLX2, FLX3,        &    !O
 421                &           SNOMLT, SNCOVR,                                  &    !O
 422                &           RUNOFF1, Q1, RIBB)
 423
 424           ! Update our 2d arrays with results from SFLX_SEAICE
 425           ALBEDO(I,J)  = ALBEDOK
 426           EMISS(I,J)   = EMISSI
 427           TSK(I,J)     = T1
 428           Z02D(I,J)    = Z0BRD
 429           SNOWH2D(I,J) = SNOWH
 430           SNOWC(I,J)   = SNCOVR
 431
 432           ! Convert snow water equivalent from (m) back to (mm)
 433           SNOW(I,J)    = SNEQV * 1000.
 434
 435           ! Update our ice temperature array with results from SFLX_SEAICE
 436           DO NS = 1,NSOIL
 437              TSLB(I,NS,J) = STC(NS)
 438           ENDDO
 439
 440           ! Intent (OUT) from SFLX_SEAICE
 441           ZNT(I,J)    = Z0
 442           LH(I,J)     = ETA
 443           HFX(I,J)    = SHEAT
 444           QFX(I,J)    = ETA_KINEMATIC
 445           POTEVP(I,J) = POTEVP(I,J) + ETP*FDTW
 446           GRDFLX(I,J) = SSOIL
 447
 448           ! Exchange Coefficients
 449           CHS2(I,J) = CQS2(I,J)
 450           IF (Q1 .GT. QSFC(I,J)) THEN
 451              CQS2(I,J) = CHS(I,J)
 452           ENDIF
 453
 454           ! Convert QSFC term back to Mixing Ratio.
 455           QSFC(I,J)   = Q1 / ( 1.0 - Q1 )
 456
 457           IF ( SEAICE_SNOWDEPTH_OPT == 1 ) THEN
 458               SNOWSI(I,J) = SNOWONSI
 459           ENDIF
 460
 461           ! Accumulated snow precipitation.
 462           IF ( FFROZP .GT. 0.5 ) THEN
 463              ACSNOW(I,J) = ACSNOW(I,J) + PRCP * DT
 464           ENDIF
 465
 466           ! Accumulated snow melt.
 467           ACSNOM(I,J) = ACSNOM(I,J) + SNOMLT * 1000.
 468
 469           ! Accumulated snow-melt energy.
 470           SNOPCX(I,J) = SNOPCX(I,J) - SNOMLT/FDTLIW
 471
 472           ! Surface runoff
 473           SFCRUNOFF(I,J) = SFCRUNOFF(I,J) + RUNOFF1 * DT * 1000.0
 474
 475           !
 476           ! Residual of surface energy balance terms
 477           !
 478           NOAHRES(I,J) = ( SOLNET + LWDN ) &
 479                &         - SHEAT + SSOIL - ETA &
 480                &         - ( EMISSI * STBOLT * (T1**4) ) &
 481                &         - FLX1 - FLX2 - FLX3
 482 #if defined(wrfmodel)
 483 #if (NMM_CORE != 1)
 484           IF ( ( SF_URBAN_PHYSICS == NOAHUCMSCHEME ) .OR. &
 485                (SF_URBAN_PHYSICS == BEPSCHEME )      .OR. &
 486                ( SF_URBAN_PHYSICS == BEP_BEMSCHEME ) ) THEN
 487              if ( PRESENT (B_T_BEP) ) then
 488                 B_T_BEP(I,1,J)=hfx(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP
 489              endif
 490              if ( PRESENT (B_Q_BEP) ) then
 491                 B_Q_BEP(I,1,J)=qfx(i,j)/dz8w(i,1,j)/rho(i,1,j)
 492              endif
 493           ENDIF
 494 #endif
 495 #endif
 496
 497        enddo SEAICE_ILOOP
 498     enddo SEAICE_JLOOP
 499
 500   end subroutine seaice_noah
 501
 502 end module module_sf_noah_seaice_drv