hydro/Routing/Noah_distr_routing.F90

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 !DJG ------------------------------------------------
  22 !DJG   SUBROUTINE RT_PARM
  23 !DJG ------------------------------------------------
  24
  25 SUBROUTINE RT_PARM(IX,JY,IXRT,JXRT,VEGTYP,RETDP,OVRGH,  &
  26         AGGFACTR)
  27 #ifdef MPP_LAND
  28     use module_mpp_land, only: left_id,down_id,right_id,&
  29         up_id,mpp_land_com_real,MPP_LAND_UB_COM, &
  30         MPP_LAND_LR_COM,mpp_land_com_integer
  31 #endif
  32
  33     IMPLICIT NONE
  34
  35     !DJG -------- DECLARATIONS -----------------------
  36
  37     INTEGER, INTENT(IN) :: IX,JY,IXRT,JXRT,AGGFACTR
  38
  39     INTEGER, INTENT(IN), DIMENSION(IX,JY)       :: VEGTYP
  40     REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)     :: RETDP
  41     REAL, INTENT(OUT), DIMENSION(IXRT,JXRT)     :: OVRGH
  42
  43
  44     !DJG Local Variables
  45
  46     INTEGER     :: I,J,IXXRT,JYYRT
  47     INTEGER :: AGGFACYRT,AGGFACXRT
  48
  49
  50     !DJG Assign RETDP and OVRGH based on VEGTYP...
  51
  52     do J=1,JY
  53         do I=1,IX
  54
  55             do AGGFACYRT=AGGFACTR-1,0,-1
  56                 do AGGFACXRT=AGGFACTR-1,0,-1
  57
  58                     IXXRT=I*AGGFACTR-AGGFACXRT
  59                     JYYRT=J*AGGFACTR-AGGFACYRT
  60 #ifdef MPP_LAND
  61                     if(left_id.ge.0) IXXRT=IXXRT+1
  62                     if(down_id.ge.0) JYYRT=JYYRT+1
  63 #else
  64                     !yw ????
  65                     !       IXXRT=IXXRT+1
  66                     !       JYYRT=JYYRT+1
  67 #endif
  68
  69                     !        if(AGGFACTR .eq. 1) then
  70                     !            IXXRT=I
  71                     !            JYYRT=J
  72                     !        endif
  73
  74
  75
  76                     !DJG Urban, rock, playa, snow/ice...
  77                     IF (VEGTYP(I,J).EQ.1.OR.VEGTYP(I,J).EQ.26.OR.   &
  78                             VEGTYP(I,J).EQ.26.OR.VEGTYP(I,J).EQ.24) THEN
  79                         RETDP(IXXRT,JYYRT)=1.3
  80                         OVRGH(IXXRT,JYYRT)=0.1
  81                         !DJG Wetlands and water bodies...
  82                     ELSE IF (VEGTYP(I,J).EQ.17.OR.VEGTYP(I,J).EQ.18.OR.  &
  83                             VEGTYP(I,J).EQ.19.OR.VEGTYP(I,J).EQ.16) THEN
  84                         RETDP(IXXRT,JYYRT)=10.0
  85                         OVRGH(IXXRT,JYYRT)=0.2
  86                         !DJG All other natural covers...
  87                     ELSE
  88                         RETDP(IXXRT,JYYRT)=5.0
  89                         OVRGH(IXXRT,JYYRT)=0.2
  90                     END IF
  91
  92                 end do
  93             end do
  94
  95         end do
  96     end do
  97 #ifdef MPP_LAND
  98     call MPP_LAND_COM_REAL(RETDP,IXRT,JXRT,99)
  99     call MPP_LAND_COM_REAL(OVRGH,IXRT,JXRT,99)
 100 #endif
 101
 102     !DJG ----------------------------------------------------------------
 103 END SUBROUTINE RT_PARM
 104 !DJG ----------------------------------------------------------------
 105
 106
 107 !DJG ----------------------------------------------------------------
 108 SUBROUTINE GETMAX8DIR(IXX0,JYY0,I,J,H,RETENT_DEP,sox,tmp_gsize,max,XX,YY)
 109     implicit none
 110     INTEGER:: IXX0,JYY0,IXX8,JYY8, XX, YY
 111     INTEGER, INTENT(IN) :: I,J
 112
 113     REAL,INTENT(IN) :: H(XX,YY),RETENT_DEP(XX,YY),sox(XX,YY,8),tmp_gsize(9)
 114     REAL  max
 115     IXX0 = -1
 116     max = 0
 117     if (h(I,J).LE.retent_dep(I,J)) return
 118
 119     IXX8 = I
 120     JYY8 = J+1
 121     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,1),IXX0,JYY0,max,tmp_gsize(1),XX,YY)
 122
 123     IXX8 = I+1
 124     JYY8 = J+1
 125     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,2),IXX0,JYY0,max,tmp_gsize(2),XX,YY)
 126
 127     IXX8 = I+1
 128     JYY8 = J
 129     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,3),IXX0,JYY0,max,tmp_gsize(3),XX,YY)
 130
 131     IXX8 = I+1
 132     JYY8 = J-1
 133     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,4),IXX0,JYY0,max,tmp_gsize(4),XX,YY)
 134
 135     IXX8 = I
 136     JYY8 = J-1
 137     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,5),IXX0,JYY0,max,tmp_gsize(5),XX,YY)
 138
 139     IXX8 = I-1
 140     JYY8 = J-1
 141     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,6),IXX0,JYY0,max,tmp_gsize(6),XX,YY)
 142
 143     IXX8 = I-1
 144     JYY8 = J
 145     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,7),IXX0,JYY0,max,tmp_gsize(7),XX,YY)
 146
 147     IXX8 = I-1
 148     JYY8 = J+1
 149     call GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox(:,:,8),IXX0,JYY0,max,tmp_gsize(8),XX,YY)
 150     RETURN
 151 END SUBROUTINE GETMAX8DIR
 152
 153 SUBROUTINE GET8DIR(IXX8,JYY8,I,J,H,RETENT_DEP,sox   &
 154         ,IXX0,JYY0,max,tmp_gsize,XX,YY)
 155     implicit none
 156     integer,INTENT(INOUT) ::IXX0,JYY0
 157     INTEGER, INTENT(IN) :: I,J,IXX8,JYY8,XX,YY
 158     REAL,INTENT(IN) :: H(XX,YY),RETENT_DEP(XX,YY),sox(XX,YY)
 159     REAL, INTENT(INOUT) ::max
 160     real, INTENT(IN) :: tmp_gsize
 161     real :: sfx
 162
 163     sfx = sox(i,j)-(h(IXX8,JYY8)-h(i,j))*0.001/tmp_gsize
 164     if(sfx .le. 0 ) return
 165     if(max < sfx ) then
 166         IXX0 = IXX8
 167         JYY0 = JYY8
 168         max = sfx
 169     end if
 170
 171 END SUBROUTINE GET8DIR
 172
 173 !DJG------------------------------------------------------------
 174
 175
 176 SUBROUTINE GETSUB8(I, J, XX, YY, wattbl, terrslpNeighbors, distNeighbors, &
 177                    maxneighI, maxneighJ, maxneighIndx, maxneighSlp)
 178    implicit none
 179    integer, intent(in)    :: I, J ! i, j indices for cell of interest
 180    integer, intent(in)    :: XX, YY ! rt domain dimensions
 181    real, intent(in)       :: wattbl(XX, YY) ! water table depth domain array (m)
 182    real, intent(in)       :: terrslpNeighbors(XX, YY, 8)
 183                              ! terrain slope 2d domain array for all 8 neighbors (m/m)
 184    real, intent(in)       :: distNeighbors(9)
 185                              ! lateral distance array for all 8 neighbors of I,J cell (m)
 186    integer, intent(inout) :: maxneighI, maxneighJ ! i, j for max terrain+head slope neighbor
 187    integer, intent(inout) :: maxneighIndx
 188                              ! neighbor cell direction index for max terrain+head slope neighbor
 189    real, intent(inout)    :: maxneighSlp ! terrain+head slope for max slope neighbor (m/m)
 190
 191    ! Local variables
 192    integer                :: neighIndx ! local neighbor direction index (1-8)
 193    integer                :: neighI(8), neighJ(8) ! arrays of neighbor i, j indices
 194
 195    ! Initialize maxslp vars to negative in case max cannot be found
 196    maxneighSlp = -1
 197
 198    ! Setup i, j arrays for neighbors, starting north and rotating clockwise
 199    neighI = (/ I,   I+1, I+1, I+1, I,   I-1, I-1, I-1 /)
 200    neighJ = (/ J+1, J+1, J,   J-1, J-1, J-1, J,   J+1 /)
 201
 202    ! Loop through all neighbors to update max elev+head slope neighbor
 203    do neighIndx = 1, 8
 204       call GETSUB8DIR(I, J, wattbl(I, J), &
 205                       neighI(neighIndx), neighJ(neighIndx), neighIndx, &
 206                       wattbl(neighI(neighIndx), neighJ(neighIndx)), &
 207                       terrslpNeighbors(I,J,neighIndx), distNeighbors(neighIndx), &
 208                       maxneighI, maxneighJ, maxneighIndx, maxneighSlp)
 209    enddo
 210    RETURN
 211 END SUBROUTINE GETSUB8
 212
 213 SUBROUTINE GETSUB8DIR(I, J, selfWattbl, &
 214                       neighI, neighJ, neighIndx, neighWattbl, &
 215                       neighTerrslp, neighDist, &
 216                       maxneighI, maxneighJ, maxneighIndx, maxneighSlp)
 217    implicit none
 218    integer, intent(in)    :: I, J ! i, j indices for cell of interest
 219    real, intent(in)       :: selfWattbl ! water table depth for cell of interest (m)
 220    integer, intent(in)    :: neighI, neighJ ! neighbor cell i, j indices
 221    integer, intent(in)    :: neighIndx ! neighbor cell direction index (1-8)
 222    real, intent(in)       :: neighWattbl ! water table depth for specified neighbor index (m)
 223    real, intent(in)       :: neighTerrslp
 224                              ! terrain slope for specified neighbor index (m/m)
 225    real, intent(in)       :: neighDist
 226                              ! lateral distance for specified neighbor index (m)
 227    integer, intent(inout) :: maxneighI, maxneighJ ! i, j for max terrain+head slope neighbor
 228    integer, intent(inout) :: maxneighIndx
 229                              ! neighbor cell direction index for max terrain+head slope neighbor
 230    real, intent(inout)    :: maxneighSlp
 231                              ! terrain+head slope for max terrain+head slope neighbor (m/m)
 232
 233    ! Local variables
 234    real                   :: dzdx ! subsurface head slope (m/m)
 235    real                   :: beta ! total terrain+head slope (m/m)
 236
 237    ! Calculate total head slope
 238    ! NOTE: wattbl is depth of water table from surface (e.g., 0m if saturated column,
 239    ! 2m if no water table present). So total head slope should be:
 240    !    beta = ( (elev - wattbl) - (elev_neigh - wattbl_neigh) ) / distance
 241    !         = neighTerrslp - dzdx
 242    dzdx = ( selfWattbl - neighWattbl ) / neighDist
 243    beta = neighTerrslp - dzdx
 244    ! Check if this is max and update tracking variables
 245    if ( maxneighSlp < beta ) then
 246       maxneighI = neighI
 247       maxneighJ = neighJ
 248       maxneighSlp = beta
 249       maxneighIndx = neighIndx
 250    end if
 251 END SUBROUTINE GETSUB8DIR
 252
 253
 254 !DJG-----------------------------------------------------------------------
 255 !DJG SUBROUTINE TER_ADJ_SOL    - Terrain adjustment of incoming solar radiation
 256 !DJG-----------------------------------------------------------------------
 257 SUBROUTINE TER_ADJ_SOL(IX,JX,SO8LD_D,TSLP,SHORT,XLAT,XLONG,olddate,DT)
 258
 259 #ifdef MPP_LAND
 260     use module_mpp_land, only:  my_id, io_id, &
 261         mpp_land_bcast_int1
 262 #endif
 263     implicit none
 264     integer,INTENT(IN)     :: IX,JX
 265     INTEGER,INTENT(in), DIMENSION(IX,JX,3)   :: SO8LD_D
 266     real,INTENT(IN), DIMENSION(IX,JX)  :: XLAT,XLONG
 267     real,INTENT(IN) :: DT
 268     real,INTENT(INOUT), DIMENSION(IX,JX)  :: SHORT
 269     character(len=19) :: olddate
 270
 271     ! Local Variables...
 272     real, dimension(IX,JX) ::TSLP,TAZI
 273     real, dimension(IX,JX) ::SOLDN
 274     real :: SOLDEC,DGRD,ITIME2,HRANGLE
 275     real :: BINSH,SOLZANG,SOLAZI,INCADJ
 276     real :: TAZIR,TSLPR,LATR,LONR,SOLDNADJ
 277     integer :: JULDAY0,HHTIME0,MMTIME0,YYYY0,MM0,DD0
 278     integer :: JULDAY,HHTIME,MMTIME,YYYY,MM,DD
 279     integer :: I,J
 280
 281
 282     !----------------------------------------------------------------------
 283     !  SPECIFY PARAMETERS and VARIABLES
 284     !----------------------------------------------------------------------
 285
 286     JULDAY = 0
 287     SOLDN = SHORT
 288     DGRD = 3.14159/180.
 289
 290     ! Set up time variables...
 291 #ifdef MPP_LAND
 292     if(my_id .eq. IO_id) then
 293 #endif
 294         read(olddate(1:4),"(I4)") YYYY0 ! real-time year (GMT)
 295         read(olddate(6:7),"(I2.2)") MM0 ! real-time month (GMT)
 296         read(olddate(9:10),"(I2.2)") DD0 ! real-time day (GMT)
 297         read(olddate(12:13),"(I2.2)") HHTIME0 ! real-time hour (GMT)
 298         read(olddate(15:16),"(I2.2)") MMTIME0 ! real-time minutes (GMT)
 299 #ifdef MPP_LAND
 300     endif
 301     call mpp_land_bcast_int1(YYYY0)
 302     call mpp_land_bcast_int1(MM0)
 303     call mpp_land_bcast_int1(DD0)
 304     call mpp_land_bcast_int1(HHTIME0)
 305     call mpp_land_bcast_int1(MMTIME0)
 306 #endif
 307
 308
 309     ! Set up terrain variables...(returns TSLP&TAZI in radians)
 310     call SLOPE_ASPECT(IX,JX,SO8LD_D,TAZI)
 311
 312     !----------------------------------------------------------------------
 313     !  BEGIN LOOP THROUGH GRID
 314     !----------------------------------------------------------------------
 315     DO J=1,JX
 316         DO I=1,IX
 317             YYYY = YYYY0
 318             MM  = MM0
 319             DD  = DD0
 320             HHTIME = HHTIME0
 321             MMTIME = MMTIME0
 322             call GMT2LOCAL(1,1,XLONG(i,j),YYYY,MM,DD,HHTIME,MMTIME,DT)
 323             call JULDAY_CALC(YYYY,MM,DD,JULDAY)
 324
 325             ! Convert to radians...
 326             LATR = XLAT(I,J)   !send solsub local lat in deg
 327             LONR = XLONG(I,J)   !send solsub local lon in deg
 328             TSLPR = TSLP(I,J)/DGRD !send solsub local slp in deg
 329             TAZIR = TAZI(I,J)/DGRD !send solsub local azim in deg
 330
 331             !Call SOLSUB to return terrain adjusted incoming solar radiation...
 332             ! SOLSUB taken from Whiteman and Allwine, 1986, Environ. Software.
 333
 334             call SOLSUB(LONR,LATR,TAZIR,TSLPR,SOLDN(I,J),YYYY,MM,         &
 335                 DD,HHTIME,MMTIME,SOLDNADJ,SOLZANG,SOLAZI,INCADJ)
 336
 337             SOLDN(I,J)=SOLDNADJ
 338
 339         ENDDO
 340     ENDDO
 341
 342     SHORT = SOLDN
 343
 344     return
 345 end SUBROUTINE TER_ADJ_SOL
 346 !DJG-----------------------------------------------------------------------
 347 !DJG END SUBROUTINE TER_ADJ_SOL
 348 !DJG-----------------------------------------------------------------------
 349
 350
 351 !DJG-----------------------------------------------------------------------
 352 !DJG SUBROUTINE GMT2LOCAL
 353 !DJG-----------------------------------------------------------------------
 354 subroutine GMT2LOCAL(IX,JX,XLONG,YY,MM,DD,HH,MIN,DT)
 355
 356     implicit none
 357
 358     !!! Declare Passed Args.
 359
 360     INTEGER, INTENT(INOUT) :: yy,mm,dd,hh,min
 361     INTEGER, INTENT(IN) :: IX,JX
 362     REAL,INTENT(IN), DIMENSION(IX,JX)  :: XLONG
 363     REAL,INTENT(IN) :: DT
 364
 365     !!! Declare local variables
 366
 367     integer :: i,j,minflag,hhflag,ddflag,mmflag,yyflag
 368     integer :: adj_min,lst_adj_min,lst_adj_hh,adj_hh
 369     real, dimension(IX,JX) :: TDIFF
 370     real :: tmp
 371     integer :: yyinit,mminit,ddinit,hhinit,mininit
 372
 373     !!! Initialize flags
 374     hhflag=0
 375     ddflag=0
 376     mmflag=0
 377     yyflag=0
 378
 379     !!! Set up constants...
 380     yyinit = yy
 381     mminit = mm
 382     ddinit = dd
 383     hhinit = hh
 384     mininit = min
 385
 386
 387     ! Loop through data...
 388     do j=1,JX
 389         do i=1,IX
 390
 391             ! Reset yy,mm,dd...
 392             yy = yyinit
 393             mm = mminit
 394             dd = ddinit
 395             hh = hhinit
 396             min = mininit
 397
 398             !!! Set up adjustments...
 399             !   - assumes +E , -W  longitude and 0.06667 hr/deg (=24/360)
 400             TDIFF(I,J) = XLONG(I,J)*0.06667   ! time offset in hr
 401             tmp = TDIFF(I,J)
 402             lst_adj_hh = INT(tmp)
 403             lst_adj_min = NINT(MOD(int(tmp),1)*60.) + int(DT/2./60.)  ! w/ 1/2 timestep adjustment...
 404
 405             !!! Process Minutes...
 406             adj_min = min+lst_adj_min
 407             if (adj_min.lt.0) then
 408                 min=60+adj_min
 409                 lst_adj_hh = lst_adj_hh - 1
 410             else if (adj_min.ge.0.AND.adj_min.lt.60) then
 411                 min=adj_min
 412             else if (adj_min.ge.60) then
 413                 min=adj_min-60
 414                 lst_adj_hh = lst_adj_hh + 1
 415             end if
 416
 417             !!! Process Hours
 418             adj_hh = hh+lst_adj_hh
 419             if (adj_hh.lt.0) then
 420                 hh = 24+adj_hh
 421                 ddflag=1
 422             else if (adj_hh.ge.0.AND.adj_hh.lt.24) then
 423                 hh=adj_hh
 424             else if (adj_hh.ge.24) then
 425                 hh=adj_hh-24
 426                 ddflag = 2
 427             end if
 428
 429
 430
 431             !!! Process Days, Months, Years
 432             ! Subtract a day
 433             if (ddflag.eq.1) then
 434                 if (dd.gt.1) then
 435                     dd=dd-1
 436                 else
 437                     if (mm.eq.1) then
 438                         mm=12
 439                         yy=yy-1
 440                     else
 441                         mm=mm-1
 442                     end if
 443                     if ((mm.eq.1).or.(mm.eq.3).or.(mm.eq.5).or. &
 444                             (mm.eq.7).or.(mm.eq.8).or.(mm.eq.10).or. &
 445                             (mm.eq.12)) then
 446                         dd=31
 447                     else
 448
 449                         !!! Adjustment for leap years!!!
 450                         if(mm.eq.2) then
 451                             if(MOD(yy,4).eq.0) then
 452                                 dd=29
 453                             else
 454                                 dd=28
 455                             end if
 456                         end if
 457                         if(mm.ne.2) dd=30
 458                     end if
 459                 end if
 460             end if
 461
 462             ! Add a day
 463             if (ddflag.eq.2) then
 464                 if ((mm.eq.1).or.(mm.eq.3).or.(mm.eq.5).or. &
 465                         (mm.eq.7).or.(mm.eq.8).or.(mm.eq.10).or. &
 466                         (mm.eq.12)) then
 467                     if (dd.eq.31) then
 468                         dd=1
 469                         if (mm.eq.12) then
 470                             mm=1
 471                             yy=yy+1
 472                         else
 473                             mm=mm+1
 474                         end if
 475                     else
 476                         dd=dd+1
 477                     end if
 478
 479                     !!! Adjustment for leap years!!!
 480                 else if (mm.eq.2) then
 481                     if(MOD(yy,4).eq.0) then
 482                         if (dd.eq.29) then
 483                             dd=1
 484                             mm=3
 485                         else
 486                             dd=dd+1
 487                         end if
 488                     else
 489                         if (dd.eq.28) then
 490                             dd=1
 491                             mm=3
 492                         else
 493                             dd=dd+1
 494                         end if
 495                     end if
 496                 else
 497                     if (dd.eq.30) then
 498                         dd=1
 499                         mm=mm+1
 500                     else
 501                         dd=dd+1
 502                     end if
 503                 end if
 504
 505             end if
 506
 507         end do   !i-loop
 508     end do   !j-loop
 509
 510     return
 511 end subroutine
 512
 513 !DJG-----------------------------------------------------------------------
 514 !DJG END SUBROUTINE GMT2LOCAL
 515 !DJG-----------------------------------------------------------------------
 516
 517
 518
 519 !DJG-----------------------------------------------------------------------
 520 !DJG SUBROUTINE JULDAY_CALC
 521 !DJG-----------------------------------------------------------------------
 522 subroutine JULDAY_CALC(YYYY,MM,DD,JULDAY)
 523
 524     implicit none
 525     integer,intent(in) :: YYYY,MM,DD
 526     integer,intent(out) :: JULDAY
 527
 528     integer :: resid
 529     integer julm(13)
 530     DATA JULM/0, 31, 59, 90, 120, 151, 181, 212, 243, 273, &
 531         304, 334, 365 /
 532
 533     integer LPjulm(13)
 534     DATA LPJULM/0, 31, 60, 91, 121, 152, 182, 213, 244, 274, &
 535         305, 335, 366 /
 536
 537     resid = MOD(YYYY,4) !Set up leap year check...
 538
 539     if (resid.ne.0) then    !If not a leap year....
 540         JULDAY = JULM(MM) + DD
 541     else                    !If a leap year...
 542         JULDAY = LPJULM(MM) + DD
 543     end if
 544
 545     RETURN
 546 END subroutine JULDAY_CALC
 547 !DJG-----------------------------------------------------------------------
 548 !DJG END SUBROUTINE JULDAY
 549 !DJG-----------------------------------------------------------------------
 550
 551 !DJG-----------------------------------------------------------------------
 552 !DJG SUBROUTINE SLOPE_ASPECT
 553 !DJG-----------------------------------------------------------------------
 554 subroutine SLOPE_ASPECT(IX,JX,SO8LD_D,TAZI)
 555
 556     implicit none
 557     integer, INTENT(IN)            :: IX,JX
 558     !   real,INTENT(in),DIMENSION(IX,JX)   :: TSLP  !terrain slope (m/m)
 559     real,INTENT(OUT),DIMENSION(IX,JX)   :: TAZI  !terrain aspect (deg)
 560
 561     INTEGER, DIMENSION(IX,JX,3)   :: SO8LD_D
 562     real :: DGRD
 563     integer :: i,j
 564
 565     !   TSLP = 0.  !Initialize as flat
 566     TAZI = 0.  !Initialize as north facing
 567
 568     ! Find steepest descent slope and direction...
 569     do j=1,JX
 570         do i=1,IX
 571             !   TSLP(I,J) = TANH(Vmax(i,j)) ! calculate slope in radians...
 572
 573             ! Convert steepest slope and aspect to radians...
 574             IF (SO8LD_D(i,j,3).eq.1) then
 575                 TAZI(I,J) = 0.0
 576             ELSEIF (SO8LD_D(i,j,3).eq.2) then
 577                 TAZI(I,J) = 45.0
 578             ELSEIF (SO8LD_D(i,j,3).eq.3) then
 579                 TAZI(I,J) = 90.0
 580             ELSEIF (SO8LD_D(i,j,3).eq.4) then
 581                 TAZI(I,J) = 135.0
 582             ELSEIF (SO8LD_D(i,j,3).eq.5) then
 583                 TAZI(I,J) = 180.0
 584             ELSEIF (SO8LD_D(i,j,3).eq.6) then
 585                 TAZI(I,J) = 225.0
 586             ELSEIF (SO8LD_D(i,j,3).eq.7) then
 587                 TAZI(I,J) = 270.0
 588             ELSEIF (SO8LD_D(i,j,3).eq.8) then
 589                 TAZI(I,J) = 315.0
 590             END IF
 591
 592             DGRD = 3.141593/180.
 593             TAZI(I,J) = TAZI(I,J)*DGRD ! convert azimuth to radians...
 594
 595         END DO
 596     END DO
 597
 598     RETURN
 599 END  subroutine SLOPE_ASPECT
 600 !DJG-----------------------------------------------------------------------
 601 !DJG END SUBROUTINE SLOPE_ASPECT
 602 !DJG-----------------------------------------------------------------------
 603
 604 !DJG----------------------------------------------------------------
 605 !DJG    SUBROUTINE SOLSUB
 606 !DJG----------------------------------------------------------------
 607 SUBROUTINE SOLSUB(LONG,LAT,AZ,IN,SC,YY,MO,IDA,IHR,MM,OUT1, &
 608         OUT2,OUT3,INCADJ)
 609
 610
 611     ! Notes....
 612
 613     implicit none
 614     logical               :: daily, first
 615     integer               :: yy,mo,ida,ihr,mm,d
 616     integer,dimension(12) :: nday
 617     real                  :: lat,long,longcor,longsun,in,inslo
 618     real :: az,sc,out1,out2,out3,cosbeta,dzero,eccent,pi,calint
 619     real :: rtod,decmax,omega,onehr,omd,omdzero,rdvecsq,sdec
 620     real :: declin,cdec,arg,declon,sr,stdmrdn,b,em,timnoon,azslo
 621     real :: slat,clat,caz,saz,sinc,cinc,hinc,h,cosz,extra,extslo
 622     real :: t1,z,cosa,a,cosbeta_flat,INCADJ
 623     integer :: HHTIME,MMTIME,i,ik
 624     real, dimension(4) :: ACOF,BCOF
 625
 626     ! Constants
 627     daily=.FALSE.
 628     ACOF(1) = 0.00839
 629     ACOF(2) = -0.05391
 630     ACOF(3) = -0.00154
 631     ACOF(4) = -0.0022
 632     BCOF(1) = -0.12193
 633     BCOF(2) = -0.15699
 634     BCOF(3) = -0.00657
 635     BCOF(4) = -0.00370
 636     DZERO = 80.
 637     ECCENT = 0.0167
 638     PI = 3.14159
 639     CALINT = 1.
 640     RTOD = PI / 180.
 641     DECMAX=(23.+26./60.)*RTOD
 642     OMEGA=2*PI/365.
 643     ONEHR=15.*RTOD
 644
 645     ! Calculate Julian Day...
 646     D = 0
 647     call JULDAY_CALC(YY,MO,IDA,D)
 648
 649     ! Ratio of radius vectors squared...
 650     OMD=OMEGA*D
 651     OMDZERO=OMEGA*DZERO
 652     !       RDVECSQ=1./(1.-ECCENT*COS(OMD))**2
 653     RDVECSQ = 1.    ! no adjustment for orbital changes when coupled to HRLDAS...
 654
 655     ! Declination of sun...
 656     LONGSUN=OMEGA*(D-Dzero)+2.*ECCENT*(SIN(OMD)-SIN(OMDZERO))
 657     DECLIN=ASIN(SIN(DECMAX)*SIN(LONGSUN))
 658     SDEC=SIN(DECLIN)
 659     CDEC=COS(DECLIN)
 660
 661     ! Check for Polar Day/night...
 662     ARG=((PI/2.)-ABS(DECLIN))/RTOD
 663     IF(ABS(LAT).GT.ARG) THEN
 664         IF((LAT.GT.0..AND.DECLIN.LT.0) .OR.       &
 665                 (LAT.LT.0..AND.DECLON.GT.0.)) THEN
 666             OUT1 = 0.
 667             OUT2 = 0.
 668             OUT3 = 0.
 669             RETURN
 670         ENDIF
 671         SR=-1.*PI
 672     ELSE
 673
 674         ! Calculate sunrise hour angle...
 675         SR=-1.*ABS(ACOS(-1.*TAN(LAT*RTOD)*TAN(DECLIN)))
 676     END IF
 677
 678     ! Find standard meridian for site
 679     STDMRDN=NINT(LONG/15.)*15.
 680     LONGCOR=(LONG-STDMRDN)/15.
 681
 682     ! Compute time correction from equation of time...
 683     B=2.*PI*(D-.4)/365
 684     EM=0.
 685     DO I=1,4
 686         EM=EM+(BCOF(I)*SIN(I*B)+ACOF(I)*COS(I*B))
 687     END DO
 688
 689     ! Compute time of solar noon...
 690     TIMNOON=12.-EM-LONGCOR
 691
 692     ! Set up a few more terms...
 693     AZSLO=AZ*RTOD
 694     INSLO=IN*RTOD
 695     SLAT=SIN(LAT*RTOD)
 696     CLAT=COS(LAT*RTOD)
 697     CAZ=COS(AZSLO)
 698     SAZ=SIN(AZSLO)
 699     SINC=SIN(INSLO)
 700     CINC=COS(INSLO)
 701
 702     ! Begin solar radiation calculations...daily first, else instantaneous...
 703     IF (DAILY) THEN   ! compute daily integrated values...(Not used in HRLDAS!)
 704         IHR=0
 705         MM=0
 706         HINC=CALINT*ONEHR/60.
 707         IK=(2.*ABS(SR)/HINC)+2.
 708         FIRST=.TRUE.
 709         OUT1=0.
 710         DO I=1,IK
 711             H=SR+HINC*FLOAT(I-1)
 712             COSZ=SLAT*SDEC+CLAT*CDEC*COS(H)
 713             COSBETA=CDEC*((SLAT*COS(H))*(-1.*CAZ*SINC)- &
 714                 SIN(H)*(SAZ*SINC)+(CLAT*COS(H))*CINC)+ &
 715                 SDEC*(CLAT*(CAZ*SINC)+SLAT*CINC)
 716             EXTRA=SC*RDVECSQ*COSZ
 717             IF(EXTRA.LE.0.) EXTRA=0.
 718             EXTSLO=SC*RDVECSQ*COSBETA
 719             IF(EXTRA.LE.0. .OR. EXTSLO.LT.0.) EXTSLO=0.
 720             IF(FIRST .AND. EXTSLO.GT.0.) THEN
 721                 OUT2=(H-HINC)/ONEHR+TIMNOON
 722                 FIRST = .FALSE.
 723             END IF
 724             IF(.NOT.FIRST .AND. EXTSLO.LE.0.) OUT3=H/ONEHR+TIMNOON
 725             OUT1=EXTSLO+OUT1
 726         END DO
 727         OUT1=OUT1*CALINT*60./1000000.
 728
 729     ELSE   ! Compute instantaneous values...(Is used in HRLDAS!)
 730
 731         T1=FLOAT(IHR)+FLOAT(MM)/60.
 732         H=ONEHR*(T1-TIMNOON)
 733         COSZ=SLAT*SDEC+CLAT*CDEC*COS(H)
 734
 735         ! Assuming HRLDAS forcing already accounts for season, time of day etc,
 736         ! subtract out the component of adjustment that would occur for
 737         ! a flat surface, this should leave only the sloped component remaining
 738
 739         COSBETA=CDEC*((SLAT*COS(H))*(-1.*CAZ*SINC)-  &
 740             SIN(H)*(SAZ*SINC)+(CLAT*COS(H))*CINC)+ &
 741             SDEC*(CLAT*(CAZ*SINC)+SLAT*CINC)
 742
 743         COSBETA_FLAT=CDEC*CLAT*COS(H)+SDEC*SLAT
 744
 745         INCADJ = COSBETA+(1-COSBETA_FLAT)
 746
 747         EXTRA=SC*RDVECSQ*COSZ
 748         IF(EXTRA.LE.0.) EXTRA=0.
 749         EXTSLO=SC*RDVECSQ*INCADJ
 750         !         IF(EXTRA.LE.0. .OR. EXTSLO.LT.0.) EXTSLO=0.  !remove check for HRLDAS.
 751         OUT1=EXTSLO
 752         Z=ACOS(COSZ)
 753         COSA=(SLAT*COSZ-SDEC)/(CLAT*SIN(Z))
 754         IF(COSA.LT.-1.) COSA=-1.
 755         IF(COSA.GT.1.) COSA=1.
 756         A=ABS(ACOS(COSA))
 757         IF(H.LT.0.) A=-A
 758         OUT2=Z/RTOD
 759         OUT3=A/RTOD+180
 760
 761     END IF    ! End if for daily vs instantaneous values...
 762
 763     !DJG-----------------------------------------------------------------------
 764     RETURN
 765 END SUBROUTINE SOLSUB
 766 !DJG-----------------------------------------------------------------------
 767
 768 subroutine seq_land_SO8(SO8LD_D,Vmax,TERR,dx,ix,jx)
 769     implicit none
 770     integer :: ix,jx,i,j
 771     REAL, DIMENSION(IX,JX,8)      :: SO8LD
 772     INTEGER, DIMENSION(IX,JX,3)   :: SO8LD_D
 773     real,DIMENSION(IX,JX)      :: TERR
 774     real                       :: dx(ix,jx,9),Vmax(ix,jx)
 775     SO8LD_D = -1
 776     do j = 2, jx -1
 777         do i = 2, ix -1
 778             SO8LD(i,j,1) = (TERR(i,j)-TERR(i,j+1))/dx(i,j,1)
 779             SO8LD_D(i,j,1) = i
 780             SO8LD_D(i,j,2) = j + 1
 781             SO8LD_D(i,j,3) = 1
 782             Vmax(i,j) = SO8LD(i,j,1)
 783
 784             SO8LD(i,j,2) = (TERR(i,j)-TERR(i+1,j+1))/DX(i,j,2)
 785             if(SO8LD(i,j,2) .gt. Vmax(i,j) ) then
 786                 SO8LD_D(i,j,1) = i + 1
 787                 SO8LD_D(i,j,2) = j + 1
 788                 SO8LD_D(i,j,3) = 2
 789                 Vmax(i,j) = SO8LD(i,j,2)
 790             end if
 791             SO8LD(i,j,3) = (TERR(i,j)-TERR(i+1,j))/DX(i,j,3)
 792             if(SO8LD(i,j,3) .gt. Vmax(i,j) ) then
 793                 SO8LD_D(i,j,1) = i + 1
 794                 SO8LD_D(i,j,2) = j
 795                 SO8LD_D(i,j,3) = 3
 796                 Vmax(i,j) = SO8LD(i,j,3)
 797             end if
 798             SO8LD(i,j,4) = (TERR(i,j)-TERR(i+1,j-1))/DX(i,j,4)
 799             if(SO8LD(i,j,4) .gt. Vmax(i,j) ) then
 800                 SO8LD_D(i,j,1) = i + 1
 801                 SO8LD_D(i,j,2) = j - 1
 802                 SO8LD_D(i,j,3) = 4
 803                 Vmax(i,j) = SO8LD(i,j,4)
 804             end if
 805             SO8LD(i,j,5) = (TERR(i,j)-TERR(i,j-1))/DX(i,j,5)
 806             if(SO8LD(i,j,5) .gt. Vmax(i,j) ) then
 807                 SO8LD_D(i,j,1) = i
 808                 SO8LD_D(i,j,2) = j - 1
 809                 SO8LD_D(i,j,3) = 5
 810                 Vmax(i,j) = SO8LD(i,j,5)
 811             end if
 812             SO8LD(i,j,6) = (TERR(i,j)-TERR(i-1,j-1))/DX(i,j,6)
 813             if(SO8LD(i,j,6) .gt. Vmax(i,j) ) then
 814                 SO8LD_D(i,j,1) = i - 1
 815                 SO8LD_D(i,j,2) = j - 1
 816                 SO8LD_D(i,j,3) = 6
 817                 Vmax(i,j) = SO8LD(i,j,6)
 818             end if
 819             SO8LD(i,j,7) = (TERR(i,j)-TERR(i-1,j))/DX(i,j,7)
 820             if(SO8LD(i,j,7) .gt. Vmax(i,j) ) then
 821                 SO8LD_D(i,j,1) = i - 1
 822                 SO8LD_D(i,j,2) = j
 823                 SO8LD_D(i,j,3) = 7
 824                 Vmax(i,j) = SO8LD(i,j,7)
 825             end if
 826             SO8LD(i,j,8) = (TERR(i,j)-TERR(i-1,j+1))/DX(i,j,8)
 827             if(SO8LD(i,j,8) .gt. Vmax(i,j) ) then
 828                 SO8LD_D(i,j,1) = i - 1
 829                 SO8LD_D(i,j,2) = j + 1
 830                 SO8LD_D(i,j,3) = 8
 831                 Vmax(i,j) = SO8LD(i,j,8)
 832             end if
 833         enddo
 834     enddo
 835     Vmax = TANH(Vmax)
 836     return
 837 end  subroutine seq_land_SO8
 838
 839 #ifdef MPP_LAND
 840 subroutine MPP_seq_land_SO8(SO8LD_D,Vmax,TERRAIN,dx,ix,jx,&
 841         global_nx,global_ny)
 842
 843     use module_mpp_land, only:  my_id, io_id, &
 844         write_io_real,decompose_data_int,decompose_data_real
 845
 846     implicit none
 847     integer,intent(in) :: ix,jx,global_nx,global_ny
 848     INTEGER, intent(inout),DIMENSION(IX,JX,3)   :: SO8LD_D
 849     !         real,intent(in), DIMENSION(IX,JX)   :: TERRAIN
 850     real,DIMENSION(IX,JX)   :: TERRAIN
 851     real,intent(out),dimension(ix,jx) ::  Vmax
 852     real,intent(in)                     :: dx(ix,jx,9)
 853     real                     :: g_dx(ix,jx,9)
 854
 855     real,DIMENSION(global_nx,global_ny)      :: g_TERRAIN
 856     real,DIMENSION(global_nx,global_ny)      :: g_Vmax
 857     integer,DIMENSION(global_nx,global_ny,3)      :: g_SO8LD_D
 858     integer :: k
 859
 860     g_SO8LD_D = 0
 861     g_Vmax    = 0
 862
 863     do k = 1, 9
 864         call write_IO_real(dx(:,:,k),g_dx(:,:,k))
 865     end do
 866
 867     call write_IO_real(TERRAIN,g_TERRAIN)
 868     if(my_id .eq. IO_id) then
 869         call seq_land_SO8(g_SO8LD_D,g_Vmax,g_TERRAIN,g_dx,global_nx,global_ny)
 870     endif
 871     call decompose_data_int(g_SO8LD_D(:,:,3),SO8LD_D(:,:,3))
 872     call decompose_data_real(g_Vmax,Vmax)
 873     return
 874 end subroutine MPP_seq_land_SO8
 875
 876 #endif
 877
 878
 879 subroutine disaggregateDomain_drv(did)
 880    use module_RT_data, only: rt_domain
 881    use config_base, only: nlst, noah_lsm
 882    integer :: did
 883    call disaggregateDomain( RT_DOMAIN(did)%IX, RT_DOMAIN(did)%JX, nlst(did)%NSOIL, &
 884                             RT_DOMAIN(did)%IXRT, RT_DOMAIN(did)%JXRT, nlst(did)%AGGFACTRT, &
 885                             RT_DOMAIN(did)%SICE, RT_DOMAIN(did)%SMC, RT_DOMAIN(did)%SH2OX, &
 886                             RT_DOMAIN(did)%INFXSRT, rt_domain(did)%dist_lsm(:,:,9), &
 887                             RT_DOMAIN(did)%SMCMAX1, RT_DOMAIN(did)%SMCREF1, &
 888                             RT_DOMAIN(did)%SMCWLT1, RT_DOMAIN(did)%VEGTYP, RT_DOMAIN(did)%LKSAT, &
 889                             RT_DOMAIN(did)%NEXP, &
 890                             rt_domain(did)%overland%properties%distance_to_neighbor, &
 891                             RT_DOMAIN(did)%INFXSWGT, &
 892                             RT_DOMAIN(did)%LKSATFAC, &
 893                             rt_domain(did)%overland%streams_and_lakes%ch_netrt, &
 894                             RT_DOMAIN(did)%SH2OWGT, &
 895                             RT_DOMAIN(did)%subsurface%grid_transform%smcrefrt, &
 896                             rt_domain(did)%overland%control%infiltration_excess, &
 897                             RT_DOMAIN(did)%subsurface%grid_transform%smcmaxrt, &
 898                             RT_DOMAIN(did)%subsurface%grid_transform%smcwltrt, &
 899                             rt_domain(did)%subsurface%grid_transform%smcrt, &
 900                             rt_domain(did)%overland%streams_and_lakes%lake_mask, &
 901                             rt_domain(did)%subsurface%properties%lksatrt, &
 902                             rt_domain(did)%subsurface%properties%nexprt, &
 903                             RT_DOMAIN(did)%subsurface%properties%sldpth, &
 904                             RT_DOMAIN(did)%soiltypRT, RT_DOMAIN(did)%soiltyp, &
 905                             rt_domain(did)%ELRT, RT_DOMAIN(did)%iswater, &
 906                             rt_domain(did)%IMPERVFRAC, nlst(did)%imperv_adj)
 907 end subroutine disaggregateDomain_drv
 908
 909 !===================================================================================================
 910 ! Subroutine Name:
 911 !   subroutine disaggregateDomain
 912 ! Author(s)/Contact(s):
 913 !   D. Gochis and W. Yu <gochis><ucar><edu>
 914 ! Abstract:
 915 !   Disaggregates states and parameters from coarse grid to fine grid
 916 ! History Log:
 917 ! Usage:
 918 ! Parameters:
 919 ! Input Files:
 920 ! Output Files:
 921 ! Condition codes:
 922 ! User controllable options: None.
 923 ! Notes:
 924 !===================================================================================================
 925 subroutine disaggregateDomain(IX, JX, NSOIL, IXRT, JXRT, AGGFACTRT, &
 926                               SICE, SMC, SH2OX, INFXSRT, area_lsm, SMCMAX1, SMCREF1, &
 927                               SMCWLT1, VEGTYP, LKSAT, NEXP, dist, INFXSWGT, &
 928                               LKSATFAC, CH_NETRT, SH2OWGT, SMCREFRT, INFXSUBRT, SMCMAXRT, &
 929                               SMCWLTRT, SMCRT, LAKE_MSKRT, LKSATRT, NEXPRT,  &
 930                               SLDPTH, soiltypRT, soiltyp, elrt, iswater, impervfrac, imperv_adj)
 931 #ifdef MPP_LAND
 932    use module_mpp_land, only: left_id,down_id,right_id, &
 933                               up_id,mpp_land_com_real, my_id, io_id, numprocs, &
 934                               mpp_land_sync,mpp_land_com_integer,mpp_land_max_int1, &
 935                               sum_real1
 936    use module_hydro_stop, only:HYDRO_stop
 937 #endif
 938
 939    implicit none
 940
 941    ! Input Variables ------------------------------------------------------------------------
 942    ! General geometry/domain info:
 943    integer, intent(in)                           :: IX,JX ! coarse grid i,j dims
 944    integer, intent(in)                           :: IXRT,JXRT ! fine grid i,j dims
 945    integer, intent(in)                           :: AGGFACTRT ! aggregation factor between grids
 946    integer, intent(in)                           :: NSOIL ! number of soil layers
 947    integer, intent(in)                           :: iswater ! water veg class (from geogrid attrib)
 948    real, intent(in), dimension(NSOIL)            :: SLDPTH ! array soil layer depth intervals (m)
 949    integer, intent(in)                           :: imperv_adj ! impervious configuration option from hydro.namelist
 950    ! LSM grid parameters:
 951    real, intent(in),  dimension(IX,JX)           :: area_lsm ! cell area on the coarse grid (m2)
 952    integer, intent(in), dimension(IX,JX)         :: VEGTYP, soiltyp ! coarse grid veg and soil types
 953    real, intent(in),  dimension(IX,JX)           :: SMCMAX1 ! coarse grid porosity
 954    real, intent(in),  dimension(IX,JX)           :: SMCREF1 ! coarse grid field capacity
 955    real, intent(in),  dimension(IX,JX)           :: SMCWLT1 ! coarse grid wilting point
 956    real, intent(in),  dimension(IX,JX)           :: LKSAT ! coarse grid lateral ksat (m/s)
 957    real, intent(in),  dimension(IX,JX)           :: NEXP ! coarse grid n exponent
 958    ! LSM states:
 959    real, intent(in),  dimension(IX,JX,NSOIL)     :: SMC ! total soil moisture (m3/m3)
 960    real, intent(in),  dimension(IX,JX,NSOIL)     :: SH2OX ! liquid soil moisture (m3/m3)
 961    real, intent(in),  dimension(IX,JX)           :: INFXSRT
 962                                                     ! infiltration excess on coarse grid (mm)
 963    ! Routing grid parameters:
 964    real, intent(in), dimension(IXRT,JXRT,9)      :: dist
 965                                                     ! routing grid cell distances (m) and area (m2)
 966                                                     ! TODO: can we just pass in area since we don't need other distances?
 967    real, intent(in), dimension(IXRT,JXRT)        :: LKSATFAC ! lateral ksat adj factor
 968    real, intent(in), dimension(IXRT,JXRT)        :: elrt ! elevation grid (m)
 969    integer, intent(in), dimension(IXRT,JXRT)     :: CH_NETRT ! channel network routing grid
 970    real, intent(in), dimension(IXRT,JXRT)        :: impervfrac ! impervious fraction
 971    ! Routing states:
 972    real, intent(in), dimension(IXRT,JXRT)        :: INFXSWGT ! infiltration excess weighting grid
 973    real, intent(in), dimension(IXRT,JXRT,NSOIL)  :: SH2OWGT ! soil moisture weighting grid
 974
 975    ! Update Variables ------------------------------------------------------------------------
 976    integer, intent(inout), dimension(IXRT,JXRT)  :: LAKE_MSKRT ! lake mask on the routing grid
 977
 978    ! Output Variables ------------------------------------------------------------------------
 979    ! Parameters:
 980    integer, intent(out), dimension(IXRT,JXRT)    :: soiltypRT ! soil type on the routing grid
 981    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCMAXRT ! porosity on routing grid
 982    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCREFRT ! field capacity on routing grid
 983    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCWLTRT ! wilting point on routing grid
 984    real, intent(out), dimension(IXRT,JXRT)       :: LKSATRT ! lateral ksat on the routing grid (m/s)
 985    real, intent(out), dimension(IXRT,JXRT)       :: NEXPRT ! n exponent on the routing grid
 986    ! States:
 987    real, intent(out), dimension(IX,JX,NSOIL)     :: SICE ! soil ice content on coarse grid (m3/m3)
 988    real, intent(out), dimension(IXRT,JXRT,NSOIL) :: SMCRT
 989                                                     ! soil moisture contant on routing grid (m3/m3)
 990    real, intent(out), dimension(IXRT,JXRT)       :: INFXSUBRT
 991                                                     ! infiltration excess on routing grid (mm)
 992
 993    ! Local Variables ------------------------------------------------------------------------
 994    integer                    :: i, j ! coarse grid loop indices
 995    integer                    :: AGGFACYRT, AGGFACXRT ! fine grid aggregation factors
 996    integer                    :: IXXRT, JYYRT ! fine grid i,j coordinates
 997    integer                    :: KRT, KF ! soil layer loop indixes
 998    real                       :: LSMVOL ! total infiltration excess volume per LSM cell (mm * m2)
 999    real                       :: SMCEXCS ! excess soil moisture (m3/m3)
1000    real                       :: WATHOLDCAP ! water holding capacity, smcmax - smcwlt (m3/m3)
1001    real                       :: smScaleFact ! soil moisture scaling factor for ksat (0-1)
1002    real, dimension(IXRT,JXRT) :: OCEAN_INFXSUBRT
1003                                  ! dummy variable to dump infiltration excess over ocean cells
1004 #ifdef HYDRO_D
1005    ! For water budget calcs
1006    integer :: ii, jj, kk ! local loop indices
1007    real    :: smctot1,smcrttot2 ! total soil moisture, start and finish
1008    real    :: sicetot1 ! total ice, start and finish
1009    real    :: suminfxs1,suminfxsrt2 ! total infiltration excess, start and finish
1010 #endif
1011
1012    ! Initialize variables
1013    SICE = SMC - SH2OX
1014    SMCREFRT = 0.0
1015    ! ADCHANGE: Initialize ocean infxsubrt var to 0. Currently just a dump
1016    ! variable but could be used for future ocean model coupling
1017    OCEAN_INFXSUBRT = 0.0
1018
1019    !DJG First, Disaggregate a few key fields for routing...
1020    !DJG Debug...
1021 #ifdef HYDRO_D
1022    print *, "Beginning Disaggregation..."
1023 #endif
1024
1025    !DJG Mass balance check for disagg...
1026 #ifdef HYDRO_D
1027    ! ADCHANGE: START Initial water balance variables
1028    ! ALL VARS in MM
1029    suminfxs1 = 0.
1030    smctot1 = 0.
1031    sicetot1 = 0.
1032    do ii=1,IX
1033       do jj=1,JX
1034          suminfxs1 = suminfxs1 + INFXSRT(ii,jj) / float(IX*JX)
1035          do kk=1,NSOIL
1036             smctot1 = smctot1 + SMC(ii,jj,kk)*SLDPTH(kk)*1000. / float(IX*JX)
1037             sicetot1 = sicetot1 + SICE(ii,jj,kk)*SLDPTH(kk)*1000. / float(IX*JX)
1038          end do
1039       end do
1040    end do
1041 #ifdef MPP_LAND
1042    ! not tested
1043    CALL sum_real1(suminfxs1)
1044    CALL sum_real1(smctot1)
1045    CALL sum_real1(sicetot1)
1046    suminfxs1 = suminfxs1/float(numprocs)
1047    smctot1 = smctot1/float(numprocs)
1048    sicetot1 = sicetot1/float(numprocs)
1049 #endif
1050    ! END Initial water balance variables
1051 #endif
1052
1053    !DJG Weighting alg. alteration...(prescribe wghts if time = 1)
1054
1055    do J=1,JX ! Start coarse grid j loop
1056       do I=1,IX ! Start coarse grid i loop
1057
1058          !DJG Weighting alg. alteration...
1059          LSMVOL = INFXSRT(I,J) * area_lsm(I,J) ! mm * m2
1060
1061          do AGGFACYRT=AGGFACTRT-1,0,-1 ! Start disagg fine grid j loop
1062             do AGGFACXRT=AGGFACTRT-1,0,-1 ! Start disagg fine grid i loop
1063
1064                IXXRT = I * AGGFACTRT - AGGFACXRT ! Define fine grid i
1065                JYYRT = J * AGGFACTRT - AGGFACYRT ! Define fine grid j
1066 #ifdef MPP_LAND
1067                if(left_id.ge.0) IXXRT = IXXRT+1
1068                if(down_id.ge.0) JYYRT = JYYRT+1
1069 #endif
1070                ! Initial redistribution of total coarse grid cell infiltration excess
1071                ! based on subgrid weight factor and current volume of water (mm * m2 / m2 = mm)
1072                INFXSUBRT(IXXRT,JYYRT) = LSMVOL *     &
1073                                         INFXSWGT(IXXRT,JYYRT) / dist(IXXRT,JYYRT,9)
1074
1075                do KRT=1,NSOIL ! Soil layer loop
1076
1077                   ! Adjustments for soil ice
1078                   IF (SICE(I,J,KRT) .gt. 0) then
1079                      !DJG Adjust SMCMAX for SICE when subsfc routing...make 3d variable
1080                      SMCMAXRT(IXXRT,JYYRT,KRT) = SMCMAX1(I,J) - SICE(I,J,KRT)
1081                      SMCREFRT(IXXRT,JYYRT,KRT) = SMCREF1(I,J) - SICE(I,J,KRT) !TODO: This can be negative! e.g., when almost saturated and fully frozen
1082                      WATHOLDCAP = SMCMAX1(I,J) - SMCWLT1(I,J)
1083                      IF (SICE(I,J,KRT) .le. WATHOLDCAP)    then
1084                         SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J)
1085                      else
1086                         if (SICE(I,J,KRT) .lt. SMCMAX1(I,J)) then
1087                            SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J) - &
1088                                                        (SICE(I,J,KRT) - WATHOLDCAP)
1089                         endif
1090                         if (SICE(I,J,KRT) .ge. SMCMAX1(I,J)) then
1091                            SMCWLTRT(IXXRT,JYYRT,KRT) = 0.
1092                         endif
1093                      endif
1094                   ELSE ! no ice
1095                      SMCMAXRT(IXXRT,JYYRT,KRT) = SMCMAX1(I,J)
1096                      SMCREFRT(IXXRT,JYYRT,KRT) = SMCREF1(I,J)
1097                      WATHOLDCAP = SMCMAX1(I,J) - SMCWLT1(I,J) !TODO: Not used again so can delete
1098                      SMCWLTRT(IXXRT,JYYRT,KRT) = SMCWLT1(I,J)
1099                   ENDIF   ! endif adjust for soil ice
1100
1101                   !Now Adjust soil moisture
1102                   !DJG Use SH2O instead of SMC for 'liquid' water...
1103                   IF (SMCMAXRT(IXXRT,JYYRT,KRT) .GT. 0) THEN !Check for smcmax data (=0 over water)
1104                      SMCRT(IXXRT,JYYRT,KRT) = SH2OX(I,J,KRT) * SH2OWGT(IXXRT,JYYRT,KRT)
1105                      !old SMCRT(IXXRT,JYYRT,KRT)=SMC(I,J,KRT)
1106                   ELSE
1107                      !TODO: Double-check the values for water cells
1108                      SMCRT(IXXRT,JYYRT,KRT) = 0.001  !will be skipped w/ landmask
1109                      SMCMAXRT(IXXRT,JYYRT,KRT) = 0.001
1110                   END IF
1111
1112                   !DJG Check/Adjust so that subgrid cells do not exceed saturation...
1113                   IF (SMCRT(IXXRT,JYYRT,KRT) .GT. SMCMAXRT(IXXRT,JYYRT,KRT)) THEN
1114                      SMCEXCS = (SMCRT(IXXRT,JYYRT,KRT) - SMCMAXRT(IXXRT,JYYRT,KRT)) &
1115                                * SLDPTH(KRT) * 1000.  !Excess soil water in units of (mm)
1116                      SMCRT(IXXRT,JYYRT,KRT) = SMCMAXRT(IXXRT,JYYRT,KRT)
1117                      DO KF = KRT-1,1,-1  !loop back upward to redistribute excess water from disagg.
1118                         SMCRT(IXXRT,JYYRT,KF) = SMCRT(IXXRT,JYYRT,KF) + SMCEXCS/(SLDPTH(KF)*1000.)
1119                         !Recheck new lyr sat.
1120                         IF (SMCRT(IXXRT,JYYRT,KF) .GT. SMCMAXRT(IXXRT,JYYRT,KF)) THEN
1121                            SMCEXCS = (SMCRT(IXXRT,JYYRT,KF) - SMCMAXRT(IXXRT,JYYRT,KF)) &
1122                                      * SLDPTH(KF) * 1000.  !Excess soil water in units of (mm)
1123                            SMCRT(IXXRT,JYYRT,KF) = SMCMAXRT(IXXRT,JYYRT,KF)
1124                         ELSE  ! Excess soil water expired
1125                            SMCEXCS = 0.
1126                            EXIT
1127                         END IF
1128                      END DO
1129                      IF (SMCEXCS .GT. 0) THEN  !If not expired by sfc then add to Infil. Excess
1130                         INFXSUBRT(IXXRT,JYYRT) = INFXSUBRT(IXXRT,JYYRT) + SMCEXCS
1131                         SMCEXCS = 0.
1132                      END IF
1133                   END IF  !End if for soil moisture saturation excess
1134
1135                end do !End do for soil profile loop
1136
1137                ! Debug loop
1138                do KRT=1,NSOIL
1139                   IF (SMCRT(IXXRT,JYYRT,KRT) .GT. SMCMAXRT(IXXRT,JYYRT,KRT)) THEN
1140                      print *, "FATAL ERROR: SMCMAX exceeded upon disaggregation3...", &
1141                               ixxrt, jyyrt, krt, &
1142                               SMCRT(IXXRT,JYYRT,KRT),SMCMAXRT(IXXRT,JYYRT,KRT)
1143                      call hydro_stop("In disaggregateDomain() - SMCMAX exceeded upon disaggregation3")
1144                   ELSE IF (SMCRT(IXXRT,JYYRT,KRT).LE.0.) THEN
1145                      print *, "SMCRT fully depleted upon disaggregation...", &
1146                          ixxrt,jyyrt,krt,&
1147                          "SMCRT=",SMCRT(IXXRT,JYYRT,KRT),"SH2OWGT=",SH2OWGT(IXXRT,JYYRT,KRT),&
1148                          "SH2O=",SH2OX(I,J,KRT)
1149                      print*, "SMC=", SMC(i,j,KRT), "SICE =", sice(i,j,KRT)
1150                      print *, "VEGTYP = ", VEGTYP(I,J)
1151                      print *, "i,j,krt, nsoil",i,j,krt,nsoil
1152                      ! ADCHANGE: If values are close but not exact, end up with a crash.
1153                      ! Force values to match.
1154                      !IF (SMC(i,j,KRT).EQ.sice(i,j,KRT)) THEN
1155                      IF (ABS(SMC(i,j,KRT) - sice(i,j,KRT)) .LE. 0.00001) THEN
1156                         print *, "SMC = SICE, soil layer totally frozen, proceeding..."
1157                         SMCRT(IXXRT,JYYRT,KRT) = 0.001
1158                         sice(i,j,KRT) = SMC(i,j,KRT)
1159                      ELSE
1160                         call hydro_stop("In disaggregateDomain() - SMCRT depleted")
1161                      END IF
1162                   END IF
1163                end do !debug loop
1164
1165                ! Now do simple grid remapping tasks
1166
1167                ! Lateral ksat
1168                !DJG 6.12.08 Map lateral hydraulic conductivity and apply distributed scaling
1169                ! ---        factor that will be read in from hires terrain file
1170                !              LKSATRT(IXXRT,JYYRT) = LKSAT(I,J) * LKSATFAC(IXXRT,JYYRT) * &  !Apply scaling factor...
1171                ! ...and scale from max to 0 when SMC decreases from SMCMAX to SMCREF...
1172                !!DJG error found from KIT,improper scaling       ((SMCMAXRT(IXXRT,JYYRT,NSOIL) - SMCRT(IXXRT,JYYRT,NSOIL)) / &
1173                !                                    (max(0.,(SMCMAXRT(IXXRT,JYYRT,NSOIL) - SMCRT(IXXRT,JYYRT,NSOIL))) / &
1174                !                                    (SMCMAXRT(IXXRT,JYYRT,NSOIL)-SMCREFRT(IXXRT,JYYRT,NSOIL)) )
1175                ! AD_CHANGE:
1176                ! Corrected to scale from 0 at SMCREF to full LKSAT*LKSATFAC at SMCMAX.
1177                ! This is a linear simplification of the true k to smc relationship to cover the
1178                ! range of conditions between smcmax and smcref. The DHSVM lateral routing
1179                ! scheme assumes saturated conditions (and therfore ksat), but we are extending
1180                ! that slightly to route water up until smcref.
1181                smScaleFact = (SMCRT(IXXRT,JYYRT,NSOIL) - SMCREFRT(IXXRT,JYYRT,NSOIL)) / &
1182                                 (SMCMAXRT(IXXRT,JYYRT,NSOIL) - SMCREFRT(IXXRT,JYYRT, NSOIL))
1183                smScaleFact = max(0., smScaleFact) !becomes 0 if less than SMCREF
1184                smScaleFact = min(1., smScaleFact) !make sure scale factor doesn't go over 1
1185                LKSATRT(IXXRT,JYYRT) = LKSAT(I,J) * LKSATFAC(IXXRT,JYYRT) * smScaleFact
1186
1187                ! n exponent for subsurface routing
1188                nexprt(ixxrt,jyyrt) = nexp(i,j)
1189
1190                ! Lake mask
1191                !DJG set up lake mask...
1192                !--- modify to make lake mask large here, but not one of the routed lakes!!!
1193                !--            IF (VEGTYP(I,J).eq.16) then
1194                IF (VEGTYP(I,J) .eq. iswater .and. CH_NETRT(IXXRT,JYYRT).le.0) then
1195                   !--LAKE_MSKRT(IXXRT,JYYRT) = 1
1196                   !yw  LAKE_MSKRT(IXXRT,JYYRT) = 9999
1197                   LAKE_MSKRT(IXXRT,JYYRT) = -9999
1198                end if
1199
1200                ! Soil type
1201                ! BF disaggregate soiltype information for gw-soil-coupling
1202                ! TODO: move this disaggregation code line to lsm_init section because soiltype is time-invariant
1203                soiltypRT(ixxrt,jyyrt) = soiltyp(i,j)
1204
1205             end do ! end disagg fine grid i loop
1206          end do ! end disagg fine grid j loop
1207
1208       end do ! end coarse grid i loop
1209    end do ! end coarse fine grid j loop
1210
1211    ! AD: Add new zeroing out of -9999 elevation cells which are ocean
1212    where (ELRT .lt. -9998)
1213       OCEAN_INFXSUBRT = INFXSUBRT
1214       INFXSUBRT = 0.0
1215    endwhere
1216
1217 #ifdef HYDRO_D
1218    ! ADCHANGE: START Final water balance variables
1219    ! ALL VARS in MM
1220    suminfxsrt2 = 0.
1221    smcrttot2 = 0.
1222    do ii=1,IXRT
1223       do jj=1,JXRT
1224          suminfxsrt2 = suminfxsrt2 + INFXSUBRT(ii,jj) / float(IXRT*JXRT)
1225          do kk=1,NSOIL
1226             smcrttot2 = smcrttot2 + SMCRT(ii,jj,kk)*SLDPTH(kk)*1000. / float(IXRT*JXRT)
1227          end do
1228       end do
1229    end do
1230 #ifdef MPP_LAND
1231    ! not tested
1232    CALL sum_real1(suminfxsrt2)
1233    CALL sum_real1(smcrttot2)
1234    suminfxsrt2 = suminfxsrt2/float(numprocs)
1235    smcrttot2 = smcrttot2/float(numprocs)
1236 #endif
1237 #ifdef MPP_LAND
1238    if (my_id .eq. IO_id) then
1239 #endif
1240       print *, "Disagg Mass Bal: "
1241       print *, "WB_DISAG!InfxsDiff", suminfxsrt2-suminfxs1
1242       print *, "WB_DISAG!Infxs1", suminfxs1
1243       print *, "WB_DISAG!Infxs2", suminfxsrt2
1244       print *, "WB_DISAG!SMCDIff", smcrttot2-(smctot1-sicetot1)
1245       print *, "WB_DISAG!SMC1", smctot1
1246       print *, "WB_DISAG!SICE1", sicetot1
1247       print *, "WB_DISAG!SMC2", smcrttot2
1248       print *, "WB_DISAG!Residual", (suminfxsrt2-suminfxs1) + (smcrttot2-(smctot1-sicetot1))
1249 #ifdef MPP_LAND
1250    endif
1251 #endif
1252    ! END Final water balance variables
1253 #endif
1254
1255 #ifdef HYDRO_D
1256    print *, "After Disaggregation..."
1257 #endif
1258 #ifdef MPP_LAND
1259    call MPP_LAND_COM_REAL(INFXSUBRT,IXRT,JXRT,99)
1260    call MPP_LAND_COM_REAL(LKSATRT,IXRT,JXRT,99)
1261    call MPP_LAND_COM_REAL(NEXPRT,IXRT,JXRT,99)
1262    call MPP_LAND_COM_INTEGER(LAKE_MSKRT,IXRT,JXRT,99)
1263    do i = 1, NSOIL
1264       call MPP_LAND_COM_REAL(SMCMAXRT(:,:,i),IXRT,JXRT,99)
1265       call MPP_LAND_COM_REAL(SMCRT(:,:,i),IXRT,JXRT,99)
1266       call MPP_LAND_COM_REAL(SMCWLTRT(:,:,i),IXRT,JXRT,99)
1267    end DO
1268 #endif
1269
1270 end subroutine disaggregateDomain
1271 !===================================================================================================
1272
1273
1274 subroutine SubsurfaceRouting_drv(did)
1275
1276     use module_RT_data, only: rt_domain
1277     use config_base, only: nlst
1278
1279     implicit none
1280     integer :: did
1281     IF (nlst(did)%SUBRTSWCRT.EQ.1) THEN
1282         call subsurfaceRouting ( rt_domain(did)%subsurface, &
1283             rt_domain(did)%subsurface_static, &
1284             rt_domain(did)%subsurface_input, &
1285             rt_domain(did)%subsurface_output, &
1286             rt_domain(did)%overland%control%infiltration_excess)
1287     endif
1288
1289 end subroutine SubsurfaceRouting_drv
1290
1291
1292
1293 subroutine OverlandRouting_drv(did)
1294     use module_RT_data, only: rt_domain
1295     use config_base, only: nlst
1296
1297     implicit none
1298     integer :: did
1299     if(nlst(did)%OVRTSWCRT .eq. 1) then
1300         !call OverlandRouting (nlst_rt(did)%DT, nlst_rt(did)%DTRT_TER, nlst_rt(did)%rt_option, &
1301             !         rt_domain(did)%ixrt, rt_domain(did)%jxrt,rt_domain(did)%overland%streams_and_lakes%lake_mask, &
1302             !         rt_domain(did)%overland%control%infiltration_excess, rt_domain(did)%overland%properties%retention_depth,rt_domain(did)%overland%properties%roughness, &
1303             !         rt_domain(did)%overland%properties%surface_slope_x, rt_domain(did)%overland%properties%surface_slope_y, rt_domain(did)%overland%control%surface_water_head_routing,  &
1304             !         rt_domain(did)%overland%control%dhrt, rt_domain(did)%overland%streams_and_lakes%ch_netrt, rt_domain(did)%overland%streams_and_lakes%surface_water_to_channel,&
1305             !         rt_domain(did)%overland%streams_and_lakes%surface_water_to_lake,rt_domain(did)%overland%control%boundary_flux, &
1306             !         rt_domain(did)%overland%streams_and_lakes%accumulated_surface_water_to_channel,rt_domain(did)%overland%control%boundary_flux_total, rt_domain(did)%overland%streams_and_lakes%accumulated_surface_water_to_lake,&
1307             !         rt_domain(did)%q_sfcflx_x,rt_domain(did)%q_sfcflx_y, &
1308             !         rt_domain(did)%overland%properties%distance_to_neighbor, rt_domain(did)%overland%properties%surface_slope, rt_domain(did)%overland%properties%max_surface_slope_index , &
1309             !         rt_domain(did)%overland%mass_balance%post_soil_moisture_content,rt_domain(did)%overland%mass_balance%pre_infiltration_excess,rt_domain(did)%overland%mass_balance%post_infiltration_excess, &
1310             !         rt_domain(did)%overland%mass_balance%pre_soil_moisture_content,rt_domain(did)%overland%mass_balance%accumulated_change_in_soil_moisture )
1311         call OverlandRouting( &
1312             rt_domain(did)%overland, &
1313             nlst(did)%DT, &
1314             nlst(did)%DTRT_TER, &
1315             nlst(did)%rt_option, &
1316             rt_domain(did)%ixrt, &
1317             rt_domain(did)%jxrt, &
1318             rt_domain(did)%q_sfcflx_x, &
1319             rt_domain(did)%q_sfcflx_y &
1320             )
1321         ! ADCHANGE: If overland routing is called, INFXSUBRT is moved to SFCHEADSUBRT, so
1322         !           zeroing out just in case
1323         rt_domain(did)%overland%control%infiltration_excess = 0.0
1324     endif
1325 end subroutine OverlandRouting_drv
1326
1327
1328
1329
1330
1331
1332 subroutine time_seconds(i3)
1333     integer time_array(8)
1334     real*8 i3
1335     call date_and_time(values=time_array)
1336     i3 = time_array(4)*24*3600+time_array(5) * 3600 + time_array(6) * 60 + &
1337         time_array(7) + 0.001 * time_array(8)
1338     return
1339 end subroutine time_seconds