hydro/Routing/module_lsm_forcing.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 module module_lsm_forcing
  22
  23 #ifdef MPP_LAND
  24     use module_mpp_land
  25 #endif
  26     use module_HYDRO_io, only: get_2d_netcdf, get_soilcat_netcdf, get2d_int
  27     use module_hydro_stop, only:HYDRO_stop
  28
  29 implicit none
  30 #include <netcdf.inc>
  31     integer :: i_forcing
  32 character(len=19) out_date
  33
  34 interface read_hydro_forcing
  35 #ifdef MPP_LAND
  36    !yw module procedure read_hydro_forcing_mpp
  37    module procedure read_hydro_forcing_mpp1
  38 #else
  39    module procedure read_hydro_forcing_seq
  40 #endif
  41 end interface
  42
  43 Contains
  44
  45   subroutine READFORC_WRF(flnm,ix,jx,target_date,t,q,u,v,p,lw,sw,pcp,lai,fpar)
  46
  47     implicit none
  48     character(len=*),                   intent(in)  :: flnm
  49     integer,                            intent(in)  :: ix
  50     integer,                            intent(in)  :: jx
  51     character(len=*),                   intent(in)  :: target_date
  52     real,             dimension(ix,jx) :: t,q,u,v,p,lw,sw,pcp,pcpc, lai,fpar
  53     integer   tlevel
  54
  55     character(len=256) :: units
  56     integer :: ierr
  57     integer :: ncid
  58
  59     tlevel = 1
  60
  61     pcp = 0
  62     pcpc = 0
  63
  64     ! Open the NetCDF file.
  65     ierr = nf_open(flnm, NF_NOWRITE, ncid)
  66     if (ierr /= 0) then
  67        write(*,'("READFORC_WRF Problem opening netcdf file: ''", A, "''")') trim(flnm)
  68        call hydro_stop("In READFORC_WRF() - Problem opening netcdf file")
  69     endif
  70
  71     call get_2d_netcdf_ruc("T2",     ncid, t,     ix, jx,tlevel, .true., ierr)
  72     call get_2d_netcdf_ruc("Q2",     ncid, q,     ix, jx,tlevel, .true., ierr)
  73     call get_2d_netcdf_ruc("U10",    ncid, u,     ix, jx,tlevel, .true., ierr)
  74     call get_2d_netcdf_ruc("V10",    ncid, v,     ix, jx,tlevel, .true., ierr)
  75     call get_2d_netcdf_ruc("PSFC",   ncid, p,     ix, jx,tlevel, .true., ierr)
  76     call get_2d_netcdf_ruc("GLW",    ncid, lw,    ix, jx,tlevel, .true., ierr)
  77     call get_2d_netcdf_ruc("SWDOWN", ncid, sw,    ix, jx,tlevel, .true., ierr)
  78     call get_2d_netcdf_ruc("RAINC",  ncid, pcpc,  ix, jx,tlevel, .true., ierr)
  79     call get_2d_netcdf_ruc("RAINNC", ncid, pcp,   ix, jx,tlevel, .true., ierr)
  80     call get_2d_netcdf_ruc("VEGFRA", ncid, fpar,  ix, jx,tlevel, .true., ierr)
  81     if(ierr == 0) then
  82         if(maxval(fpar) .gt. 10 .and. (maxval(fpar) .lt. 10000) ) fpar = fpar*0.01
  83     endif
  84     call get_2d_netcdf_ruc("LAI", ncid, lai,  ix, jx,tlevel, .true., ierr)
  85
  86     ierr = nf_close(ncid)
  87
  88 !DJG  Add the convective and non-convective rain components (note: conv. comp=0
  89 !for cloud resolving runs...)
  90 !DJG  Note that for WRF these are accumulated values to be adjusted to rates in
  91 !driver...
  92
  93     pcp=pcp+pcpc   ! assumes pcpc=0 for resolved convection...
  94
  95   end subroutine READFORC_WRF
  96
  97   subroutine read_hrldas_hdrinfo(geo_static_flnm, ix, jx, land_cat, soil_cat)
  98     ! Simply return the dimensions of the grid.
  99     implicit none
 100     character(len=*),          intent(in)  :: geo_static_flnm
 101     integer, intent(out) :: ix, jx, land_cat, soil_cat ! dimensions
 102
 103     integer :: iret, ncid, dimid
 104
 105     ! Open the NetCDF file.
 106     iret = nf_open(geo_static_flnm, NF_NOWRITE, ncid)
 107     if (iret /= 0) then
 108        write(*,'("Problem opening geo_static file: ''", A, "''")') &
 109             trim(geo_static_flnm)
 110        call hydro_stop("In read_hrldas_hdrinfo() - Problem opening geo_static file")
 111     endif
 112
 113     iret = nf_inq_dimid(ncid, "west_east", dimid)
 114
 115     if (iret /= 0) then
 116 !       print*, "nf_inq_dimid:  west_east"
 117        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimid:  west_east problem")
 118     endif
 119
 120     iret = nf_inq_dimlen(ncid, dimid, ix)
 121     if (iret /= 0) then
 122 !       print*, "nf_inq_dimlen:  west_east"
 123        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimlen:  west_east problem")
 124     endif
 125
 126     iret = nf_inq_dimid(ncid, "south_north", dimid)
 127     if (iret /= 0) then
 128 !       print*, "nf_inq_dimid:  south_north"
 129        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimid:  south_north problem")
 130     endif
 131
 132     iret = nf_inq_dimlen(ncid, dimid, jx)
 133     if (iret /= 0) then
 134  !      print*, "nf_inq_dimlen:  south_north"
 135        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimlen:  south_north problem")
 136     endif
 137
 138     iret = nf_inq_dimid(ncid, "land_cat", dimid)
 139     if (iret /= 0) then
 140  !      print*, "nf_inq_dimid:  land_cat"
 141        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimid:  land_cat problem")
 142     endif
 143
 144     iret = nf_inq_dimlen(ncid, dimid, land_cat)
 145     if (iret /= 0) then
 146        print*, "nf_inq_dimlen:  land_cat"
 147        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimlen:  land_cat problem")
 148     endif
 149
 150     iret = nf_inq_dimid(ncid, "soil_cat", dimid)
 151     if (iret /= 0) then
 152  !      print*, "nf_inq_dimid:  soil_cat"
 153        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimid:  soil_cat problem")
 154     endif
 155
 156     iret = nf_inq_dimlen(ncid, dimid, soil_cat)
 157     if (iret /= 0) then
 158  !      print*, "nf_inq_dimlen:  soil_cat"
 159        call hydro_stop("In read_hrldas_hdrinfo() - nf_inq_dimlen:  soil_cat problem")
 160     endif
 161
 162     iret = nf_close(ncid)
 163
 164   end subroutine read_hrldas_hdrinfo
 165
 166
 167
 168   subroutine readland_hrldas(geo_static_flnm,ix,jx,land_cat,soil_cat,vegtyp,soltyp, &
 169                   terrain,latitude,longitude,SOLVEG_INITSWC)
 170     implicit none
 171     character(len=*),          intent(in)  :: geo_static_flnm
 172     integer,                   intent(in)  :: ix, jx, land_cat, soil_cat,SOLVEG_INITSWC
 173     integer, dimension(ix,jx), intent(out) :: vegtyp, soltyp
 174     real,    dimension(ix,jx), intent(out) :: terrain, latitude, longitude
 175
 176     character(len=256) :: units
 177     integer :: ierr,i,j,jj
 178     integer :: ncid,varid
 179     real, dimension(ix,jx) :: xdum
 180     integer, dimension(ix,jx) :: vegtyp_inv, soiltyp_inv,xdum_int
 181     integer flag ! flag = 1 from wrfsi, flag =2 from WPS.
 182     CHARACTER(len=256)       :: var_name
 183     integer :: islake, iswater, isoilwater
 184
 185     ! Open the NetCDF file.
 186     ierr = nf_open(geo_static_flnm, NF_NOWRITE, ncid)
 187
 188     if (ierr /= 0) then
 189        write(*,'("Problem opening geo_static file: ''", A, "''")') trim(geo_static_flnm)
 190        call hydro_stop("In readland_hrldas() - Problem opening geo_static file")
 191     endif
 192
 193     flag = -99
 194     ierr = nf_inq_varid(ncid,"XLAT", varid)
 195     flag = 1
 196     if(ierr .ne. 0) then
 197         ierr = nf_inq_varid(ncid,"XLAT_M", varid)
 198         if(ierr .ne. 0) then
 199 !            write(6,*) "XLAT not found from wrfstatic file. "
 200             call hydro_stop("In readland_hrldas() - XLAT not found from wrfstatic file")
 201         endif
 202         flag = 2
 203     endif
 204
 205     ! Get Latitude (lat)
 206     if(flag .eq. 1) then
 207        call get_2d_netcdf("XLAT", ncid, latitude,  units, ix, jx, .TRUE., ierr)
 208     else
 209       call get_2d_netcdf("XLAT_M", ncid, latitude,  units, ix, jx, .TRUE., ierr)
 210     endif
 211
 212     ! Get Longitude (lon)
 213     if(flag .eq. 1) then
 214         call get_2d_netcdf("XLONG", ncid, longitude, units, ix, jx, .TRUE., ierr)
 215     else
 216         call get_2d_netcdf("XLONG_M", ncid, longitude, units, ix, jx, .TRUE., ierr)
 217     endif
 218
 219     ! Get Terrain (avg)
 220     if(flag .eq. 1) then
 221        call get_2d_netcdf("HGT", ncid, terrain,   units, ix, jx, .TRUE., ierr)
 222     else
 223         call get_2d_netcdf("HGT_M", ncid, terrain,   units, ix, jx, .TRUE., ierr)
 224     endif
 225
 226
 227     if (SOLVEG_INITSWC.eq.0) then
 228 !      ! Get Dominant Land Use categories (use)
 229 !      call get_landuse_netcdf(ncid, xdum ,   units, ix, jx, land_cat)
 230 !      vegtyp = nint(xdum)
 231
 232      var_name = "LU_INDEX"
 233          call get2d_int(var_name,xdum_int,ix,jx,&
 234                trim(geo_static_flnm))
 235          vegtyp = xdum_int
 236
 237       ! Get Dominant Soil Type categories in the top layer (stl)
 238       call get_soilcat_netcdf(ncid, xdum ,   units, ix, jx, soil_cat)
 239       soltyp = nint(xdum)
 240
 241     else if (SOLVEG_INITSWC.eq.1) then
 242        var_name = "VEGTYP"
 243        call get2d_int(var_name,VEGTYP_inv,ix,jx,&
 244               trim(geo_static_flnm))
 245
 246        var_name = "SOILTYP"
 247        call get2d_int(var_name,SOILTYP_inv,ix,jx,&
 248               trim(geo_static_flnm))
 249        do i=1,ix
 250          jj=jx
 251          do j=1,jx
 252            VEGTYP(i,j)=VEGTYP_inv(i,jj)
 253            SOLTYP(i,j)=SOILTYP_inv(i,jj)
 254            jj=jx-j
 255          end do
 256        end do
 257
 258     endif
 259
 260     ierr = NF_GET_ATT_INT(ncid, NF_GLOBAL, 'ISWATER', iswater)
 261     if (ierr /= 0) then
 262        call hydro_stop("In readland_hrldas() - Attribute ISWATER unable to be read from geo_static_flnm")
 263     endif
 264
 265     ierr = NF_GET_ATT_INT(ncid, NF_GLOBAL, 'ISOILWATER', isoilwater)
 266     if (ierr /= 0) then
 267        call hydro_stop("In readland_hrldas() - Attribute ISOILWATER unable to be read from geo_static_flnm")
 268     endif
 269
 270     ierr = NF_GET_ATT_INT(ncid, NF_GLOBAL, 'ISLAKE', islake)
 271     if (ierr /= 0) then
 272        call hydro_stop("In readland_hrldas() - Attribute ISLAKE unable to be read from geo_static_flnm")
 273     endif
 274
 275     ! Close the NetCDF file
 276     ierr = nf_close(ncid)
 277     if (ierr /= 0) then
 278        print*, "MODULE_NOAHLSM_HRLDAS_INPUT:  READLAND_HRLDAS:  NF_CLOSE"
 279        call hydro_stop("In readland_hrldas() - NF_CLOSE problem")
 280     endif
 281
 282  write(6, *) "readland_hrldas: ISLAKE ISWATER ISOILWATER", islake, iswater, isoilwater
 283
 284     ! Make sure vegtyp and soltyp are consistent when it comes to water points,
 285     ! by setting soil category to water when vegetation category is water, and
 286     ! vice-versa.
 287     where (vegtyp == islake) vegtyp = iswater
 288     where (vegtyp == iswater) soltyp = isoilwater
 289     where (soltyp == isoilwater) vegtyp = iswater
 290
 291 !DJG test for deep gw function...
 292 !    where (soltyp <> 14) soltyp = 1
 293
 294   end subroutine readland_hrldas
 295
 296
 297       subroutine get_2d_netcdf_ruc(var_name,ncid,var, &
 298             ix,jx,tlevel,fatal_if_error,ierr)
 299           character(len=*), intent(in) :: var_name
 300           integer,intent(in) ::  ncid,ix,jx,tlevel
 301           real, intent(out):: var(ix,jx)
 302           logical, intent(in) :: fatal_if_error
 303           integer dims(4), dim_len(4)
 304           integer ierr,iret
 305           integer varid
 306            integer start(4),count(4)
 307            data count /1,1,1,1/
 308            data start /1,1,1,1/
 309           count(1) = ix
 310           count(2) = jx
 311           start(4) = tlevel
 312       ierr = nf_inq_varid(ncid,  var_name,  varid)
 313
 314       if (ierr /= 0) then
 315         if (fatal_IF_ERROR) then
 316            print*, "MODULE_NOAHLSM_HRLDAS_INPUT: get_2d_netcdf_ruc:nf_inq_varid ", trim(var_name)
 317            call hydro_stop("In get_2d_netcdf_ruc() - nf_inq_varid problem")
 318         else
 319           return
 320         endif
 321       endif
 322
 323       ierr = nf_get_vara_real(ncid, varid, start,count,var)
 324
 325
 326       return
 327       end subroutine get_2d_netcdf_ruc
 328
 329
 330       subroutine get_2d_netcdf_cows(var_name,ncid,var, &
 331             ix,jx,tlevel,fatal_if_error,ierr)
 332           character(len=*), intent(in) :: var_name
 333           integer,intent(in) ::  ncid,ix,jx,tlevel
 334           real, intent(out):: var(ix,jx)
 335           logical, intent(in) :: fatal_if_error
 336           integer ierr, iret
 337           integer varid
 338           integer start(4),count(4)
 339           data count /1,1,1,1/
 340           data start /1,1,1,1/
 341           count(1) = ix
 342           count(2) = jx
 343           start(4) = tlevel
 344       iret = nf_inq_varid(ncid,  var_name,  varid)
 345
 346       if (iret /= 0) then
 347         if (fatal_IF_ERROR) then
 348            print*, "MODULE_NOAHLSM_HRLDAS_INPUT: get_2d_netcdf_cows:nf_inq_varid"
 349            call hydro_stop("In get_2d_netcdf_cows() - nf_inq_varid problem")
 350         else
 351           ierr = iret
 352           return
 353         endif
 354       endif
 355       iret = nf_get_vara_real(ncid, varid, start,count,var)
 356
 357       return
 358       end subroutine get_2d_netcdf_cows
 359
 360
 361
 362
 363
 364   subroutine readinit_hrldas(netcdf_flnm, ix, jx, nsoil, target_date, &
 365        smc, stc, sh2o, cmc, t1, weasd, snodep)
 366     implicit none
 367     character(len=*),                intent(in)  :: netcdf_flnm
 368     integer,                         intent(in)  :: ix
 369     integer,                         intent(in)  :: jx
 370     integer,                         intent(in)  :: nsoil
 371     character(len=*),                intent(in)  :: target_date
 372     real,    dimension(ix,jx,nsoil), intent(out) :: smc
 373     real,    dimension(ix,jx,nsoil), intent(out) :: stc
 374     real,    dimension(ix,jx,nsoil), intent(out) :: sh2o
 375     real,    dimension(ix,jx),       intent(out) :: cmc
 376     real,    dimension(ix,jx),       intent(out) :: t1
 377     real,    dimension(ix,jx),       intent(out) :: weasd
 378     real,    dimension(ix,jx),       intent(out) :: snodep
 379
 380     character(len=256) :: units
 381     character(len=8) :: name
 382     integer :: ix_read, jx_read,i,j
 383
 384     integer :: ierr, ncid, ierr_snodep
 385     integer :: idx
 386
 387     logical :: found_canwat, found_skintemp, found_weasd, found_stemp, found_smois
 388
 389     ! Open the NetCDF file.
 390     ierr = nf_open(netcdf_flnm, NF_NOWRITE, ncid)
 391     if (ierr /= 0) then
 392        write(*,'("READINIT Problem opening netcdf file: ''", A, "''")') &
 393             trim(netcdf_flnm)
 394        call hydro_stop("In readinit_hrldas()- Problem opening netcdf file")
 395     endif
 396
 397     call get_2d_netcdf("CANWAT",     ncid, cmc,     units, ix, jx, .TRUE., ierr)
 398     call get_2d_netcdf("SKINTEMP",   ncid, t1,      units, ix, jx, .TRUE., ierr)
 399     call get_2d_netcdf("WEASD",      ncid, weasd,   units, ix, jx, .TRUE., ierr)
 400
 401     if (trim(units) == "m") then
 402        ! No conversion necessary
 403     else if (trim(units) == "mm") then
 404        ! convert WEASD from mm to m
 405        weasd = weasd * 1.E-3
 406     else
 407        print*, 'units = "'//trim(units)//'"'
 408 !       print*, "Unrecognized units on WEASD"
 409        call hydro_stop("In readinit_hrldas() - Unrecognized units on WEASD")
 410     endif
 411
 412     call get_2d_netcdf("SNODEP",     ncid, snodep,   units, ix, jx, .FALSE., ierr_snodep)
 413     call get_2d_netcdf("STEMP_1",    ncid, stc(:,:,1), units,  ix, jx, .TRUE., ierr)
 414     call get_2d_netcdf("STEMP_2",    ncid, stc(:,:,2), units,  ix, jx, .TRUE., ierr)
 415     call get_2d_netcdf("STEMP_3",    ncid, stc(:,:,3), units,  ix, jx, .TRUE., ierr)
 416     call get_2d_netcdf("STEMP_4",    ncid, stc(:,:,4), units,  ix, jx, .TRUE., ierr)
 417     call get_2d_netcdf("SMOIS_1",    ncid, smc(:,:,1), units,  ix, jx, .TRUE., ierr)
 418     call get_2d_netcdf("SMOIS_2",    ncid, smc(:,:,2), units,  ix, jx, .TRUE., ierr)
 419     call get_2d_netcdf("SMOIS_3",    ncid, smc(:,:,3), units,  ix, jx, .TRUE., ierr)
 420     call get_2d_netcdf("SMOIS_4",    ncid, smc(:,:,4), units,  ix, jx, .TRUE., ierr)
 421
 422
 423     if (ierr_snodep /= 0) then
 424        ! Quick assumption regarding snow depth.
 425        snodep = weasd * 10.
 426     endif
 427
 428
 429 !DJG check for erroneous neg WEASD or SNOWD due to offline interpolation...
 430        do i=1,ix
 431          do j=1,jx
 432            if (WEASD(i,j).lt.0.) WEASD(i,j)=0.0  !set lower bound to correct bi-lin interp err...
 433            if (snodep(i,j).lt.0.) snodep(i,j)=0.0  !set lower bound to correct bi-lin interp err...
 434          end do
 435        end do
 436
 437
 438     sh2o = smc
 439
 440     ierr = nf_close(ncid)
 441   end subroutine readinit_hrldas
 442
 443
 444
 445
 446   subroutine READFORC_HRLDAS(flnm,ix,jx,target_date, t,q,u,v,p,lw,sw,pcp,lai,fpar)
 447     implicit none
 448
 449     character(len=*),                   intent(in)  :: flnm
 450     integer,                            intent(in)  :: ix
 451     integer,                            intent(in)  :: jx
 452     character(len=*),                   intent(in)  :: target_date
 453     real,             dimension(ix,jx), intent(out) :: t
 454     real,             dimension(ix,jx), intent(out) :: q
 455     real,             dimension(ix,jx), intent(out) :: u
 456     real,             dimension(ix,jx), intent(out) :: v
 457     real,             dimension(ix,jx), intent(out) :: p
 458     real,             dimension(ix,jx), intent(out) :: lw
 459     real,             dimension(ix,jx), intent(out) :: sw
 460     real,             dimension(ix,jx), intent(out) :: pcp
 461     real,             dimension(ix,jx), intent(inout) :: lai
 462     real,             dimension(ix,jx), intent(inout) :: fpar
 463
 464     character(len=256) :: units
 465     integer :: ierr
 466     integer :: ncid
 467
 468     ! Open the NetCDF file.
 469     ierr = nf_open(trim(flnm), NF_NOWRITE, ncid)
 470     if (ierr /= 0) then
 471        write(*,'("READFORC Problem opening netcdf file: ''", A, "''")') trim(flnm)
 472        call hydro_stop("In READFORC_HRLDAS() - Problem opening netcdf file")
 473     endif
 474
 475     call get_2d_netcdf("T2D",     ncid, t,     units, ix, jx, .TRUE., ierr)
 476     call get_2d_netcdf("Q2D",     ncid, q,     units, ix, jx, .TRUE., ierr)
 477     call get_2d_netcdf("U2D",     ncid, u,     units, ix, jx, .TRUE., ierr)
 478     call get_2d_netcdf("V2D",     ncid, v,     units, ix, jx, .TRUE., ierr)
 479     call get_2d_netcdf("PSFC",    ncid, p,     units, ix, jx, .TRUE., ierr)
 480     call get_2d_netcdf("LWDOWN",  ncid, lw,    units, ix, jx, .TRUE., ierr)
 481     call get_2d_netcdf("SWDOWN",  ncid, sw,    units, ix, jx, .TRUE., ierr)
 482     call get_2d_netcdf("RAINRATE",ncid, pcp,   units, ix, jx, .TRUE., ierr)
 483     call get_2d_netcdf("VEGFRA",  ncid, fpar,  units, ix, jx, .FALSE., ierr)
 484     if (ierr == 0) then
 485       if(maxval(fpar) .gt. 10 .and. maxval(fpar) .lt. 10000)  fpar = fpar * 1.E-2
 486     endif
 487     call get_2d_netcdf("LAI",     ncid, lai,   units, ix, jx, .FALSE., ierr)
 488
 489     ierr = nf_close(ncid)
 490
 491   end subroutine READFORC_HRLDAS
 492
 493
 494
 495   subroutine READFORC_DMIP(flnm,ix,jx,var)
 496     implicit none
 497
 498     character(len=*),                   intent(in)  :: flnm
 499     integer,                            intent(in)  :: ix
 500     integer,                            intent(in)  :: jx
 501     real,       dimension(ix,jx), intent(out)       :: var
 502     character(len=13)                               :: head
 503     integer                          :: ncols, nrows, cellsize
 504     real                             :: xllc, yllc, no_data
 505     integer                          :: i,j
 506     character(len=256)                              ::junk
 507
 508     open (77,file=trim(flnm),form="formatted",status="old")
 509
 510 !    read(77,732) head,ncols
 511 !    read(77,732) head,nrows
 512 !732        FORMAT(A13,I4)
 513 !    read(77,733) head,xllc
 514 !    read(77,733) head,yllc
 515 !733        FORMAT(A13,F16.9)
 516 !    read(77,732) head,cellsize
 517 !    read(77,732) head,no_data
 518
 519     read(77,*) junk
 520     read(77,*) junk
 521     read(77,*) junk
 522     read(77,*) junk
 523     read(77,*) junk
 524     read(77,*) junk
 525
 526     do j=jx,1,-1
 527       read(77,*) (var(I,J),I=1,ix)
 528     end do
 529     close(77)
 530
 531   end subroutine READFORC_DMIP
 532
 533
 534
 535   subroutine READFORC_MDV(flnm,ix,jx,pcp,mmflag,ierr_flg)
 536     implicit none
 537
 538     character(len=*),                   intent(in)  :: flnm
 539     integer,                            intent(in)  :: ix
 540     integer,                            intent(in)  :: jx
 541     integer,                            intent(out)  :: ierr_flg
 542     integer :: it,jew,zsn
 543     real,             dimension(ix,jx), intent(out) :: pcp
 544
 545     character(len=256) :: units
 546     integer :: ierr,i,j,i2,j2,varid
 547     integer :: ncid,mmflag
 548     real, dimension(ix,jx) :: temp
 549
 550     mmflag = 0   ! flag for units spec. (0=mm, 1=mm/s)
 551
 552
 553 !open NetCDF file...
 554         ierr_flg = nf_open(flnm, NF_NOWRITE, ncid)
 555         if (ierr_flg /= 0) then
 556 #ifdef HYDRO_D
 557           write(*,'("READFORC_MDV Problem opening netcdf file: ''",A,"''")') &
 558                 trim(flnm)
 559 #endif
 560            return
 561         end if
 562
 563         ierr = nf_inq_varid(ncid,  "precip",  varid)
 564         if(ierr /= 0) ierr_flg = ierr
 565         if (ierr /= 0) then
 566           ierr = nf_inq_varid(ncid,  "precip_rate",  varid)   !recheck variable name...
 567           if (ierr /= 0) then
 568             ierr = nf_inq_varid(ncid,  "RAINRATE",  varid)   !recheck variable name...
 569             if (ierr /= 0) then
 570 #ifdef HYDRO_D
 571               write(*,'("READFORC_MDV Problem reading precip netcdf file: ''", A,"''")') &
 572                  trim(flnm)
 573 #endif
 574             end if
 575           end if
 576           ierr_flg = ierr
 577           mmflag = 1
 578         end if
 579         ierr = nf_get_var_real(ncid, varid, pcp)
 580         ierr = nf_close(ncid)
 581
 582         if (ierr /= 0) then
 583 #ifdef HYDRO_D
 584           write(*,'("READFORC_MDV Problem reading netcdf file: ''", A,"''")') trim(flnm)
 585 #endif
 586         end if
 587
 588   end subroutine READFORC_MDV
 589
 590
 591
 592   subroutine READFORC_NAMPCP(flnm,ix,jx,pcp,k,product)
 593     implicit none
 594
 595     character(len=*),                   intent(in)  :: flnm
 596     integer,                            intent(in)  :: ix
 597     integer,                            intent(in)  :: jx
 598     integer,                            intent(in)  :: k
 599     character(len=*),                   intent(in)  :: product
 600     integer :: it,jew,zsn
 601     parameter(it =  496,jew = 449, zsn = 499)   ! domain 1
 602 !    parameter(it =  496,jew = 74, zsn = 109)   ! domain 2
 603     real,             dimension(it,jew,zsn) :: buf
 604     real,             dimension(ix,jx), intent(out) :: pcp
 605
 606     character(len=256) :: units
 607     integer :: ierr,i,j,i2,j2,varid
 608     integer :: ncid
 609     real, dimension(ix,jx) :: temp
 610
 611 !      varname = trim(product)
 612
 613 !open NetCDF file...
 614       if (k.eq.1.) then
 615         ierr = nf_open(flnm, NF_NOWRITE, ncid)
 616         if (ierr /= 0) then
 617           write(*,'("READFORC_NAMPCP1 Problem opening netcdf file: ''",A, "''")') &
 618               trim(flnm)
 619           call hydro_stop("In READFORC_NAMPCP() - Problem opening netcdf file")
 620         end if
 621
 622         ierr = nf_inq_varid(ncid,  trim(product),  varid)
 623         ierr = nf_get_var_real(ncid, varid, buf)
 624         ierr = nf_close(ncid)
 625
 626         if (ierr /= 0) then
 627           write(*,'("READFORC_NAMPCP2 Problem reading netcdf file: ''", A,"''")') &
 628              trim(flnm)
 629           call hydro_stop("In READFORC_NAMPCP() - Problem reading netcdf file")
 630         end if
 631       endif
 632 #ifdef HYDRO_D
 633       print *, "Data read in...",it,ix,jx,k
 634 #endif
 635
 636 ! Extract single time slice from dataset...
 637
 638       do i=1,ix
 639         do j=1,jx
 640           pcp(i,j) = buf(k,i,j)
 641         end do
 642       end do
 643
 644 !      call get_2d_netcdf_ruc("trmm",ncid, pcp, jx, ix,k, .true., ierr)
 645
 646   end subroutine READFORC_NAMPCP
 647
 648
 649
 650
 651   subroutine READFORC_COWS(flnm,ix,jx,target_date, t,q,u,p,lw,sw,pcp,tlevel)
 652     implicit none
 653
 654     character(len=*),                   intent(in)  :: flnm
 655     integer,                            intent(in)  :: ix
 656     integer,                            intent(in)  :: jx
 657     character(len=*),                   intent(in)  :: target_date
 658     real,             dimension(ix,jx), intent(out) :: t
 659     real,             dimension(ix,jx), intent(out) :: q
 660     real,             dimension(ix,jx), intent(out) :: u
 661     real,             dimension(ix,jx) :: v
 662     real,             dimension(ix,jx), intent(out) :: p
 663     real,             dimension(ix,jx), intent(out) :: lw
 664     real,             dimension(ix,jx), intent(out) :: sw
 665     real,             dimension(ix,jx), intent(out) :: pcp
 666     integer   tlevel
 667
 668     character(len=256) :: units
 669     integer :: ierr
 670     integer :: ncid
 671
 672     ! Open the NetCDF file.
 673     ierr = nf_open(flnm, NF_NOWRITE, ncid)
 674     if (ierr /= 0) then
 675        write(*,'("READFORC_COWS Problem opening netcdf file: ''", A, "''")') trim(flnm)
 676        call hydro_stop("In READFORC_COWS() - Problem opening netcdf file")
 677     endif
 678
 679     call get_2d_netcdf_cows("TA2",     ncid, t,     ix, jx,tlevel, .TRUE., ierr)
 680     call get_2d_netcdf_cows("QV2",     ncid, q,     ix, jx,tlevel, .TRUE., ierr)
 681     call get_2d_netcdf_cows("WSPD10",  ncid, u,     ix, jx,tlevel, .TRUE., ierr)
 682     call get_2d_netcdf_cows("PRES",    ncid, p,     ix, jx,tlevel, .TRUE., ierr)
 683     call get_2d_netcdf_cows("GLW",     ncid, lw,    ix, jx,tlevel, .TRUE., ierr)
 684     call get_2d_netcdf_cows("RSD",     ncid, sw,    ix, jx,tlevel, .TRUE., ierr)
 685     call get_2d_netcdf_cows("RAIN",    ncid, pcp,   ix, jx,tlevel, .TRUE., ierr)
 686 !yw   call get_2d_netcdf_cows("V2D",     ncid, v,     ix, jx,tlevel, .TRUE., ierr)
 687
 688     ierr = nf_close(ncid)
 689
 690   end subroutine READFORC_COWS
 691
 692
 693
 694
 695   subroutine READFORC_RUC(flnm,ix,jx,target_date,t,q,u,v,p,lw,sw,pcp)
 696
 697     implicit none
 698
 699     character(len=*),                   intent(in)  :: flnm
 700     integer,                            intent(in)  :: ix
 701     integer,                            intent(in)  :: jx
 702     character(len=*),                   intent(in)  :: target_date
 703     real,             dimension(ix,jx) :: t,q,u,v,p,lw,sw,pcp,pcpc
 704     integer   tlevel
 705
 706     character(len=256) :: units
 707     integer :: ierr
 708     integer :: ncid
 709
 710     tlevel = 1
 711
 712     ! Open the NetCDF file.
 713     ierr = nf_open(flnm, NF_NOWRITE, ncid)
 714     if (ierr /= 0) then
 715        write(*,'("READFORC_RUC Problem opening netcdf file: ''", A, "''")') trim(flnm)
 716        call hydro_stop("In READFORC_RUC() - Problem opening netcdf file")
 717     endif
 718
 719     call get_2d_netcdf_ruc("T2",     ncid, t,     ix, jx,tlevel, .true., ierr)
 720     call get_2d_netcdf_ruc("Q2",     ncid, q,     ix, jx,tlevel, .true., ierr)
 721     call get_2d_netcdf_ruc("U10",    ncid, u,     ix, jx,tlevel, .true., ierr)
 722     call get_2d_netcdf_ruc("V10",    ncid, v,     ix, jx,tlevel, .true., ierr)
 723     call get_2d_netcdf_ruc("PSFC",   ncid, p,     ix, jx,tlevel, .true., ierr)
 724     call get_2d_netcdf_ruc("GLW",    ncid, lw,    ix, jx,tlevel, .true., ierr)
 725     call get_2d_netcdf_ruc("SWDOWN", ncid, sw,    ix, jx,tlevel, .true., ierr)
 726     call get_2d_netcdf_ruc("RAINC",  ncid, pcpc,  ix, jx,tlevel, .true., ierr)
 727     call get_2d_netcdf_ruc("RAINNC", ncid, pcp,   ix, jx,tlevel, .true., ierr)
 728
 729     ierr = nf_close(ncid)
 730
 731
 732 !DJG  Add the convective and non-convective rain components (note: conv. comp=0
 733 !for cloud resolving runs...)
 734 !DJG  Note that for RUC these are accumulated values to be adjusted to rates in
 735 !driver...
 736
 737     pcp=pcp+pcpc   ! assumes pcpc=0 for resolved convection...
 738
 739   end subroutine READFORC_RUC
 740
 741
 742
 743
 744   subroutine READSNOW_FORC(flnm,ix,jx,weasd,snodep)
 745     implicit none
 746
 747     character(len=*),                   intent(in)  :: flnm
 748     integer,                            intent(in)  :: ix
 749     integer,                            intent(in)  :: jx
 750     real,             dimension(ix,jx), intent(out) :: weasd
 751     real,             dimension(ix,jx), intent(out) :: snodep
 752     real, dimension(ix,jx) :: tmp
 753
 754     character(len=256) :: units
 755     integer :: ierr
 756     integer :: ncid,i,j
 757
 758     ! Open the NetCDF file.
 759
 760     ierr = nf_open(flnm, NF_NOWRITE, ncid)
 761     if (ierr /= 0) then
 762        write(*,'("READSNOW Problem opening netcdf file: ''", A, "''")') trim(flnm)
 763        call hydro_stop("In READSNOW_FORC() - Problem opening netcdf file")
 764     endif
 765
 766     call get_2d_netcdf("WEASD",  ncid, tmp,   units, ix, jx, .FALSE., ierr)
 767     if (ierr /= 0) then
 768          call get_2d_netcdf("SNOW",  ncid, tmp,   units, ix, jx, .FALSE., ierr)
 769          if (ierr == 0) then
 770             units = "mm"
 771             print *, "read WEASD from wrfoutput ...... "
 772             weasd = tmp * 1.E-3
 773          endif
 774     else
 775          weasd = tmp
 776          if (trim(units) == "m") then
 777             ! No conversion necessary
 778          else if (trim(units) == "mm") then
 779             ! convert WEASD from mm to m
 780             weasd = weasd * 1.E-3
 781          endif
 782     endif
 783
 784     if (ierr /= 0) then
 785        print *, "!!!!! NO WEASD present in input file...initialize to 0."
 786     endif
 787
 788     call get_2d_netcdf("SNODEP",     ncid, tmp,   units, ix, jx, .FALSE., ierr)
 789     if (ierr /= 0) then
 790        ! Quick assumption regarding snow depth.
 791        call get_2d_netcdf("SNOWH",     ncid, tmp,   units, ix, jx, .FALSE., ierr)
 792        if(ierr .eq. 0) then
 793             print *, "read snow depth from wrfoutput ... "
 794             snodep = tmp
 795        endif
 796     else
 797        snodep = tmp
 798     endif
 799
 800     if (ierr /= 0) then
 801        ! Quick assumption regarding snow depth.
 802 !yw       snodep = weasd * 10.
 803        where(snodep .lt. weasd) snodep = weasd*10  !set lower bound to correct bi-lin interp err...
 804     endif
 805
 806 !DJG check for erroneous neg WEASD or SNOWD due to offline interpolation...
 807        where(snodep .lt. 0) snodep = 0
 808        where(weasd .lt. 0) weasd = 0
 809     ierr = nf_close(ncid)
 810
 811   end subroutine READSNOW_FORC
 812
 813     subroutine get2d_hrldas(inflnm,ix,jx,nsoil,smc,stc,sh2ox,cmc,t1,weasd,snodep)
 814           implicit none
 815           integer :: iret,varid,ncid,ix,jx,nsoil,ierr
 816           real,dimension(ix,jx):: weasd,snodep,cmc,t1
 817           real,dimension(ix,jx,nsoil):: smc,stc,sh2ox
 818           character(len=*), intent(in) :: inflnm
 819           character(len=256)::   units
 820           iret = nf_open(trim(inflnm), NF_NOWRITE, ncid)
 821           if(iret .ne. 0 )then
 822               write(6,*) "Error: failed to open file :",trim(inflnm)
 823              call hydro_stop("In get2d_hrldas() - failed to open file")
 824           endif
 825
 826           call get2d_hrldas_real("CMC",     ncid, cmc,     ix, jx)
 827           call get2d_hrldas_real("TSKIN",   ncid, t1,      ix, jx)
 828           call get2d_hrldas_real("SWE",      ncid, weasd,   ix, jx)
 829           call get2d_hrldas_real("SNODEP",     ncid, snodep,   ix, jx)
 830
 831     call get2d_hrldas_real("SOIL_T_1",    ncid, stc(:,:,1),  ix, jx)
 832     call get2d_hrldas_real("SOIL_T_2",    ncid, stc(:,:,2),  ix, jx)
 833     call get2d_hrldas_real("SOIL_T_3",    ncid, stc(:,:,3),  ix, jx)
 834     call get2d_hrldas_real("SOIL_T_4",    ncid, stc(:,:,4),  ix, jx)
 835     call get2d_hrldas_real("SOIL_T_5",    ncid, stc(:,:,5),  ix, jx)
 836     call get2d_hrldas_real("SOIL_T_6",    ncid, stc(:,:,6),  ix, jx)
 837     call get2d_hrldas_real("SOIL_T_7",    ncid, stc(:,:,7),  ix, jx)
 838     call get2d_hrldas_real("SOIL_T_8",    ncid, stc(:,:,8),  ix, jx)
 839
 840     call get2d_hrldas_real("SOIL_M_1",    ncid, SMC(:,:,1),  ix, jx)
 841     call get2d_hrldas_real("SOIL_M_2",    ncid, SMC(:,:,2),  ix, jx)
 842     call get2d_hrldas_real("SOIL_M_3",    ncid, SMC(:,:,3),  ix, jx)
 843     call get2d_hrldas_real("SOIL_M_4",    ncid, SMC(:,:,4),  ix, jx)
 844     call get2d_hrldas_real("SOIL_M_5",    ncid, SMC(:,:,5),  ix, jx)
 845     call get2d_hrldas_real("SOIL_M_6",    ncid, SMC(:,:,6),  ix, jx)
 846     call get2d_hrldas_real("SOIL_M_7",    ncid, SMC(:,:,7),  ix, jx)
 847     call get2d_hrldas_real("SOIL_M_8",    ncid, SMC(:,:,8),  ix, jx)
 848
 849     call get2d_hrldas_real("SOIL_W_1",    ncid, SH2OX(:,:,1),  ix, jx)
 850     call get2d_hrldas_real("SOIL_W_2",    ncid, SH2OX(:,:,2),  ix, jx)
 851     call get2d_hrldas_real("SOIL_W_3",    ncid, SH2OX(:,:,3),  ix, jx)
 852     call get2d_hrldas_real("SOIL_W_4",    ncid, SH2OX(:,:,4),  ix, jx)
 853     call get2d_hrldas_real("SOIL_W_5",    ncid, SH2OX(:,:,5),  ix, jx)
 854     call get2d_hrldas_real("SOIL_W_6",    ncid, SH2OX(:,:,6),  ix, jx)
 855     call get2d_hrldas_real("SOIL_W_7",    ncid, SH2OX(:,:,7),  ix, jx)
 856     call get2d_hrldas_real("SOIL_W_8",    ncid, SH2OX(:,:,8),  ix, jx)
 857
 858           iret = nf_close(ncid)
 859          return
 860       end subroutine get2d_hrldas
 861
 862       subroutine get2d_hrldas_real(var_name,ncid,out_buff,ix,jx)
 863           implicit none
 864           integer ::iret,varid,ncid,ix,jx
 865           real out_buff(ix,jx)
 866           character(len=*), intent(in) :: var_name
 867           iret = nf_inq_varid(ncid,trim(var_name),  varid)
 868           iret = nf_get_var_real(ncid, varid, out_buff)
 869          return
 870       end subroutine get2d_hrldas_real
 871
 872     subroutine read_stage4(flnm,IX,JX,pcp)
 873         integer IX,JX,ierr,ncid,i,j
 874         real pcp(IX,JX),buf(ix,jx)
 875         character(len=*),  intent(in)  :: flnm
 876         character(len=256) :: units
 877
 878         ierr = nf_open(flnm, NF_NOWRITE, ncid)
 879
 880         if(ierr .ne. 0) then
 881             call hydro_stop("In read_stage4() - failed to open stage4 file.")
 882         endif
 883
 884         call get_2d_netcdf("RAINRATE",ncid, buf,   units, ix, jx, .TRUE., ierr)
 885         ierr = nf_close(ncid)
 886         do j = 1, jx
 887         do i = 1, ix
 888             if(buf(i,j) .lt. 0) then
 889                  buf(i,j) = pcp(i,j)
 890             end if
 891         end do
 892         end do
 893         pcp = buf
 894         return
 895     END subroutine read_stage4
 896
 897
 898
 899
 900  subroutine read_hydro_forcing_seq( &
 901        indir,olddate,hgrid, &
 902        ix,jx,forc_typ,snow_assim,  &
 903        T2,q2x,u,v,pres,xlong,short,prcp1,&
 904        lai,fpar,snodep,dt,k,prcp_old)
 905 ! This subrouting is going to read different forcing.
 906    implicit none
 907    ! in variable
 908    character(len=*) :: olddate,hgrid,indir
 909    character(len=256) :: filename
 910    integer :: ix,jx,forc_typ,k,snow_assim  ! k is time loop
 911    real,dimension(ix,jx):: T2,q2x,u,v,pres,xlong,short,prcp1,&
 912           prcpnew,weasd,snodep,prcp0,prcp2,prcp_old
 913    real ::  dt, wrf_dt
 914    ! tmp variable
 915    character(len=256) :: inflnm, inflnm2, product
 916    integer  :: i,j,mmflag,ierr_flg
 917    real,dimension(ix,jx):: lai,fpar
 918    character(len=4) nwxst_t
 919    logical :: fexist
 920
 921         inflnm = trim(indir)//"/"//&
 922              olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
 923              ".LDASIN_DOMAIN"//hgrid
 924
 925 !!!DJG... Call READFORC_(variable) Subroutine for forcing data...
 926 !!!DJG HRLDAS Format Forcing with hour format filename (NOTE: precip must be in mm/s!!!)
 927    if(FORC_TYP.eq.1) then
 928 !!Create forcing data filename...
 929         call geth_newdate(out_date,olddate,nint(dt))
 930         inflnm = trim(indir)//"/"//&
 931              out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
 932              ".LDASIN_DOMAIN"//hgrid
 933
 934         inquire (file=trim(inflnm), exist=fexist)
 935         if ( .not. fexist ) then
 936            print*, "no forcing data found", inflnm
 937            call hydro_stop("In read_hydro_forcing_seq")
 938         endif
 939
 940       CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
 941           PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
 942    end if
 943
 944
 945
 946
 947 !!!DJG HRLDAS Forcing with minute format filename (NOTE: precip must be in mm/s!!!)
 948    if(FORC_TYP.eq.2) then
 949 !!Create forcing data filename...
 950         call geth_newdate(out_date,olddate,nint(dt))
 951         inflnm = trim(indir)//"/"//&
 952              out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
 953              out_date(15:16)//".LDASIN_DOMAIN"//hgrid
 954         inquire (file=trim(inflnm), exist=fexist)
 955         if ( .not. fexist ) then
 956            print*, "no forcing data found", inflnm
 957            call hydro_stop("In read_hydro_forcing_seq() - no forcing data found")
 958         endif
 959       CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
 960           PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
 961    end if
 962
 963
 964
 965
 966
 967 !!!DJG WRF Output File Direct Ingest Forcing...
 968      if(FORC_TYP.eq.3) then
 969 !!Create forcing data filename...
 970         inflnm = trim(indir)//"/"//&
 971              "wrfout_d0"//hgrid//"_"//&
 972              olddate(1:4)//"-"//olddate(6:7)//"-"//olddate(9:10)//&
 973              "_"//olddate(12:13)//":00:00"
 974
 975         inquire (file=trim(inflnm), exist=fexist)
 976         if ( .not. fexist ) then
 977            print*, "no forcing data found", inflnm
 978            call hydro_stop("In read_hydro_forcing_seq() - no forcing data found")
 979         endif
 980
 981         do i_forcing = 1, int(24*3600/dt)
 982            wrf_dt = i_forcing*dt
 983            call geth_newdate(out_date,olddate,nint(wrf_dt))
 984            inflnm2 = trim(indir)//"/"//&
 985              "wrfout_d0"//hgrid//"_"//&
 986              out_date(1:4)//"-"//out_date(6:7)//"-"//out_date(9:10)//&
 987              "_"//out_date(12:13)//":00:00"
 988            inquire (file=trim(inflnm2), exist=fexist)
 989            if (fexist ) goto 991
 990         end do
 991 991     continue
 992
 993         if(.not. fexist) then
 994            write(6,*) "FATAL ERROR: could not find file ",trim(inflnm2)
 995            call hydro_stop("In read_hydro_forcing_seq() - could not find file ")
 996         endif
 997 #ifdef HYDRO_D
 998            print*, "read WRF forcing data: ", trim(inflnm)
 999            print*, "read WRF forcing data: ", trim(inflnm2)
1000 #endif
1001        CALL READFORC_WRF(inflnm2,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1002           PRES,XLONG,SHORT,PRCPnew,lai,fpar)
1003        CALL READFORC_WRF(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1004           PRES,XLONG,SHORT,prcp0,lai,fpar)
1005         PRCP1=(PRCPnew-prcp0)/wrf_dt   !Adjustment to convert accum to rate...(mm/s)
1006
1007      end if
1008
1009 !!!DJG CONSTant, idealized forcing...
1010      if(FORC_TYP.eq.4) then
1011 ! Impose a fixed diurnal cycle...
1012 ! assumes model timestep is 1 hr
1013 ! assumes K=1 is 12z (Ks or ~ sunrise)
1014 ! First Precip...
1015        IF (K.EQ.2) THEN
1016        PRCP1 =25.4/3600.0      !units mm/s  (Simulates 1"/hr for first time step...)
1017 !       PRCP1 =0./3600.0      !units mm/s  (Simulates <1"/hr for first 10 hours...)
1018        ELSEIF (K.GT.1) THEN
1019 !        PRCP1 =0./3600.0      !units mm/s
1020 !       ELSE
1021          PRCP1 = 0.
1022        END IF
1023 !       PRCP1 = 0.
1024 !       PRCP1 =10./3600.0      !units mm/s
1025 ! Other Met. Vars...
1026        T2=290.0 + 3.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1027        Q2X = 0.01
1028        U = 1.0
1029        V = 1.0
1030        PRES = 100000.0
1031        XLONG=400.0 + 25.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1032        SHORT=450.0 + 450.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1033
1034 !      print *, "PCP", PRCP1
1035
1036     end if
1037
1038 !!!DJG  Idealized Met. w/ Specified Precip. Forcing Data...(Note: input precip units here are in 'mm/hr')
1039 !   This option uses hard-wired met forcing EXCEPT precipitation which is read in
1040 !   from a single, separate input file called 'YYYYMMDDHHMM.PRECIP_FORCING.nc'
1041 !
1042     if(FORC_TYP.eq.5) then
1043 ! Standard Met. Vars...
1044        T2=290.0 + 3.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1045        Q2X = 0.01
1046        U = 1.0
1047        V = 1.0
1048        PRES = 100000.0
1049        XLONG=400.0 + 25.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1050        SHORT=450.0 + 450.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
1051
1052 !Get specified precip....
1053 !!!VIP, dimensions of grid are currently hardwired in input subroutine!!!
1054 !       product = "trmm"
1055 !       inflnm = trim(indir)//"/"//"sat_domain1.nc"
1056 !!Create forcing data filename...
1057         inflnm = trim(indir)//"/"//&
1058                 olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1059                 olddate(15:16)//".PRECIP_FORCING.nc"
1060         inquire (file=trim(inflnm), exist=fexist)
1061         if ( .not. fexist ) then
1062            print*, "no specified precipitation data found", inflnm
1063            call hydro_stop("In read_hydro_forcing_seq() - no specified precipitation data found")
1064         endif
1065
1066        PRCP1 = 0.
1067        PRCP_old = PRCP1
1068
1069 #ifdef HYDRO_D
1070       print *, "Opening supplemental precipitation forcing file...",inflnm
1071 #endif
1072        CALL READFORC_MDV(inflnm,IX,JX,   &
1073           PRCP2,mmflag,ierr_flg)
1074
1075 !If radar or spec. data is ok use if not, skip to original NARR data...
1076       IF (ierr_flg.eq.0) then   ! use spec. precip
1077 !Convert units if necessary
1078         IF (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
1079            PRCP1=PRCP2/DT     !convert from mm to mm/s
1080 #ifdef HYDRO_D
1081            print*, "Supplemental pcp is accumulated pcp/dt. "
1082 #endif
1083         else
1084            PRCP1=PRCP2   !assumes PRCP2 is in mm/s
1085 #ifdef HYDRO_D
1086            print*, "Supplemental pcp is rate. "
1087 #endif
1088         END IF  ! Endif mmflag
1089       ELSE   ! either stop or default to original forcing data...
1090 #ifdef HYDRO_D
1091         print *,"Current RADAR precip data not found !!! Using previous available file..."
1092 #endif
1093         PRCP1 = PRCP_old
1094       END IF  ! Endif ierr_flg
1095
1096 ! Loop through data to screen for plausible values
1097        do i=1,ix
1098          do j=1,jx
1099            if (PRCP1(i,j).lt.0.) PRCP1(i,j)= PRCP_old(i,j)
1100            if (PRCP1(i,j).gt.0.138889) PRCP1(i,j)=0.138889  !set max pcp intens = 500 mm/h
1101          end do
1102        end do
1103
1104     end if
1105
1106
1107
1108
1109
1110 !!!DJG HRLDAS Forcing with hourly format filename with specified precipitation forcing...
1111 !   This option uses HRLDAS-formatted met forcing EXCEPT precipitation which is read in
1112 !   from a single, separate input file called 'YYYYMMDDHHMM.PRECIP_FORCING.nc'
1113
1114    if(FORC_TYP.eq.6) then
1115
1116 !!Create forcing data filename...
1117         inflnm = trim(indir)//"/"//&
1118              olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1119              ".LDASIN_DOMAIN"//hgrid
1120
1121         inquire (file=trim(inflnm), exist=fexist)
1122
1123         if ( .not. fexist ) then
1124           do i_forcing = 1, nint(12*3600/dt)
1125            call geth_newdate(out_date,olddate,nint(i_forcing*dt))
1126            inflnm = trim(indir)//"/"//&
1127               olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1128               olddate(15:16)//".LDASIN_DOMAIN"//hgrid
1129            inquire (file=trim(inflnm), exist=fexist)
1130             if(fexist) goto 201
1131           end do
1132 201       continue
1133         endif
1134
1135
1136         if ( .not. fexist ) then
1137 #ifdef HYDRO_D
1138            print*, "no ATM forcing data found at this time", inflnm
1139 #endif
1140         else
1141 #ifdef HYDRO_D
1142            print*, "reading forcing data at this time", inflnm
1143 #endif
1144
1145            CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1146                 PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
1147            PRCP_old = PRCP1  ! This assigns new precip to last precip as a fallback for missing data...
1148         endif
1149
1150
1151 !Get specified precip....
1152 !!!VIP, dimensions of grid are currently hardwired in input subroutine!!!
1153 !!Create forcing data filename...
1154         inflnm = trim(indir)//"/"//&
1155                  olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1156                  olddate(15:16)//".PRECIP_FORCING.nc"
1157         inquire (file=trim(inflnm), exist=fexist)
1158 #ifdef HYDRO_D
1159         if(fexist) then
1160             print*, "using specified pcp forcing: ",trim(inflnm)
1161         else
1162             print*, "no specified pcp forcing: ",trim(inflnm)
1163         endif
1164 #endif
1165         if ( .not. fexist ) then
1166            prcp1 = PRCP_old ! for missing pcp data use analysis/model input
1167         else
1168            CALL READFORC_MDV(inflnm,IX,JX,   &
1169               PRCP2,mmflag,ierr_flg)
1170 !If radar or spec. data is ok use if not, skip to original NARR data...
1171            if(ierr_flg .ne. 0) then
1172 #ifdef HYDRO_D
1173                print*, "WARNING: pcp reading problem: ", trim(inflnm)
1174 #endif
1175                PRCP1=PRCP_old
1176            else
1177                PRCP1=PRCP2   !assumes PRCP2 is in mm/s
1178                IF (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
1179                 PRCP1=PRCP2/DT     !convert from mm to mm/s
1180                END IF  ! Endif mmflag
1181 #ifdef HYDRO_D
1182                print*, "replace pcp successfully! ",trim(inflnm)
1183 #endif
1184            endif
1185         endif
1186
1187
1188 ! Loop through data to screen for plausible values
1189        where(PRCP1 .lt. 0) PRCP1=PRCP_old
1190        where(PRCP1 .gt. 10 ) PRCP1= PRCP_old
1191        do i=1,ix
1192          do j=1,jx
1193            if (PRCP1(i,j).lt.0.) PRCP1(i,j)=0.0
1194            if (PRCP1(i,j).gt.0.138889) PRCP1(i,j)=0.138889  !set max pcp intens = 500 mm/h
1195          end do
1196        end do
1197
1198    end if
1199
1200
1201 !!!! FORC_TYP 7: uses WRF forcing data plus additional pcp forcing.
1202
1203    if(FORC_TYP.eq.7) then
1204
1205 !!Create forcing data filename...
1206         inflnm = trim(indir)//"/"//&
1207              "wrfout_d0"//hgrid//"_"//&
1208              olddate(1:4)//"-"//olddate(6:7)//"-"//olddate(9:10)//&
1209              "_"//olddate(12:13)//":00:00"
1210
1211         inquire (file=trim(inflnm), exist=fexist)
1212
1213
1214         if ( .not. fexist ) then
1215 #ifdef HYDRO_D
1216            print*, "no forcing data found", inflnm
1217 #endif
1218         else
1219            do i_forcing = 1, int(24*3600/dt)
1220               wrf_dt = i_forcing*dt
1221               call geth_newdate(out_date,olddate,nint(wrf_dt))
1222               inflnm2 = trim(indir)//"/"//&
1223                 "wrfout_d0"//hgrid//"_"//&
1224                 out_date(1:4)//"-"//out_date(6:7)//"-"//out_date(9:10)//&
1225                 "_"//out_date(12:13)//":00:00"
1226               inquire (file=trim(inflnm2), exist=fexist)
1227               if (fexist ) goto 992
1228            end do
1229 992        continue
1230
1231 #ifdef HYDRO_D
1232            print*, "read WRF forcing data: ", trim(inflnm)
1233            print*, "read WRF forcing data: ", trim(inflnm2)
1234 #endif
1235            CALL READFORC_WRF(inflnm2,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1236                    PRES,XLONG,SHORT,PRCPnew,lai,fpar)
1237            CALL READFORC_WRF(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1238                    PRES,XLONG,SHORT,prcp0,lai,fpar)
1239            PRCP1=(PRCPnew-prcp0)/wrf_dt   !Adjustment to convert accum to rate...(mm/s)
1240            PRCP_old = PRCP1
1241         endif
1242
1243 !Get specified precip....
1244 !!!VIP, dimensions of grid are currently hardwired in input subroutine!!!
1245 !!Create forcing data filename...
1246         inflnm = trim(indir)//"/"//&
1247                  olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1248                  olddate(15:16)//".PRECIP_FORCING.nc"
1249         inquire (file=trim(inflnm), exist=fexist)
1250 #ifdef HYDRO_D
1251         if(fexist) then
1252             print*, "using specified pcp forcing: ",trim(inflnm)
1253         else
1254             print*, "no specified pcp forcing: ",trim(inflnm)
1255         endif
1256 #endif
1257         if ( .not. fexist ) then
1258            prcp1 = PRCP_old ! for missing pcp data use analysis/model input
1259         else
1260            CALL READFORC_MDV(inflnm,IX,JX,   &
1261               PRCP2,mmflag,ierr_flg)
1262 !If radar or spec. data is ok use if not, skip to original NARR data...
1263            if(ierr_flg .ne. 0) then
1264 #ifdef HYDRO_D
1265                print*, "WARNING: pcp reading problem: ", trim(inflnm)
1266 #endif
1267                PRCP1=PRCP_old
1268            else
1269                PRCP1=PRCP2   !assumes PRCP2 is in mm/s
1270                IF (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
1271 #ifdef HYDRO_D
1272                  write(6,*) "using supplemental pcp time interval ", DT
1273 #endif
1274                 PRCP1=PRCP2/DT     !convert from mm to mm/s
1275                else
1276 #ifdef HYDRO_D
1277                  write(6,*) "using supplemental pcp rates "
1278 #endif
1279                END IF  ! Endif mmflag
1280 #ifdef HYDRO_D
1281                print*, "replace pcp successfully! ",trim(inflnm)
1282 #endif
1283            endif
1284         endif
1285
1286
1287 ! Loop through data to screen for plausible values
1288        where(PRCP1 .lt. 0) PRCP1=PRCP_old
1289        where(PRCP1 .gt. 10 ) PRCP1= PRCP_old ! set maximum to be 500 mm/h
1290        where(PRCP1 .gt. 0.138889) PRCP1= 0.138889 ! set maximum to be 500 mm/h
1291    end if
1292
1293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1294 !!!The other forcing data types below here are obsolete and left for reference...
1295 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1296
1297 !!!DJG HRLDAS Single Input with Multiple Input Times File Forcing...
1298 !     if(FORC_TYP.eq.6) then
1299 !!Create forcing data filename...
1300 !     if (len_trim(range) == 0) then
1301 !      inflnm = trim(indir)//"/"//&
1302 !             startdate(1:4)//startdate(6:7)//startdate(9:10)//startdate(12:13)//&
1303 !             olddate(15:16)//".LDASIN_DOMAIN"//hgrid//"_multiple"
1304 !!        "MET_LIS_CRO_2D_SANTEE_LU_1KM."//&
1305 !!        ".156hrfcst.radar"
1306 !     else
1307 !     endif
1308 !     CALL READFORC_HRLDAS_mult(inflnm,IX,JX,OLDDATE,T2,Q2X,U,   &
1309 !          PRES,XLONG,SHORT,PRCP1,K)
1310 !
1311 !!       IF (K.GT.0.AND.K.LT.10) THEN
1312 !!         PRCP1 = 10.0/3600.0            ! units mm/s
1313 !!          PRCP1 = 0.254/3600.0
1314 !!       ELSE
1315 !!         PRCP1 = 0.
1316 !!       END IF
1317 !      endif
1318
1319
1320
1321 !!!!!DJG  NARR Met. w/ NARR Precip. Forcing Data...
1322 !! Assumes standard 3-hrly NARR data has been resampled to NDHMS grid...
1323 !! Assumes one 3hrly time-step per forcing data file
1324 !! Input precip units here are in 'mm' accumulated over 3 hrs...
1325 !    if(FORC_TYP.eq.7) then  !NARR Met. w/ NARR Precip.
1326 !!!Create forcing data filename...
1327 !      if (len_trim(range) == 0) then
1328 !        inflnm = trim(indir)//"/"//&
1329 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1330 !             ".LDASIN_DOMAIN"//hgrid
1331 !      else
1332 !        inflnm = trim(indir)//"/"//&
1333 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1334 !             ".LDASIN_DOMAIN"//hgrid//"."//trim(range)
1335 !      endif
1336 !      CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1337 !          PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
1338 !      PRCP1=PRCP1/(3.0*3600.0)  ! convert from 3hr accum to mm/s which is what NDHMS expects
1339 !    end if  !NARR Met. w/ NARR Precip.
1340
1341
1342
1343
1344
1345
1346 !!!!DJG  NARR Met. w/ Specified Precip. Forcing Data...
1347 !    if(FORC_TYP.eq.8) then !NARR Met. w/ Specified Precip.
1348 !
1349 !!Check to make sure if Noah time step is 3 hrs as is NARR...
1350 !
1351 !        PRCP_old = PRCP1
1352 !
1353 !     if(K.eq.1.OR.(MOD((K-1)*INT(DT),10800)).eq.0) then   !if/then 3 hr check
1354 !!!Create forcing data filename...
1355 !      if (len_trim(range) == 0) then
1356 !        inflnm = trim(indir)//"/"//&
1357 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1358 !             ".LDASIN_DOMAIN"//hgrid
1359 !!        startdate(1:4)//startdate(6:7)//startdate(9:10)//startdate(12:13)//&
1360 !!        ".48hrfcst.ncf"
1361 !      else
1362 !        inflnm = trim(indir)//"/"//&
1363 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1364 !             ".LDASIN_DOMAIN"//hgrid//"."//trim(range)
1365 !      endif
1366 !      CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1367 !          PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
1368 !!       PRCP1=PRCP1/(3.0*3600.0)     !NARR 3hrly precip product in mm
1369 !       PRCP1=PRCP1     !NAM model data in mm/s
1370 !    end if    !3 hr check
1371 !
1372 !
1373 !!Get spec. precip....
1374 !! NAM Remote sensing...
1375 !!!!VIP, dimensions of grid are currently hardwired in input subroutine!!!
1376 !!       product = "trmm"
1377 !!       inflnm = trim(indir)//"/"//"sat_domain1.nc"
1378 !!!       inflnm = trim(indir)//"/"//"sat_domain2.nc"
1379 !!       PRCP1 = 0.
1380 !!       CALL READFORC_NAMPCP(inflnm,IX,JX,   &
1381 !!          PRCP2,K,product)
1382 !!       ierr_flg = 0
1383 !!       mmflag = 0
1384 !!!Convert pcp grid to units of mm/s...
1385 !!       PRCP1=PRCP1/(3.0*3600.0)     !3hrly precip product
1386 !
1387 !!Read from filelist (NAME HE...,others)...
1388 !!        if (K.eq.1) then
1389 !!          open(unit=93,file="filelist.txt",form="formatted",status="old")
1390 !!        end if
1391 !!        read (93,*) filename
1392 !!        inflnm = trim(indir)//"/"//trim(filename)
1393 !!
1394 !!
1395 !!Front Range MDV Radar...
1396 !
1397 !!         inflnm = "/ptmp/weiyu/rt_2008/radar_obs/"//&
1398 !!             inflnm = "/d3/hydrolab/HRLDAS_forcing/FRNG_research/20080809/"//&
1399 !!              olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1400 !!              olddate(15:16)//"_radar.nc"
1401 !!              olddate(15:16)//"_chill.nc"
1402 !
1403 !!        inflnm = "/d2/hydrolab/HRLDAS/forcing/FRNG/Big_Thomp_04/"//&
1404 !!       inflnm = "/d2/hydrolab/HRLDAS/forcing/FRNG/RT_2008/radar_obs/"//&
1405 !!             inflnm = "/d3/hydrolab/HRLDAS_forcing/FRNG_research/20080809/"//&
1406 !        inflnm = trim(indir)//"/"//&
1407 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//"_"//olddate(12:13)//&
1408 !!             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1409 !!             olddate(15:16)//"00_Pcp60min.nc"
1410 !!             olddate(15:16)//"00_Pcp30min.nc"
1411 !!             olddate(15:16)//"00_30min.nc"
1412 !             olddate(15:16)//"00_Pcp5min.nc"
1413 !!              olddate(15:16)//"_chill.nc"
1414 !
1415 !!         inflnm = "/d2/hydrolab/HRLDAS/forcing/COWS/"//&
1416 !!             olddate(1:4)//olddate(6:7)//olddate(9:10)//"_"//olddate(12:13)//&
1417 !!             olddate(15:16)//"00_Pcp5min.nc"
1418 !!              olddate(15:16)//"00_5.nc"
1419 !
1420 !!         inflnm = ""     ! use this for NAM frxst runs with 30 min time-step
1421 !!
1422 !
1423 !
1424 !!        if (K.le.6) then   ! use for 30min nowcast...
1425 !!          if (K.eq.1) then
1426 !!             open(unit=94,file="start_file.txt",form="formatted",status="replace")
1427 !!!             inflnm2 = "/d2/hydrolab/HRLDAS/forcing/FRNG/RT_2008/radar_obs/"//&
1428 !!             inflnm2 = "/d3/hydrolab/HRLDAS_forcing/FRNG_research/"//&
1429 !!             olddate(1:4)//olddate(6:7)//olddate(9:10)//"_"//olddate(12:13)//&
1430 !!             olddate(15:16)//"00_"
1431 !!             close(94)
1432 !!             nwxst_t = "5"! calc minutes from timestep and convert to char...
1433 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1434 !!          end if
1435 !!          if (K.eq.2) then
1436 !!             nwxst_t = "10" ! calc minutes from timestep and convert to char...
1437 !!             open(unit=94,file="start_file.txt",form="formatted",status="old")
1438 !!             read (94,*) inflnm2
1439 !!             close(94)
1440 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1441 !!          end if
1442 !!          if (K.eq.3) then
1443 !!             nwxst_t = "15" ! calc minutes from timestep and convert to char...
1444 !!             open(unit=94,file="start_file.txt",form="formatted",status="old")
1445 !!             read (94,*) inflnm
1446 !!             close(94)
1447 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1448 !!          end if
1449 !!          if (K.eq.4) then
1450 !!             nwxst_t = "20" ! calc minutes from timestep and convert to char...
1451 !!             open(unit=94,file="start_file.txt",form="formatted",status="old")
1452 !!             read (94,*) inflnm
1453 !!             close(94)
1454 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1455 !!          end if
1456 !!          if (K.eq.5) then
1457 !!             nwxst_t = "25" ! calc minutes from timestep and convert to char...
1458 !!             open(unit=94,file="start_file.txt",form="formatted",status="old")
1459 !!             read (94,*) inflnm
1460 !!             close(94)
1461 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1462 !!          end if
1463 !!          if (K.eq.6) then
1464 !!             nwxst_t = "30" ! calc minutes from timestep and convert to char...
1465 !!             open(unit=94,file="start_file.txt",form="formatted",status="old")
1466 !!             read (94,*) inflnm
1467 !!             close(94)
1468 !!             inflnm = trim(inflnm2)//trim(nwxst_t)//".nc"
1469 !!          end if
1470 !!        else
1471 !!          inflnm = ""     ! use this for NAM frxst runs with 30 min time-step
1472 !!        end if
1473 !
1474 !!             olddate(1:4)//olddate(6:7)//olddate(9:10)//"_"//olddate(12:13)//&
1475 !!             olddate(15:16)//"00_Pcp30minMerge.nc"
1476 !
1477 !       CALL READFORC_MDV(inflnm,IX,JX,   &
1478 !          PRCP2,mmflag,ierr_flg)
1479 !
1480 !!If radar or spec. data is ok use if not, skip to original NARR data...
1481 !      IF (ierr_flg.eq.0) then   ! use spec. precip
1482 !         PRCP1=PRCP2   !assumes PRCP2 is in mm/s
1483 !!Convert units if necessary
1484 !        IF (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
1485 !          PRCP1=PRCP2/DT     !convert from mm to mm/s
1486 !        END IF  ! Endif mmflag
1487 !      ELSE   ! either stop or default to original forcing data...
1488 !        PRCP1 = PRCP_old
1489 !      END IF  ! Endif ierr_flg
1490 !
1491 !! Loop through data to screen for plausible values
1492 !       do i=1,ix
1493 !         do j=1,jx
1494 !           if (PRCP1(i,j).lt.0.) PRCP1(i,j)=0.0
1495 !           if (PRCP1(i,j).gt.0.0555) PRCP1(i,j)=0.0555  !set max pcp intens = 200 mm/h
1496 !!          PRCP1(i,j) = 0.
1497 !!          PRCP1(i,j) = 0.02   !override w/ const. precip for gw testing only...
1498 !         end do
1499 !       end do
1500 !
1501 !!        if (K.eq.1) then  ! quick dump for site specific precip...
1502 !          open(unit=94,file="Christman_accumpcp.txt",form="formatted",status="new")
1503 !        end if
1504 !
1505 !
1506 !    end if  !NARR Met. w/ Specified Precip.
1507
1508
1509
1510
1511
1512 !!!!DJG  NLDAS Met. w/ NLDAS Precip. Forcing Data...
1513 !! Assumes standard hrly NLDAS data has been resampled to NDHMS grid...
1514 !! Assumes one 1-hrly time-step per forcing data file
1515 !! Input precip units here are in 'mm' accumulated over 1 hr...
1516 !    if(FORC_TYP.eq.9) then  !NLDAS Met. w/ NLDAS Precip.
1517 !!!Create forcing data filename...
1518 !      if (len_trim(range) == 0) then
1519 !        inflnm = trim(indir)//"/"//&
1520 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1521 !!Use this for minute forcing...             olddate(15:16)//".LDASIN_DOMAIN"//hgrid
1522 !             ".LDASIN_DOMAIN"//hgrid
1523 !      else
1524 !        inflnm = trim(indir)//"/"//&
1525 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1526 !             ".LDASIN_DOMAIN"//hgrid//"."//trim(range)
1527 !      endif
1528 !      CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1529 !          PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
1530 !      PRCP1=PRCP1/(1.0*3600.0)  ! convert hourly NLDAS hourly accum pcp to mm/s which is what NDHMS expects
1531 !    end if  !NLDAS Met. w/ NLDAS Precip.
1532
1533
1534
1535
1536
1537 !!!!DJG  NARR Met. w/ DMIP Precip. & Temp. Forcing Data...
1538 !    if(FORC_TYP.eq.10) then  ! If/Then for DMIP forcing data...
1539 !!Check to make sure if Noah time step is 3 hrs as is NARR...
1540 !
1541 !     if(K.eq.1.OR.(MOD((K-1)*INT(DT),10800)).eq.0) then   !if/then 3 hr check
1542 !!!Create forcing data filename...
1543 !      if (len_trim(range) == 0) then
1544 !        inflnm = trim(indir)//"/"//&
1545 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1546 !             ".LDASIN_DOMAIN"//hgrid
1547 !!        startdate(1:4)//startdate(6:7)//startdate(9:10)//startdate(12:13)//&
1548 !!        ".48hrfcst.ncf"
1549 !      else
1550 !        inflnm = trim(indir)//"/"//&
1551 !             olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1552 !             ".LDASIN_DOMAIN"//hgrid//"."//trim(range)
1553 !      endif
1554 !      CALL READFORC_HRLDAS(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
1555 !          PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
1556 !          PRCP1=PRCP1/(3.0*3600.0)  ! convert to mm/s which is what HRLDAS expects
1557 !    end if    !3 hr check
1558 !
1559 !!Get DMIP Precip...
1560 !!       inflnm = "/d3/gochis/HRLDAS/forcing/DMIP_II/PRECIP_HRAP/precip_finished"//"/"//&
1561 !       inflnm = "/d2/hydrolab/HRLDAS/forcing/DMIP_II_AmerR/PRECIP_HRAP"//"/"//&
1562 !           "proj.xmrg"//&
1563 !           olddate(6:7)//olddate(9:10)//olddate(1:4)//olddate(12:13)//&
1564 !           "z.asc"
1565 !        PRCP1 = 0.
1566 !        CALL READFORC_DMIP(inflnm,IX,JX,PRCP1)
1567 !          PRCP1 = PRCP1 / 100.0    ! Convert from native hundreths of mm to mm
1568 !!       IF (K.LT.34) THEN
1569 !!        PRCP1 = 5.0/3600.0            ! units mm/s
1570 !!!       ELSE
1571 !!!         PRCP1 = 0.
1572 !!       END IF
1573 !
1574 !!Get DMIP Temp...
1575 !!       inflnm = "/d3/gochis/HRLDAS/forcing/DMIP_II/TEMP_HRAP/tair_finished"//"/"//&
1576 !       inflnm = "/d2/hydrolab/HRLDAS/forcing/DMIP_II_AmerR/TEMP_HRAP"//"/"//&
1577 !           "proj.tair"//&
1578 !           olddate(6:7)//olddate(9:10)//olddate(1:4)//olddate(12:13)//&
1579 !           "z.asc"
1580 !        CALL READFORC_DMIP(inflnm,IX,JX,T2)
1581 !          T2 = (5./9.)*(T2-32.0) + 273.15         !Convert from deg F to deg K
1582 !
1583 !    end if  !End if for DMIP forcing data...
1584 !
1585 !
1586 !
1587 !! : add reading forcing precipitation data
1588 !!       ywinflnm = "/ptmp/weiyu/hrldas/v2/st4"//"/"//&
1589 !!            olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
1590 !!            ".LDASIN_DOMAIN2"
1591 !!       call read_stage4(ywinflnm,IX,JX,PRCP1)
1592 !!end yw
1593 !
1594 !
1595 !!!!DJG Check for snow data assimilation...
1596
1597    if (SNOW_ASSIM .eq. 1) then
1598
1599 ! Every 24 hours, update the snow field from analyses.
1600      if(forc_typ .ne. 3 .or. forc_typ .ne. 6) then
1601          if ( OLDDATE(12:13) == "00") then
1602             CALL READSNOW_FORC(inflnm,IX,JX,WEASD,SNODEP)
1603          endif
1604      else
1605         CALL READSNOW_FORC(inflnm,IX,JX,WEASD,SNODEP)
1606      endif
1607
1608    end if
1609
1610 #ifdef PRECIP_DOUBLE
1611 #ifdef HYDRO_D
1612    print*,'PRECIP DOUBLE'
1613 #endif
1614    PRCP1 = PRCP1 * 2.0
1615 #endif
1616
1617  end subroutine read_hydro_forcing_seq
1618
1619
1620 #ifdef MPP_LAND
1621     subroutine mpp_readland_hrldas(geo_static_flnm,&
1622           ix,jx,land_cat,soil_cat,&
1623           vegtyp,soltyp,terrain,latitude,longitude,&
1624           global_nx,global_ny,SOLVEG_INITSWC)
1625     implicit none
1626     character(len=*),          intent(in)  :: geo_static_flnm
1627     integer,                   intent(in)  :: ix, jx, land_cat, soil_cat, &
1628               global_nx,global_ny,SOLVEG_INITSWC
1629     integer, dimension(ix,jx), intent(out) :: vegtyp, soltyp
1630     real,    dimension(ix,jx), intent(out) :: terrain, latitude, longitude
1631     real, dimension(global_nx,global_ny) ::g_terrain, g_latitude, g_longitude
1632     integer, dimension(global_nx,global_ny) :: g_vegtyp, g_soltyp
1633
1634     character(len=256) :: units
1635     integer :: ierr
1636     integer :: ncid,varid
1637     real, dimension(ix,jx) :: xdum
1638     integer flag ! flag = 1 from wrfsi, flag =2 from WPS.
1639      if(my_id.eq.IO_id) then
1640         CALL READLAND_HRLDAS(geo_static_flnm,global_nx,  &
1641                global_ny,LAND_CAT,SOIL_CAT,      &
1642                g_VEGTYP,g_SOLTYP,g_TERRAIN,g_LATITUDE,g_LONGITUDE, SOLVEG_INITSWC)
1643      end if
1644   ! distribute the data to computation node.
1645      call mpp_land_bcast_int1(LAND_CAT)
1646      call mpp_land_bcast_int1(SOIL_CAT)
1647      call decompose_data_int(g_VEGTYP,VEGTYP)
1648      call decompose_data_int(g_SOLTYP,SOLTYP)
1649      call decompose_data_real(g_TERRAIN,TERRAIN)
1650      call decompose_data_real(g_LATITUDE,LATITUDE)
1651      call decompose_data_real(g_LONGITUDE,LONGITUDE)
1652       return
1653       end subroutine mpp_readland_hrldas
1654
1655
1656       subroutine MPP_READSNOW_FORC(flnm,ix,jx,OLDDATE,weasd,snodep,&
1657                  global_nX, global_ny)
1658         implicit none
1659
1660         character(len=*),                   intent(in)  :: flnm,OLDDATE
1661         integer,  intent(in)  :: ix, global_nx,global_ny
1662         integer,                            intent(in)  :: jx
1663         real,             dimension(ix,jx), intent(out) :: weasd
1664         real,             dimension(ix,jx), intent(out) :: snodep
1665
1666         real,dimension(global_nX, global_ny):: g_weasd, g_snodep
1667
1668         character(len=256) :: units
1669         integer :: ierr
1670         integer :: ncid,i,j
1671
1672         if(my_id .eq. IO_id) then
1673           CALL READSNOW_FORC(trim(flnm),global_nX, global_ny,g_WEASD,g_SNODEP)
1674        endif
1675        call decompose_data_real(g_WEASD,WEASD)
1676        call decompose_data_real(g_SNODEP,SNODEP)
1677
1678         return
1679         end  subroutine MPP_READSNOW_FORC
1680
1681       subroutine MPP_DEEPGW_HRLDAS(ix,jx,in_SMCMAX,&
1682                  global_nX, global_ny,nsoil,out_SMC,out_SH2OX)
1683         implicit none
1684
1685         integer,  intent(in)  :: ix,global_nx,global_ny
1686         integer,  intent(in)  :: jx,nsoil
1687         real,             dimension(ix,jx), intent(in) :: in_smcmax
1688         real,             dimension(ix,jx,nsoil), intent(out) :: out_smc,out_sh2ox
1689
1690         real,dimension(global_nX, global_ny,nsoil):: g_smc, g_sh2ox
1691         real,dimension(global_nX, global_ny):: g_smcmax
1692         integer   :: i,j,k
1693
1694
1695           call write_IO_real(in_smcmax,g_smcmax)  ! get global grid of smcmax
1696
1697 #ifdef HYDRO_D
1698           write (*,*) "In deep GW...", nsoil
1699 #endif
1700
1701 !loop to overwrite soils to saturation...
1702         do i=1,global_nx
1703          do j=1,global_ny
1704             g_smc(i,j,1:NSOIL) = g_smcmax(i,j)
1705             g_sh2ox(i,j,1:NSOIL) = g_smcmax(i,j)
1706          end do
1707         end do
1708
1709 !decompose global grid to parallel tiles...
1710        do k=1,nsoil
1711         call decompose_data_real(g_smc(:,:,k),out_smc(:,:,k))
1712         call decompose_data_real(g_sh2ox(:,:,k),out_sh2ox(:,:,k))
1713        end do
1714
1715         return
1716         end  subroutine MPP_DEEPGW_HRLDAS
1717
1718
1719  subroutine read_hydro_forcing_mpp( &
1720        indir,olddate,hgrid, &
1721        ix,jx,forc_typ,snow_assim,  &
1722        T2,q2x,u,v,pres,xlong,short,prcp1,&
1723        lai,fpar,snodep,dt,k,prcp_old)
1724 ! This subrouting is going to read different forcing.
1725
1726
1727    implicit none
1728    ! in variable
1729    character(len=*) :: olddate,hgrid,indir
1730    character(len=256) :: filename
1731    integer :: ix,jx,forc_typ,k,snow_assim  ! k is time loop
1732    real,dimension(ix,jx):: T2,q2x,u,v,pres,xlong,short,prcp1,&
1733           prcpnew,lai,fpar,snodep,prcp_old
1734    real ::  dt
1735    ! tmp variable
1736    character(len=256) :: inflnm, product
1737    integer  :: i,j,mmflag
1738    real,dimension(global_nx,global_ny):: g_T2,g_Q2X,g_U,g_V,g_XLONG, &
1739              g_SHORT,g_PRCP1,g_PRES,g_lai,g_snodep,g_prcp_old, g_fpar
1740    integer flag
1741
1742
1743
1744      call write_io_real(T2,g_T2)
1745      call write_io_real(Q2X,g_Q2X)
1746      call write_io_real(U,g_U)
1747      call write_io_real(V,g_V)
1748      call write_io_real(XLONG,g_XLONG)
1749      call write_io_real(SHORT,g_SHORT)
1750      call write_io_real(PRCP1,g_PRCP1)
1751      call write_io_real(PRES,g_PRES)
1752      call write_io_real(prcp_old,g_PRCP_old)
1753
1754      call write_io_real(lai,g_lai)
1755      call write_io_real(fpar,g_fpar)
1756      call write_io_real(snodep,g_snodep)
1757
1758
1759
1760    if(my_id .eq. IO_id) then
1761       call read_hydro_forcing_seq( &
1762         indir,olddate,hgrid,&
1763         global_nx,global_ny,forc_typ,snow_assim,  &
1764         g_T2,g_q2x,g_u,g_v,g_pres,g_xlong,g_short,g_prcp1,&
1765         g_lai,g_fpar,g_snodep,dt,k,g_prcp_old)
1766 #ifdef HYDRO_D
1767      write(6,*) "finish read forcing,olddate ",olddate
1768 #endif
1769    end if
1770
1771      call decompose_data_real(g_T2,T2)
1772      call decompose_data_real(g_Q2X,Q2X)
1773      call decompose_data_real(g_U,U)
1774      call decompose_data_real(g_V,V)
1775      call decompose_data_real(g_XLONG,XLONG)
1776      call decompose_data_real(g_SHORT,SHORT)
1777      call decompose_data_real(g_PRCP1,PRCP1)
1778      call decompose_data_real(g_prcp_old,prcp_old)
1779      call decompose_data_real(g_PRES,PRES)
1780
1781      call decompose_data_real(g_lai,lai)
1782      call decompose_data_real(g_fpar,fpar)
1783      call decompose_data_real(g_snodep,snodep)
1784
1785      return
1786    end subroutine read_hydro_forcing_mpp
1787 #endif
1788
1789   integer function nfeb_yw(year)
1790     !
1791     ! Compute the number of days in February for the given year.
1792     !
1793     implicit none
1794     integer, intent(in) :: year ! Four-digit year
1795
1796     nfeb_yw = 28 ! By default, February has 28 days ...
1797     if (mod(year,4).eq.0) then
1798        nfeb_yw = 29  ! But every four years, it has 29 days ...
1799        if (mod(year,100).eq.0) then
1800           nfeb_yw = 28  ! Except every 100 years, when it has 28 days ...
1801           if (mod(year,400).eq.0) then
1802              nfeb_yw = 29  ! Except every 400 years, when it has 29 days ...
1803              if (mod(year,3600).eq.0) then
1804                 nfeb_yw = 28  ! Except every 3600 years, when it has 28 days.
1805              endif
1806           endif
1807        endif
1808     endif
1809   end function nfeb_yw
1810
1811   subroutine geth_newdate (ndate, odate, idt)
1812     implicit none
1813
1814     !  From old date ("YYYY-MM-DD HH:MM:SS.ffff" or "YYYYMMDDHHMMSSffff") and
1815     !  delta-time, compute the new date.
1816
1817     !  on entry     -  odate  -  the old hdate.
1818     !                  idt    -  the change in time
1819
1820     !  on exit      -  ndate  -  the new hdate.
1821
1822     integer, intent(in)           :: idt
1823     character (len=*), intent(out) :: ndate
1824     character (len=*), intent(in)  :: odate
1825
1826     !  Local Variables
1827
1828     !  yrold    -  indicates the year associated with "odate"
1829     !  moold    -  indicates the month associated with "odate"
1830     !  dyold    -  indicates the day associated with "odate"
1831     !  hrold    -  indicates the hour associated with "odate"
1832     !  miold    -  indicates the minute associated with "odate"
1833     !  scold    -  indicates the second associated with "odate"
1834
1835     !  yrnew    -  indicates the year associated with "ndate"
1836     !  monew    -  indicates the month associated with "ndate"
1837     !  dynew    -  indicates the day associated with "ndate"
1838     !  hrnew    -  indicates the hour associated with "ndate"
1839     !  minew    -  indicates the minute associated with "ndate"
1840     !  scnew    -  indicates the second associated with "ndate"
1841
1842     !  mday     -  a list assigning the number of days in each month
1843
1844     !  i        -  loop counter
1845     !  nday     -  the integer number of days represented by "idt"
1846     !  nhour    -  the integer number of hours in "idt" after taking out
1847     !              all the whole days
1848     !  nmin     -  the integer number of minutes in "idt" after taking out
1849     !              all the whole days and whole hours.
1850     !  nsec     -  the integer number of minutes in "idt" after taking out
1851     !              all the whole days, whole hours, and whole minutes.
1852
1853     integer :: newlen, oldlen
1854     integer :: yrnew, monew, dynew, hrnew, minew, scnew, frnew
1855     integer :: yrold, moold, dyold, hrold, miold, scold, frold
1856     integer :: nday, nhour, nmin, nsec, nfrac, i, ifrc
1857     logical :: opass
1858     character (len=10) :: hfrc
1859     character (len=1) :: sp
1860     logical :: punct
1861     integer :: yrstart, yrend, mostart, moend, dystart, dyend
1862     integer :: hrstart, hrend, mistart, miend, scstart, scend, frstart
1863     integer :: units
1864     integer, dimension(12) :: mday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
1865 !yw    integer nfeb_yw
1866
1867     ! Determine if odate is "YYYY-MM-DD_HH ... " or "YYYYMMDDHH...."
1868     if (odate(5:5) == "-") then
1869        punct = .TRUE.
1870     else
1871        punct = .FALSE.
1872     endif
1873
1874     !  Break down old hdate into parts
1875
1876     hrold = 0
1877     miold = 0
1878     scold = 0
1879     frold = 0
1880     oldlen = LEN(odate)
1881     if (punct) then
1882        yrstart = 1
1883        yrend = 4
1884        mostart = 6
1885        moend = 7
1886        dystart = 9
1887        dyend = 10
1888        hrstart = 12
1889        hrend = 13
1890        mistart = 15
1891        miend = 16
1892        scstart = 18
1893        scend = 19
1894        frstart = 21
1895        select case (oldlen)
1896        case (10)
1897           ! Days
1898           units = 1
1899        case (13)
1900           ! Hours
1901           units = 2
1902        case (16)
1903           ! Minutes
1904           units = 3
1905        case (19)
1906           ! Seconds
1907           units = 4
1908        case (21)
1909           ! Tenths
1910           units = 5
1911        case (22)
1912           ! Hundredths
1913           units = 6
1914        case (23)
1915           ! Thousandths
1916           units = 7
1917        case (24)
1918           ! Ten thousandths
1919           units = 8
1920        case default
1921           write(*,*) 'ERROR: geth_newdate:  odd length: #'//trim(odate)//'#'
1922           call hydro_stop("In geth_newdate() - error odd length")
1923        end select
1924
1925        if (oldlen.ge.11) then
1926           sp = odate(11:11)
1927        else
1928           sp = ' '
1929        end if
1930
1931     else
1932
1933        yrstart = 1
1934        yrend = 4
1935        mostart = 5
1936        moend = 6
1937        dystart = 7
1938        dyend = 8
1939        hrstart = 9
1940        hrend = 10
1941        mistart = 11
1942        miend = 12
1943        scstart = 13
1944        scend = 14
1945        frstart = 15
1946
1947        select case (oldlen)
1948        case (8)
1949           ! Days
1950           units = 1
1951        case (10)
1952           ! Hours
1953           units = 2
1954        case (12)
1955           ! Minutes
1956           units = 3
1957        case (14)
1958           ! Seconds
1959           units = 4
1960        case (15)
1961           ! Tenths
1962           units = 5
1963        case (16)
1964           ! Hundredths
1965           units = 6
1966        case (17)
1967           ! Thousandths
1968           units = 7
1969        case (18)
1970           ! Ten thousandths
1971           units = 8
1972        case default
1973           write(*,*) 'ERROR: geth_newdate:  odd length: #'//trim(odate)//'#'
1974            call hydro_stop("In geth_newdate() - error odd length")
1975        end select
1976     endif
1977
1978     !  Use internal READ statements to convert the CHARACTER string
1979     !  date into INTEGER components.
1980
1981     read(odate(yrstart:yrend),  '(i4)') yrold
1982     read(odate(mostart:moend),  '(i2)') moold
1983     read(odate(dystart:dyend), '(i2)') dyold
1984     if (units.ge.2) then
1985        read(odate(hrstart:hrend),'(i2)') hrold
1986        if (units.ge.3) then
1987           read(odate(mistart:miend),'(i2)') miold
1988           if (units.ge.4) then
1989              read(odate(scstart:scend),'(i2)') scold
1990              if (units.ge.5) then
1991                 read(odate(frstart:oldlen),*) frold
1992              end if
1993           end if
1994        end if
1995     end if
1996
1997     !  Set the number of days in February for that year.
1998
1999     mday(2) = nfeb_yw(yrold)
2000
2001     !  Check that ODATE makes sense.
2002
2003     opass = .TRUE.
2004
2005     !  Check that the month of ODATE makes sense.
2006
2007     if ((moold.gt.12).or.(moold.lt.1)) then
2008 #ifdef HYDRO_D
2009        write(*,*) 'GETH_NEWDATE:  Month of ODATE = ', moold
2010 #endif
2011        opass = .FALSE.
2012     end if
2013
2014     !  Check that the day of ODATE makes sense.
2015
2016     if ((dyold.gt.mday(moold)).or.(dyold.lt.1)) then
2017 #ifdef HYDRO_D
2018        write(*,*) 'GETH_NEWDATE:  Day of ODATE = ', dyold
2019 #endif
2020        opass = .FALSE.
2021     end if
2022
2023     !  Check that the hour of ODATE makes sense.
2024
2025     if ((hrold.gt.23).or.(hrold.lt.0)) then
2026 #ifdef HYDRO_D
2027        write(*,*) 'GETH_NEWDATE:  Hour of ODATE = ', hrold
2028 #endif
2029        opass = .FALSE.
2030     end if
2031
2032     !  Check that the minute of ODATE makes sense.
2033
2034     if ((miold.gt.59).or.(miold.lt.0)) then
2035 #ifdef HYDRO_D
2036        write(*,*) 'GETH_NEWDATE:  Minute of ODATE = ', miold
2037 #endif
2038        opass = .FALSE.
2039     end if
2040
2041     !  Check that the second of ODATE makes sense.
2042
2043     if ((scold.gt.59).or.(scold.lt.0)) then
2044 #ifdef HYDRO_D
2045        write(*,*) 'GETH_NEWDATE:  Second of ODATE = ', scold
2046 #endif
2047        opass = .FALSE.
2048     end if
2049
2050     !  Check that the fractional part  of ODATE makes sense.
2051     if (.not.opass) then
2052        write(*,*) 'Crazy ODATE: ', odate(1:oldlen), oldlen
2053        call hydro_stop("In geth_newdate")
2054     end if
2055
2056     !  Date Checks are completed.  Continue.
2057
2058
2059     !  Compute the number of days, hours, minutes, and seconds in idt
2060
2061     if (units.ge.5) then !idt should be in fractions of seconds
2062        ifrc = oldlen-(frstart)+1
2063        ifrc = 10**ifrc
2064        nday   = abs(idt)/(86400*ifrc)
2065        nhour  = mod(abs(idt),86400*ifrc)/(3600*ifrc)
2066        nmin   = mod(abs(idt),3600*ifrc)/(60*ifrc)
2067        nsec   = mod(abs(idt),60*ifrc)/(ifrc)
2068        nfrac = mod(abs(idt), ifrc)
2069     else if (units.eq.4) then  !idt should be in seconds
2070        ifrc = 1
2071        nday   = abs(idt)/86400 ! integer number of days in delta-time
2072        nhour  = mod(abs(idt),86400)/3600
2073        nmin   = mod(abs(idt),3600)/60
2074        nsec   = mod(abs(idt),60)
2075        nfrac  = 0
2076     else if (units.eq.3) then !idt should be in minutes
2077        ifrc = 1
2078        nday   = abs(idt)/1440 ! integer number of days in delta-time
2079        nhour  = mod(abs(idt),1440)/60
2080        nmin   = mod(abs(idt),60)
2081        nsec   = 0
2082        nfrac  = 0
2083     else if (units.eq.2) then !idt should be in hours
2084        ifrc = 1
2085        nday   = abs(idt)/24 ! integer number of days in delta-time
2086        nhour  = mod(abs(idt),24)
2087        nmin   = 0
2088        nsec   = 0
2089        nfrac  = 0
2090     else if (units.eq.1) then !idt should be in days
2091        ifrc = 1
2092        nday   = abs(idt)    ! integer number of days in delta-time
2093        nhour  = 0
2094        nmin   = 0
2095        nsec   = 0
2096        nfrac  = 0
2097     else
2098        write(*,'(''GETH_NEWDATE: Strange length for ODATE: '', i3)') &
2099             oldlen
2100        write(*,*) '#'//odate(1:oldlen)//'#'
2101        call hydro_stop("In geth_newdate")
2102     end if
2103
2104     if (idt.ge.0) then
2105
2106        frnew = frold + nfrac
2107        if (frnew.ge.ifrc) then
2108           frnew = frnew - ifrc
2109           nsec = nsec + 1
2110        end if
2111
2112        scnew = scold + nsec
2113        if (scnew .ge. 60) then
2114           scnew = scnew - 60
2115           nmin  = nmin + 1
2116        end if
2117
2118        minew = miold + nmin
2119        if (minew .ge. 60) then
2120           minew = minew - 60
2121           nhour  = nhour + 1
2122        end if
2123
2124        hrnew = hrold + nhour
2125        if (hrnew .ge. 24) then
2126           hrnew = hrnew - 24
2127           nday  = nday + 1
2128        end if
2129
2130        dynew = dyold
2131        monew = moold
2132        yrnew = yrold
2133        do i = 1, nday
2134           dynew = dynew + 1
2135           if (dynew.gt.mday(monew)) then
2136              dynew = dynew - mday(monew)
2137              monew = monew + 1
2138              if (monew .gt. 12) then
2139                 monew = 1
2140                 yrnew = yrnew + 1
2141                 ! If the year changes, recompute the number of days in February
2142                 mday(2) = nfeb_yw(yrnew)
2143              end if
2144           end if
2145        end do
2146
2147     else if (idt.lt.0) then
2148
2149        frnew = frold - nfrac
2150        if (frnew .lt. 0) then
2151           frnew = frnew + ifrc
2152           nsec = nsec + 1
2153        end if
2154
2155        scnew = scold - nsec
2156        if (scnew .lt. 00) then
2157           scnew = scnew + 60
2158           nmin  = nmin + 1
2159        end if
2160
2161        minew = miold - nmin
2162        if (minew .lt. 00) then
2163           minew = minew + 60
2164           nhour  = nhour + 1
2165        end if
2166
2167        hrnew = hrold - nhour
2168        if (hrnew .lt. 00) then
2169           hrnew = hrnew + 24
2170           nday  = nday + 1
2171        end if
2172
2173        dynew = dyold
2174        monew = moold
2175        yrnew = yrold
2176        do i = 1, nday
2177           dynew = dynew - 1
2178           if (dynew.eq.0) then
2179              monew = monew - 1
2180              if (monew.eq.0) then
2181                 monew = 12
2182                 yrnew = yrnew - 1
2183                 ! If the year changes, recompute the number of days in February
2184                 mday(2) = nfeb_yw(yrnew)
2185              end if
2186              dynew = mday(monew)
2187           end if
2188        end do
2189     end if
2190
2191     !  Now construct the new mdate
2192
2193     newlen = LEN(ndate)
2194
2195     if (punct) then
2196
2197        if (newlen.gt.frstart) then
2198           write(ndate(1:scend),19) yrnew, monew, dynew, hrnew, minew, scnew
2199           write(hfrc,'(i10)') frnew+1000000000
2200           ndate = ndate(1:scend)//'.'//hfrc(31-newlen:10)
2201
2202        else if (newlen.eq.scend) then
2203           write(ndate(1:scend),19) yrnew, monew, dynew, hrnew, minew, scnew
2204 19        format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2,':',i2.2)
2205
2206        else if (newlen.eq.miend) then
2207           write(ndate,16) yrnew, monew, dynew, hrnew, minew
2208 16        format(i4,'-',i2.2,'-',i2.2,'_',i2.2,':',i2.2)
2209
2210        else if (newlen.eq.hrend) then
2211           write(ndate,13) yrnew, monew, dynew, hrnew
2212 13        format(i4,'-',i2.2,'-',i2.2,'_',i2.2)
2213
2214        else if (newlen.eq.dyend) then
2215           write(ndate,10) yrnew, monew, dynew
2216 10        format(i4,'-',i2.2,'-',i2.2)
2217
2218        end if
2219
2220     else
2221
2222        if (newlen.gt.frstart) then
2223           write(ndate(1:scend),119) yrnew, monew, dynew, hrnew, minew, scnew
2224           write(hfrc,'(i10)') frnew+1000000000
2225           ndate = ndate(1:scend)//'.'//hfrc(31-newlen:10)
2226
2227        else if (newlen.eq.scend) then
2228           write(ndate(1:scend),119) yrnew, monew, dynew, hrnew, minew, scnew
2229 119       format(i4,i2.2,i2.2,i2.2,i2.2,i2.2)
2230
2231        else if (newlen.eq.miend) then
2232           write(ndate,116) yrnew, monew, dynew, hrnew, minew
2233 116       format(i4,i2.2,i2.2,i2.2,i2.2)
2234
2235        else if (newlen.eq.hrend) then
2236           write(ndate,113) yrnew, monew, dynew, hrnew
2237 113       format(i4,i2.2,i2.2,i2.2)
2238
2239        else if (newlen.eq.dyend) then
2240           write(ndate,110) yrnew, monew, dynew
2241 110       format(i4,i2.2,i2.2)
2242
2243        end if
2244
2245     endif
2246
2247     if (punct .and. (oldlen.ge.11) .and. (newlen.ge.11)) ndate(11:11) = sp
2248
2249   end subroutine geth_newdate
2250
2251
2252 subroutine read_hydro_forcing_mpp1( &
2253      indir,olddate,hgrid, &
2254      ix,jx,forc_typ,snow_assim,  &
2255      T2,q2x,u,v,pres,xlong,short,prcp1,&
2256      lai,fpar,snodep,dt,k,prcp_old)
2257 ! This subrouting is going to read different forcing.
2258 implicit none
2259 ! in variable
2260 character(len=*) :: olddate,hgrid,indir
2261 character(len=256) :: filename
2262 integer :: ix,jx,forc_typ,k,snow_assim  ! k is time loop
2263 real,dimension(ix,jx):: T2,q2x,u,v,pres,xlong,short,prcp1,&
2264      prcpnew,weasd,snodep,prcp0,prcp2,prcp_old
2265 real ::  dt, wrf_dt
2266 ! tmp variable
2267 character(len=256) :: inflnm, inflnm2, product
2268 integer  :: i,j,mmflag,ierr_flg
2269 real,dimension(ix,jx):: lai,fpar
2270 character(len=4) nwxst_t
2271 logical :: fexist
2272
2273 inflnm = trim(indir)//"/"//&
2274          olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
2275          ".LDASIN_DOMAIN"//hgrid
2276
2277 !!!DJG... Call READFORC_(variable) Subroutine for forcing data...
2278 !!!DJG HRLDAS Format Forcing with hour format filename (NOTE: precip must be in mm/s!!!)
2279
2280
2281 !!! FORC_TYPE 1 ============================================================================
2282 if(FORC_TYP.eq.1) then
2283    !!Create forcing data filename...
2284    call geth_newdate(out_date,olddate,nint(dt))
2285    inflnm = trim(indir)//"/"//&
2286         out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
2287         ".LDASIN_DOMAIN"//hgrid
2288
2289    inquire (file=trim(inflnm), exist=fexist)
2290
2291 #ifdef MPP_LAND
2292    call mpp_land_bcast_logical(fexist)
2293 #endif
2294    if ( .not. fexist ) then
2295       print*, "no forcing data found", inflnm
2296       call hydro_stop("In read_hydro_forcing_mpp1() - no forcing data found")
2297    endif
2298
2299 #ifdef HYDRO_D
2300    print*, "read forcing data at ", OLDDATE,  trim(inflnm)
2301 #endif
2302    call READFORC_HRLDAS_mpp(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2303         PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
2304
2305    where(PRCP1 .lt. 0) PRCP1= 0  ! set minimum to be 0
2306    where(PRCP1 .gt. 0.138889) PRCP1= 0.138889 ! set maximum to be 500 mm/h
2307
2308 end if
2309
2310
2311 !!! FORC_TYPE 2 ============================================================================
2312 !!!DJG HRLDAS Forcing with minute format filename (NOTE: precip must be in mm/s!!!)
2313 if(FORC_TYP.eq.2) then
2314
2315 !!Create forcing data filename...
2316    call geth_newdate(out_date,olddate,nint(dt))
2317    inflnm = trim(indir)//"/"//&
2318         out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
2319         out_date(15:16)//".LDASIN_DOMAIN"//hgrid
2320    inquire (file=trim(inflnm), exist=fexist)
2321 #ifdef MPP_LAND
2322    call mpp_land_bcast_logical(fexist)
2323 #endif
2324    if ( .not. fexist ) then
2325       print*, "no forcing data found", inflnm
2326       call hydro_stop("In read_hydro_forcing_mpp1() - no forcing data found")
2327    endif
2328    call READFORC_HRLDAS_mpp(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2329         PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
2330
2331    where(PRCP1 .lt. 0) PRCP1= 0  ! set minimum to be 0
2332    where(PRCP1 .gt. 0.138889) PRCP1= 0.138889 ! set maximum to be 500 mm/h
2333
2334 end if
2335
2336
2337 !!! FORC_TYPE 3 ============================================================================
2338 !!!DJG WRF Output File Direct Ingest Forcing...
2339 if(FORC_TYP.eq.3) then
2340    !!Create forcing data filename...
2341    inflnm = trim(indir)//"/"//&
2342         "wrfout_d0"//hgrid//"_"//&
2343         olddate(1:4)//"-"//olddate(6:7)//"-"//olddate(9:10)//&
2344         "_"//olddate(12:13)//":00:00"
2345
2346    inquire (file=trim(inflnm), exist=fexist)
2347 #ifdef MPP_LAND
2348    call mpp_land_bcast_logical(fexist)
2349 #endif
2350    if ( .not. fexist ) then
2351       print*, "no forcing data found", inflnm
2352       call hydro_stop("read_hydro_forcing_seq")
2353    endif
2354
2355    do i_forcing = 1, int(24*3600/dt)
2356       wrf_dt = i_forcing*dt
2357       call geth_newdate(out_date,olddate,nint(wrf_dt))
2358       inflnm2 = trim(indir)//"/"//&
2359            "wrfout_d0"//hgrid//"_"//&
2360            out_date(1:4)//"-"//out_date(6:7)//"-"//out_date(9:10)//&
2361            "_"//out_date(12:13)//":00:00"
2362       inquire (file=trim(inflnm2), exist=fexist)
2363 #ifdef MPP_LAND
2364       call mpp_land_bcast_logical(fexist)
2365 #endif
2366       if (fexist ) goto 991
2367    end do
2368 991 continue
2369
2370    if(.not. fexist) then
2371       write(6,*) "Error: could not find file ",trim(inflnm2)
2372       call hydro_stop("In read_hydro_forcing_mpp1() - could not find WRF forcing file")
2373    endif
2374 #ifdef HYDRO_D
2375    print*, "read WRF forcing data: ", trim(inflnm)
2376    print*, "read WRF forcing data: ", trim(inflnm2)
2377 #endif
2378
2379
2380    call READFORC_WRF_mpp(inflnm2,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2381         PRES,XLONG,SHORT,PRCPnew,lai,fpar)
2382    call READFORC_WRF_mpp(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2383         PRES,XLONG,SHORT,prcp0,lai,fpar)
2384    PRCP1=(PRCPnew-prcp0)/wrf_dt   !Adjustment to convert accum to rate...(mm/s)
2385
2386 end if
2387
2388
2389 !!! FORC_TYPE4 ============================================================================
2390 !!!DJG CONSTant, idealized forcing...
2391 if(FORC_TYP.eq.4) then
2392    ! Impose a fixed diurnal cycle...
2393    ! assumes model timestep is 1 hr
2394    ! assumes K=1 is 12z (Ks or ~ sunrise)
2395    ! First Precip...
2396    !       IF (K.GE.1 .and. K.LE.2) THEN
2397    if (K.eq.1) then
2398       PRCP1 =25.4/3600.0      !units mm/s  (Simulates 1"/hr for first time step...)
2399    elseif (K.gt.1) then
2400       PRCP1 = 0.
2401    end if
2402    ! Other Met. Vars...
2403    T2=290.0 + 3.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2404    Q2X = 0.01
2405    U = 1.0
2406    V = 1.0
2407    PRES = 100000.0
2408    XLONG=400.0 + 25.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2409    SHORT=450.0 + 450.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2410 end if
2411
2412
2413 !!! FORC_TYPE 5 ============================================================================
2414 !!!DJG  Idealized Met. w/ Specified Precip. Forcing Data...(Note: input precip units here are in 'mm/hr')
2415 !   This option uses hard-wired met forcing EXCEPT precipitation which is read in
2416 !   from a single, separate input file called 'YYYYMMDDHHMM.PRECIP_FORCING.nc'
2417 !
2418 if(FORC_TYP.eq.5) then
2419    ! Standard Met. Vars...
2420    T2=290.0 + 3.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2421    Q2X = 0.01
2422    U = 1.0
2423    V = 1.0
2424    PRES = 100000.0
2425    XLONG=400.0 + 25.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2426    SHORT=450.0 + 450.0*(cos((2*3.1416*K/24.0)-12.0*2*3.1416/24.0))
2427
2428    !Get specified precip....
2429 !!!VIP, dimensions of grid are currently hardwired in input subroutine!!!
2430    !       product = "trmm"
2431    !       inflnm = trim(indir)//"/"//"sat_domain1.nc"
2432    !!Create forcing data filename...
2433    inflnm = trim(indir)//"/"//&
2434         olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
2435         olddate(15:16)//".PRECIP_FORCING.nc"
2436    inquire (file=trim(inflnm), exist=fexist)
2437 #ifdef MPP_LAND
2438    call mpp_land_bcast_logical(fexist)
2439 #endif
2440    if ( .not. fexist ) then
2441       print*, "no specified precipitation data found", inflnm
2442       call hydro_stop("In read_hydro_forcing_mpp1() - no specified precipitation data found")
2443    endif
2444
2445    PRCP1 = 0.
2446    PRCP_old = PRCP1
2447
2448 #ifdef HYDRO_D
2449    print *, "Opening supplemental precipitation forcing file...",inflnm
2450 #endif
2451    call READFORC_MDV_mpp(inflnm,IX,JX,   &
2452         PRCP2,mmflag,ierr_flg)
2453
2454    !If radar or spec. data is ok use if not, skip to original NARR data...
2455    if (ierr_flg.eq.0) then   ! use spec. precip
2456       !Convert units if necessary
2457       if (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
2458          PRCP1=PRCP2/DT     !convert from mm to mm/s
2459 #ifdef HYDRO_D
2460          print*, "Supplemental pcp is accumulated pcp/dt. "
2461 #endif
2462       else
2463          PRCP1=PRCP2   !assumes PRCP2 is in mm/s
2464 #ifdef HYDRO_D
2465          print*, "Supplemental pcp is rate. "
2466 #endif
2467       end if  ! Endif mmflag
2468    else   ! either stop or default to original forcing data...
2469 #ifdef HYDRO_D
2470       print *,"Current RADAR precip data not found !!! Using previous available file..."
2471 #endif
2472       PRCP1 = PRCP_old
2473    end if  ! Endif ierr_flg
2474
2475    ! Loop through data to screen for plausible values
2476    do i=1,ix
2477       do j=1,jx
2478          if (PRCP1(i,j).lt.0.) PRCP1(i,j)= PRCP_old(i,j)
2479          if (PRCP1(i,j).gt.0.138889) PRCP1(i,j)=0.138889  !set max pcp intens = 500 mm/h
2480       end do
2481    end do
2482
2483 end if
2484
2485
2486 !!! FORC_TYPE 6 ============================================================================
2487 !!!DJG HRLDAS Forcing with hourly format filename with specified precipitation forcing...
2488 !   This option uses HRLDAS-formatted met forcing EXCEPT precipitation which is read in
2489 !   from a single, separate input file called 'YYYYMMDDHHMM.PRECIP_FORCING.nc'
2490
2491 if(FORC_TYP.eq.6) then
2492
2493    !!Create forcing data filename...
2494
2495 #ifdef MPP_LAND
2496    if(my_id .eq. io_id) then
2497 #endif
2498       do i_forcing = 1, nint(3600*12/dt)
2499          call geth_newdate(out_date,olddate,nint(dt*i_forcing))
2500          inflnm = trim(indir)//"/"//&
2501               out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
2502               ".LDASIN_DOMAIN"//hgrid
2503
2504          inquire (file=trim(inflnm), exist=fexist)
2505          if(fexist) goto 101
2506       enddo
2507 101   continue
2508 #ifdef MPP_LAND
2509    endif
2510    call mpp_land_bcast_logical(fexist)
2511 #endif
2512
2513    if ( .not. fexist ) then
2514 #ifdef MPP_LAND
2515       if(my_id .eq. io_id) then
2516 #endif
2517          do i_forcing = 1, nint(3600*12/dt)
2518             call geth_newdate(out_date,olddate,nint(dt*i_forcing))
2519             inflnm = trim(indir)//"/"//&
2520                  out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
2521                  out_date(15:16)//".LDASIN_DOMAIN"//hgrid
2522             inquire (file=trim(inflnm), exist=fexist)
2523             if(fexist) goto 102
2524          end do
2525 102      continue
2526 #ifdef MPP_LAND
2527       endif
2528       call mpp_land_bcast_logical(fexist)
2529 #endif
2530    endif
2531
2532
2533    if ( .not. fexist ) then
2534 #ifdef HYDRO_D
2535       print*, "no ATM forcing data found at this time", inflnm
2536 #endif
2537    else
2538 #ifdef HYDRO_D
2539       print*, "reading forcing data at this time", inflnm
2540 #endif
2541
2542       call READFORC_HRLDAS_mpp(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2543            PRES,XLONG,SHORT,PRCP1,LAI,FPAR)
2544       PRCP_old = PRCP1  ! This assigns new precip to last precip as a fallback for missing data...
2545    endif
2546
2547
2548    !Get specified precip....
2549    !VIP, dimensions of grid are currently hardwired in input subroutine!!!
2550    !!Create forcing data filename...
2551    call geth_newdate(out_date,olddate,nint(dt))
2552    inflnm = trim(indir)//"/"//&
2553         out_date(1:4)//out_date(6:7)//out_date(9:10)//out_date(12:13)//&
2554         out_date(15:16)//".PRECIP_FORCING.nc"
2555    inquire (file=trim(inflnm), exist=fexist)
2556 #ifdef MPP_LAND
2557    call mpp_land_bcast_logical(fexist)
2558 #endif
2559 #ifdef HYDRO_D
2560    if(my_id .eq. io_id) then
2561       if(fexist) then
2562          print*, "using specified pcp forcing: ",trim(inflnm)
2563       else
2564          print*, "no specified pcp forcing: ",trim(inflnm)
2565       endif
2566    endif
2567 #endif
2568    if ( .not. fexist ) then
2569       prcp1 = PRCP_old ! for missing pcp data use analysis/model input
2570    else
2571       call READFORC_MDV_mpp(inflnm,IX,JX,   &
2572            PRCP2,mmflag,ierr_flg)
2573       !If radar or spec. data is ok use if not, skip to original NARR data...
2574       if(ierr_flg .ne. 0) then
2575 #ifdef HYDRO_D
2576          print*, "WARNING: pcp reading problem: ", trim(inflnm)
2577 #endif
2578          PRCP1=PRCP_old
2579       else
2580          PRCP1=PRCP2   !assumes PRCP2 is in mm/s
2581          if (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
2582             PRCP1=PRCP2/DT     !convert from mm to mm/s
2583          end if  ! Endif mmflag
2584 #ifdef HYDRO_D
2585          if(my_id .eq. io_id) then
2586             print*, "replace pcp successfully! ",trim(inflnm)
2587          endif
2588 #endif
2589       endif
2590    endif
2591
2592
2593    ! Loop through data to screen for plausible values
2594    where(PRCP1 .lt. 0) PRCP1=PRCP_old
2595    where(PRCP1 .gt. 10 ) PRCP1= PRCP_old
2596    do i=1,ix
2597       do j=1,jx
2598          if (PRCP1(i,j).lt.0.) PRCP1(i,j)=0.0
2599          if (PRCP1(i,j).gt.0.138889) PRCP1(i,j)=0.138889  !set max pcp intens = 500 mm/h
2600       end do
2601    end do
2602    !       write(80,*) prcp1
2603
2604 end if
2605
2606
2607 !!! FORC_TYPE 7 ============================================================================
2608 !!!! FORC_TYP 7: uses WRF forcing data plus additional pcp forcing.
2609
2610 if(FORC_TYP.eq.7) then
2611
2612    !!Create forcing data filename...
2613    inflnm = trim(indir)//"/"//&
2614         "wrfout_d0"//hgrid//"_"//&
2615         olddate(1:4)//"-"//olddate(6:7)//"-"//olddate(9:10)//&
2616         "_"//olddate(12:13)//":00:00"
2617
2618    inquire (file=trim(inflnm), exist=fexist)
2619 #ifdef MPP_LAND
2620    call mpp_land_bcast_logical(fexist)
2621 #endif
2622
2623
2624    if ( .not. fexist ) then
2625 #ifdef HYDRO_D
2626       print*, "no forcing data found", inflnm
2627 #endif
2628    else
2629       do i_forcing = 1, int(24*3600/dt)
2630          wrf_dt = i_forcing*dt
2631          call geth_newdate(out_date,olddate,nint(wrf_dt))
2632          inflnm2 = trim(indir)//"/"//&
2633               "wrfout_d0"//hgrid//"_"//&
2634               out_date(1:4)//"-"//out_date(6:7)//"-"//out_date(9:10)//&
2635               "_"//out_date(12:13)//":00:00"
2636          inquire (file=trim(inflnm2), exist=fexist)
2637 #ifdef MPP_LAND
2638          call mpp_land_bcast_logical(fexist)
2639 #endif
2640          if (fexist ) goto 992
2641       end do
2642 992   continue
2643
2644 #ifdef HYDRO_D
2645       print*, "read WRF forcing data: ", trim(inflnm)
2646       print*, "read WRF forcing data: ", trim(inflnm2)
2647 #endif
2648       call READFORC_WRF_mpp(inflnm2,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2649            PRES,XLONG,SHORT,PRCPnew,lai,fpar)
2650       call READFORC_WRF_mpp(inflnm,IX,JX,OLDDATE,T2,Q2X,U,V,   &
2651            PRES,XLONG,SHORT,prcp0,lai,fpar)
2652       PRCP1=(PRCPnew-prcp0)/wrf_dt   !Adjustment to convert accum to rate...(mm/s)
2653       PRCP_old = PRCP1
2654    endif
2655
2656    !Get specified precip....
2657    !VIP, dimensions of grid are currently hardwired in input subroutine!!!
2658    !!Create forcing data filename...
2659    inflnm = trim(indir)//"/"//&
2660         olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
2661         olddate(15:16)//".PRECIP_FORCING.nc"
2662    inquire (file=trim(inflnm), exist=fexist)
2663 #ifdef MPP_LAND
2664    call mpp_land_bcast_logical(fexist)
2665 #endif
2666 #ifdef HYDRO_D
2667    if(fexist) then
2668       print*, "using specified pcp forcing: ",trim(inflnm)
2669    else
2670       print*, "no specified pcp forcing: ",trim(inflnm)
2671    endif
2672 #endif
2673    if ( .not. fexist ) then
2674       prcp1 = PRCP_old ! for missing pcp data use analysis/model input
2675    else
2676       call READFORC_MDV_mpp(inflnm,IX,JX,   &
2677            PRCP2,mmflag,ierr_flg)
2678       !If radar or spec. data is ok use if not, skip to original NARR data...
2679       if(ierr_flg .ne. 0) then
2680 #ifdef HYDRO_D
2681          print*, "WARNING: pcp reading problem: ", trim(inflnm)
2682 #endif
2683          PRCP1=PRCP_old
2684       else
2685          PRCP1=PRCP2   !assumes PRCP2 is in mm/s
2686          if (mmflag.eq.0) then    !Convert pcp grid to units of mm/s...
2687             write(6,*) "using supplemental pcp time interval ", DT
2688             PRCP1=PRCP2/DT     !convert from mm to mm/s
2689          else
2690             write(6,*) "using supplemental pcp rates "
2691          end if  ! Endif mmflag
2692 #ifdef HYDRO_D
2693          print*, "replace pcp successfully! ",trim(inflnm)
2694 #endif
2695       endif
2696    endif
2697
2698
2699    ! Loop through data to screen for plausible values
2700    where(PRCP1 .lt. 0) PRCP1=PRCP_old
2701    where(PRCP1 .gt. 10 ) PRCP1= PRCP_old ! set maximum to be 500 mm/h
2702    where(PRCP1 .gt. 0.138889) PRCP1= 0.138889 ! set maximum to be 500 mm/h
2703 end if
2704
2705
2706
2707 !!!!DJG Check for snow data assimilation...
2708 if (SNOW_ASSIM .eq. 1) then
2709
2710    ! Every 24 hours, update the snow field from analyses.
2711    if(forc_typ .ne. 3 .or. forc_typ .ne. 6) then
2712       if ( OLDDATE(12:13) == "00") then
2713          call READSNOW_FORC_mpp(inflnm,IX,JX,WEASD,SNODEP)
2714       endif
2715    else
2716       call READSNOW_FORC_mpp(inflnm,IX,JX,WEASD,SNODEP)
2717    endif
2718
2719 end if
2720
2721 #ifdef PRECIP_DOUBLE
2722 #ifdef HYDRO_D
2723 print*,'PRECIP DOUBLE'
2724 #endif
2725 PRCP1 = PRCP1 * 2.0
2726 #endif
2727
2728 end subroutine read_hydro_forcing_mpp1
2729
2730
2731
2732   subroutine READFORC_HRLDAS_mpp(flnm,ix,jx,target_date, t,q,u,v,p,lw,sw,pcp,lai,fpar)
2733     implicit none
2734
2735     character(len=*),                   intent(in)  :: flnm
2736     integer,                            intent(in)  :: ix
2737     integer,                            intent(in)  :: jx
2738     character(len=*),                   intent(in)  :: target_date
2739     real,             dimension(ix,jx), intent(out) :: t
2740     real,             dimension(ix,jx), intent(out) :: q
2741     real,             dimension(ix,jx), intent(out) :: u
2742     real,             dimension(ix,jx), intent(out) :: v
2743     real,             dimension(ix,jx), intent(out) :: p
2744     real,             dimension(ix,jx), intent(out) :: lw
2745     real,             dimension(ix,jx), intent(out) :: sw
2746     real,             dimension(ix,jx), intent(out) :: pcp
2747     real,             dimension(ix,jx), intent(inout) :: lai
2748     real,             dimension(ix,jx), intent(inout) :: fpar
2749
2750     character(len=256) :: units
2751     integer :: ierr
2752     integer :: ncid
2753
2754     ! Open the NetCDF file.
2755 #ifdef MPP_LAND
2756     real, allocatable, dimension(:,:):: buf2
2757     if(my_id .eq. io_id) then
2758         allocate(buf2(global_nx,global_ny))
2759     else
2760         allocate(buf2(1,1))
2761     endif
2762     if(my_id .eq. io_id) then
2763         ierr = nf_open(trim(flnm), NF_NOWRITE, ncid)
2764     endif
2765     call mpp_land_bcast_int1(ierr)
2766     if (ierr /= 0) then
2767        write(*,'("READFORC Problem opening netcdf file: ''", A, "''")') trim(flnm)
2768        call hydro_stop("In READFORC_HRLDAS_mpp() - Problem opening netcdf file")
2769     endif
2770
2771
2772     if(my_id .eq. io_id ) call get_2d_netcdf("T2D", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2773     call decompose_data_real (buf2,t)
2774     if(my_id .eq. io_id ) call get_2d_netcdf("Q2D", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2775     call decompose_data_real (buf2,q)
2776     if(my_id .eq. io_id ) call get_2d_netcdf("U2D", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2777     call decompose_data_real (buf2,u)
2778     if(my_id .eq. io_id ) call get_2d_netcdf("V2D", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2779     call decompose_data_real (buf2,v)
2780     if(my_id .eq. io_id ) call get_2d_netcdf("PSFC", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2781     call decompose_data_real (buf2,p)
2782     if(my_id .eq. io_id ) call get_2d_netcdf("LWDOWN", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2783     call decompose_data_real (buf2,lw)
2784     if(my_id .eq. io_id ) call get_2d_netcdf("SWDOWN", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2785     call decompose_data_real (buf2,sw)
2786     if(my_id .eq. io_id ) call get_2d_netcdf("RAINRATE", ncid,buf2, units, global_nx, global_ny, .TRUE., ierr)
2787     call decompose_data_real (buf2,pcp)
2788     if(my_id .eq. io_id ) then
2789           call get_2d_netcdf("VEGFRA", ncid,buf2, units, global_nx, global_ny, .FALSE., ierr)
2790           if (ierr == 0) then
2791             if(maxval(buf2) .gt. 10 .and. maxval(buf2) .lt. 10000)  buf2 = buf2 * 1.E-2
2792           endif
2793     endif
2794     call mpp_land_bcast_int1(ierr)
2795     if(ierr == 0) call decompose_data_real (buf2,fpar)
2796     if(my_id .eq. io_id ) call get_2d_netcdf("LAI",     ncid, buf2,   units, global_nx, global_ny, .FALSE., ierr)
2797     call mpp_land_bcast_int1(ierr)
2798     if(ierr == 0) call decompose_data_real (buf2,lai)
2799
2800     deallocate(buf2)
2801 #else
2802     ierr = nf_open(trim(flnm), NF_NOWRITE, ncid)
2803     if (ierr /= 0) then
2804        write(*,'("READFORC Problem opening netcdf file: ''", A, "''")') trim(flnm)
2805        call hydro_stop("READFORC_HRLDAS")
2806     endif
2807     call get_2d_netcdf("T2D",     ncid, t,     units, ix, jx, .TRUE., ierr)
2808     call get_2d_netcdf("Q2D",     ncid, q,     units, ix, jx, .TRUE., ierr)
2809     call get_2d_netcdf("U2D",     ncid, u,     units, ix, jx, .TRUE., ierr)
2810     call get_2d_netcdf("V2D",     ncid, v,     units, ix, jx, .TRUE., ierr)
2811     call get_2d_netcdf("PSFC",    ncid, p,     units, ix, jx, .TRUE., ierr)
2812     call get_2d_netcdf("LWDOWN",  ncid, lw,    units, ix, jx, .TRUE., ierr)
2813     call get_2d_netcdf("SWDOWN",  ncid, sw,    units, ix, jx, .TRUE., ierr)
2814     call get_2d_netcdf("RAINRATE",ncid, pcp,   units, ix, jx, .TRUE., ierr)
2815     call get_2d_netcdf("VEGFRA",  ncid, fpar,  units, ix, jx, .FALSE., ierr)
2816
2817     if (ierr == 0) then
2818       if(maxval(fpar) .gt. 10 .and. maxval(fpar) .lt. 10000)  fpar = fpar * 1.E-2
2819     endif
2820     call get_2d_netcdf("LAI",     ncid, lai,   units, ix, jx, .FALSE., ierr)
2821 #endif
2822
2823     ierr = nf_close(ncid)
2824
2825   end subroutine READFORC_HRLDAS_mpp
2826
2827   subroutine READFORC_WRF_mpp(flnm,ix,jx,target_date,t,q,u,v,p,lw,sw,pcp,lai,fpar)
2828
2829     implicit none
2830     character(len=*),                   intent(in)  :: flnm
2831     integer,                            intent(in)  :: ix
2832     integer,                            intent(in)  :: jx
2833     character(len=*),                   intent(in)  :: target_date
2834     real,             dimension(ix,jx) :: t,q,u,v,p,lw,sw,pcp,pcpc, lai,fpar
2835     integer   tlevel
2836
2837     character(len=256) :: units
2838     integer :: ierr
2839     integer :: ncid
2840 #ifdef MPP_LAND
2841     real, allocatable, dimension(:,:) :: buf2
2842 #endif
2843
2844     tlevel = 1
2845
2846     pcpc = 0
2847
2848 #ifdef MPP_LAND
2849     if(my_id .eq. io_id) then
2850           allocate(buf2(global_nx, global_ny) )
2851     else
2852           allocate(buf2(1, 1) )
2853     endif
2854
2855     ! Open the NetCDF file.
2856
2857     if(my_id .eq. io_id) ierr = nf_open(flnm, NF_NOWRITE, ncid)
2858     call mpp_land_bcast_int1(ierr)
2859     if (ierr /= 0) then
2860        write(*,'("READFORC_WRF Problem opening netcdf file: ''", A, "''")') trim(flnm)
2861        call hydro_stop("In READFORC_WRF_mpp() - Problem opening netcdf file")
2862     endif
2863     if(my_id .eq. io_id) call get_2d_netcdf_ruc("T2",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2864     call decompose_data_real (buf2,t)
2865     if(my_id .eq. io_id) call get_2d_netcdf_ruc("Q2",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2866     call decompose_data_real (buf2,q)
2867     if(my_id .eq. io_id) call get_2d_netcdf_ruc("U10",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2868     call decompose_data_real (buf2,u)
2869     if(my_id .eq. io_id) call get_2d_netcdf_ruc("V10",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2870     call decompose_data_real (buf2,v)
2871     if(my_id .eq. io_id) call get_2d_netcdf_ruc("PSFC",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2872     call decompose_data_real (buf2,p)
2873     if(my_id .eq. io_id) call get_2d_netcdf_ruc("GLW",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2874     call decompose_data_real (buf2,lw)
2875     if(my_id .eq. io_id) call get_2d_netcdf_ruc("SWDOWN",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2876     call decompose_data_real (buf2,sw)
2877     if(my_id .eq. io_id) call get_2d_netcdf_ruc("RAINC",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2878     call decompose_data_real (buf2,pcpc)
2879     if(my_id .eq. io_id) call get_2d_netcdf_ruc("RAINNC",     ncid, buf2, global_nx, global_ny,tlevel, .true., ierr)
2880     call decompose_data_real (buf2,pcp)
2881     if(my_id .eq. io_id) call get_2d_netcdf_ruc("LAI",     ncid, buf2, global_nx, global_ny,tlevel, .false., ierr)
2882     call mpp_land_bcast_int1(ierr)
2883     if(ierr == 0) call decompose_data_real (buf2,lai)
2884     if(my_id .eq. io_id) then
2885        call get_2d_netcdf_ruc("VEGFRA", ncid, buf2, global_nx, global_ny, tlevel, .true., ierr)
2886        if(ierr == 0) then
2887           if(maxval(buf2) .gt. 10 .and. maxval(buf2) .lt. 10000)  buf2 = buf2 * 1.E-2
2888        endif
2889     endif
2890     call mpp_land_bcast_int1(ierr)
2891     if(ierr == 0) call decompose_data_real (buf2,fpar)
2892     deallocate(buf2)
2893 #else
2894
2895     ! Open the NetCDF file.
2896     ierr = nf_open(flnm, NF_NOWRITE, ncid)
2897     if (ierr /= 0) then
2898        write(*,'("READFORC_WRF Problem opening netcdf file: ''", A, "''")') trim(flnm)
2899        call hydro_stop("In READFORC_WRF_mpp() - Problem opening netcdf file")
2900     endif
2901     call get_2d_netcdf_ruc("T2",     ncid, t,     ix, jx,tlevel, .true., ierr)
2902     call get_2d_netcdf_ruc("Q2",     ncid, q,     ix, jx,tlevel, .true., ierr)
2903     call get_2d_netcdf_ruc("U10",    ncid, u,     ix, jx,tlevel, .true., ierr)
2904     call get_2d_netcdf_ruc("V10",    ncid, v,     ix, jx,tlevel, .true., ierr)
2905     call get_2d_netcdf_ruc("PSFC",   ncid, p,     ix, jx,tlevel, .true., ierr)
2906     call get_2d_netcdf_ruc("GLW",    ncid, lw,    ix, jx,tlevel, .true., ierr)
2907     call get_2d_netcdf_ruc("SWDOWN", ncid, sw,    ix, jx,tlevel, .true., ierr)
2908     call get_2d_netcdf_ruc("RAINC",  ncid, pcpc,  ix, jx,tlevel, .true., ierr)
2909     call get_2d_netcdf_ruc("RAINNC", ncid, pcp,   ix, jx,tlevel, .true., ierr)
2910     call get_2d_netcdf_ruc("VEGFRA", ncid, fpar,  ix, jx,tlevel, .false., ierr)
2911     if(ierr == 0) then
2912         if(maxval(fpar) .gt. 10 .and. (maxval(fpar) .lt. 10000) ) fpar = fpar * 0.01
2913     endif
2914     call get_2d_netcdf_ruc("LAI", ncid, lai,  ix, jx,tlevel, .false., ierr)
2915
2916 #endif
2917
2918
2919     pcp=pcp+pcpc   ! assumes pcpc=0 for resolved convection...
2920     ierr = nf_close(ncid)
2921
2922
2923   end subroutine READFORC_WRF_mpp
2924
2925   subroutine READFORC_MDV_mpp(flnm,ix,jx,pcp,mmflag,ierr_flg)
2926     implicit none
2927
2928     character(len=*),                   intent(in)  :: flnm
2929     integer,                            intent(in)  :: ix
2930     integer,                            intent(in)  :: jx
2931     integer,                            intent(out)  :: ierr_flg
2932     integer :: it,jew,zsn
2933     real,             dimension(ix,jx), intent(out) :: pcp
2934
2935     character(len=256) :: units
2936     integer :: ierr,i,j,i2,j2,varid
2937     integer :: ncid,mmflag
2938     real, dimension(ix,jx) :: temp
2939 #ifdef MPP_LAND
2940     real, allocatable, dimension(:,:) :: buf2
2941     if(my_id .eq. io_id) then
2942        allocate(buf2(global_nx, global_ny))
2943     else
2944        allocate(buf2(1,1))
2945     endif
2946 #endif
2947
2948     mmflag = 0   ! flag for units spec. (0=mm, 1=mm/s)
2949
2950
2951 !open NetCDF file...
2952 #ifdef MPP_LAND
2953       if(my_id .eq. io_id) then
2954 #endif
2955         ierr_flg = nf_open(flnm, NF_NOWRITE, ncid)
2956 #ifdef MPP_LAND
2957       endif
2958       call mpp_land_bcast_int1(ierr_flg)
2959 #endif
2960         if (ierr_flg /= 0) then
2961           write(*,'("READFORC_MDV Problem opening netcdf file: ''",A,"''")') &
2962                 trim(flnm)
2963 #ifdef MPP_LAND
2964            deallocate(buf2)
2965 #endif
2966            return
2967         end if
2968
2969 #ifdef MPP_LAND
2970       if(my_id .eq. io_id) then
2971 #endif
2972         ierr = nf_inq_varid(ncid,  "precip",  varid)
2973 #ifdef MPP_LAND
2974       endif
2975       call mpp_land_bcast_int1(ierr)
2976 #endif
2977         if(ierr /= 0) ierr_flg = ierr
2978         if (ierr /= 0) then
2979 #ifdef MPP_LAND
2980       if(my_id .eq. io_id) then
2981 #endif
2982           ierr = nf_inq_varid(ncid,  "precip_rate",  varid)   !recheck variable name...
2983 #ifdef MPP_LAND
2984       endif
2985       call mpp_land_bcast_int1(ierr)
2986 #endif
2987           if (ierr /= 0) then
2988 #ifdef MPP_LAND
2989           if(my_id .eq. io_id) then
2990 #endif
2991              ierr = nf_inq_varid(ncid,  "RAINRATE",  varid)   !recheck variable name...
2992 #ifdef MPP_LAND
2993           endif
2994           call mpp_land_bcast_int1(ierr)
2995 #endif
2996             if (ierr /= 0) then
2997               write(*,'("READFORC_MDV Problem reading precip netcdf file: ''", A,"''")') &
2998                  trim(flnm)
2999 #ifdef MPP_LAND
3000                 deallocate(buf2)
3001 #endif
3002                 return
3003             end if
3004           end if
3005           ierr_flg = ierr
3006           mmflag = 1
3007         end if
3008 #ifdef MPP_LAND
3009       if(my_id .eq. io_id) then
3010           ierr = nf_get_var_real(ncid, varid, buf2)
3011       endif
3012       call mpp_land_bcast_int1(ierr)
3013       if(ierr ==0) call decompose_data_real (buf2,pcp)
3014       deallocate(buf2)
3015 #else
3016         ierr = nf_get_var_real(ncid, varid, pcp)
3017 #endif
3018         if (ierr /= 0) then
3019            write(*,'("READFORC_MDV Problem reading netcdf file: ''", A,"''")') trim(flnm)
3020         end if
3021         ierr = nf_close(ncid)
3022
3023   end subroutine READFORC_MDV_mpp
3024
3025   subroutine READSNOW_FORC_mpp(flnm,ix,jx,weasd,snodep)
3026     implicit none
3027
3028     character(len=*),                   intent(in)  :: flnm
3029     integer,                            intent(in)  :: ix
3030     integer,                            intent(in)  :: jx
3031     real,             dimension(ix,jx), intent(out) :: weasd
3032     real,             dimension(ix,jx), intent(out) :: snodep
3033     real, dimension(ix,jx) :: tmp
3034
3035     character(len=256) :: units
3036     integer :: ierr
3037     integer :: ncid,i,j
3038 #ifdef MPP_LAND
3039     real, allocatable, dimension(:,:) :: buf2
3040     if(my_id .eq. io_id) then
3041        allocate(buf2(global_nx, global_ny))
3042     else
3043        allocate(buf2(1,1))
3044     endif
3045 #endif
3046
3047     ! Open the NetCDF file.
3048 #ifdef MPP_LAND
3049       if(my_id .eq. io_id) then
3050 #endif
3051     ierr = nf_open(flnm, NF_NOWRITE, ncid)
3052 #ifdef MPP_LAND
3053       endif
3054       call mpp_land_bcast_int1(ierr)
3055 #endif
3056     if (ierr /= 0) then
3057        write(*,'("READSNOW Problem opening netcdf file: ''", A, "''")') trim(flnm)
3058        call hydro_stop("In READSNOW_FORC_mpp() - Problem opening netcdf file")
3059     endif
3060
3061 #ifdef MPP_LAND
3062       if(my_id .eq. io_id) then
3063           call get_2d_netcdf("WEASD",  ncid, buf2,   units, ix, jx, .FALSE., ierr)
3064       endif
3065       call mpp_land_bcast_int1(ierr)
3066       if(ierr == 0) call decompose_data_real (buf2,tmp)
3067 #else
3068     call get_2d_netcdf("WEASD",  ncid, tmp,   units, ix, jx, .FALSE., ierr)
3069 #endif
3070     if (ierr /= 0) then
3071          call get_2d_netcdf("SNOW",  ncid, tmp,   units, ix, jx, .FALSE., ierr)
3072          if (ierr == 0) then
3073             units = "mm"
3074 #ifdef HYDRO_D
3075             print *, "read WEASD from wrfoutput ...... "
3076 #endif
3077             weasd = tmp * 1.E-3
3078          endif
3079     else
3080          weasd = tmp
3081          if (trim(units) == "m") then
3082             ! No conversion necessary
3083          else if (trim(units) == "mm") then
3084             ! convert WEASD from mm to m
3085             weasd = weasd * 1.E-3
3086          endif
3087     endif
3088
3089     if (ierr /= 0) then
3090        print *, "!!!!! NO WEASD present in input file...initialize to 0."
3091     endif
3092 #ifdef MPP_LAND
3093       if(my_id .eq. io_id) then
3094          call get_2d_netcdf("SNODEP",     ncid, buf2,   units, ix, jx, .FALSE., ierr)
3095       endif
3096       call mpp_land_bcast_int1(ierr)
3097       if(ierr == 0) call decompose_data_real (buf2,tmp)
3098 #else
3099     call get_2d_netcdf("SNODEP",     ncid, tmp,   units, ix, jx, .FALSE., ierr)
3100 #endif
3101     if (ierr /= 0) then
3102        ! Quick assumption regarding snow depth.
3103
3104 #ifdef MPP_LAND
3105       if(my_id .eq. io_id) then
3106          call get_2d_netcdf("SNOWH",     ncid, buf2,   units, ix, jx, .FALSE., ierr)
3107       endif
3108       call mpp_land_bcast_int1(ierr)
3109       if(ierr == 0) call decompose_data_real (buf2,tmp)
3110 #else
3111          call get_2d_netcdf("SNOWH",     ncid, tmp,   units, ix, jx, .FALSE., ierr)
3112 #endif
3113        if(ierr .eq. 0) then
3114 #ifdef HYDRO_D
3115             print *, "read snow depth from wrfoutput ... "
3116 #endif
3117             snodep = tmp
3118        endif
3119     else
3120        snodep = tmp
3121     endif
3122
3123     if (ierr /= 0) then
3124        ! Quick assumption regarding snow depth.
3125 !yw       snodep = weasd * 10.
3126        where(snodep .lt. weasd) snodep = weasd*10  !set lower bound to correct bi-lin interp err...
3127     endif
3128
3129 !DJG check for erroneous neg WEASD or SNOWD due to offline interpolation...
3130        where(snodep .lt. 0) snodep = 0
3131        where(weasd .lt. 0) weasd = 0
3132     ierr = nf_close(ncid)
3133
3134   end subroutine READSNOW_FORC_mpp
3135
3136   subroutine read_ldasout(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3137
3138       implicit none
3139       logical :: fexist
3140       integer :: ix,jx
3141       character(len=*) :: olddate,hgrid,indir
3142       character(len=19) :: outdate
3143       character(len=256) :: inflnm, inflnm2
3144       real :: dt
3145       real, dimension(ix,jx):: infxsrt,infxsrt2,soldrain,soldrain2
3146       integer :: ncid, ierr
3147       character(len=256) :: units
3148 #ifdef MPP_LAND
3149       real, dimension(global_nx,global_ny) :: gArr
3150 #endif
3151
3152         ! check for file with hours first
3153         inflnm = trim(indir)//"/"//&
3154              olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
3155              ".LDASOUT_DOMAIN"//hgrid
3156         inquire (file=trim(inflnm), exist=fexist)
3157
3158         if(.not. fexist) then
3159            ! check for file with minutes
3160              inflnm = trim(indir)//"/"//&
3161                   olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//olddate(15:16)//&
3162                   ".LDASOUT_DOMAIN"//hgrid
3163              inquire (file=trim(inflnm), exist=fexist)
3164         endif
3165         if(.not. fexist) then
3166             write(6,*) "Error: input file does not exist. Check ", trim(olddate)
3167             call hydro_stop( "LDASOUT input Error")
3168         endif
3169
3170         call geth_newdate(outdate,olddate,nint(dt))
3171         ! check file for next date
3172         ! check for file with hours first
3173         inflnm2 = trim(indir)//"/"//&
3174              outdate(1:4)//outdate(6:7)//outdate(9:10)//outdate(12:13)//&
3175              ".LDASOUT_DOMAIN"//hgrid
3176         inquire (file=trim(inflnm2), exist=fexist)
3177
3178         if(.not. fexist) then
3179            ! check for file with minutes
3180              inflnm2 = trim(indir)//"/"//&
3181                   outdate(1:4)//outdate(6:7)//outdate(9:10)//outdate(12:13)//outdate(15:16)//&
3182                   ".LDASOUT_DOMAIN"//hgrid
3183              inquire (file=trim(inflnm2), exist=fexist)
3184         endif
3185         if(.not. fexist) then
3186             write(6,*) "FATAL ERROR: input file does not exist. Check ", trim(outdate)
3187             call hydro_stop( "LDASOUT input Error")
3188         endif
3189 !       read file1
3190 #ifdef MPP_LAND
3191         if(my_id .eq. io_id) then
3192            ierr = nf_open(trim(inflnm), NF_NOWRITE, ncid)
3193            call get_2d_netcdf("SFCRNOFF",    ncid, gArr, units,  global_nx, global_ny, .TRUE., ierr)
3194         endif
3195         call decompose_data_real (gArr,infxsrt)
3196         if(my_id .eq. io_id) then
3197            call get_2d_netcdf("UGDRNOFF",    ncid, gArr, units, global_nx, global_ny, .TRUE., ierr)
3198         endif
3199         call decompose_data_real (gArr,soldrain)
3200         if(my_id .eq. io_id) then
3201             ierr = nf_close(ncid)
3202         endif
3203 #else
3204         ierr = nf_open(trim(inflnm), NF_NOWRITE, ncid)
3205         call get_2d_netcdf("SFCRNOFF",    ncid, infxsrt, units,  ix, jx, .TRUE., ierr)
3206         call get_2d_netcdf("UGDRNOFF",    ncid, soldrain, units,  ix, jx, .TRUE., ierr)
3207         ierr = nf_close(ncid)
3208 #endif
3209 !       read file2
3210 #ifdef MPP_LAND
3211        if(my_id .eq. io_id) then
3212            ierr = nf_open(trim(inflnm2), NF_NOWRITE, ncid)
3213            call get_2d_netcdf("SFCRNOFF",    ncid, gArr, units,  global_nx, global_ny, .TRUE., ierr)
3214         endif
3215         call decompose_data_real (gArr,infxsrt2)
3216         if(my_id .eq. io_id) then
3217            call get_2d_netcdf("UGDRNOFF",    ncid, gArr, units, global_nx, global_ny, .TRUE., ierr)
3218         endif
3219         call decompose_data_real (gArr,soldrain2)
3220         if(my_id .eq. io_id) then
3221            ierr = nf_close(ncid)
3222         endif
3223 #else
3224         ierr = nf_open(trim(inflnm2), NF_NOWRITE, ncid)
3225         call get_2d_netcdf("SFCRNOFF",    ncid, infxsrt2, units,  ix, jx, .TRUE., ierr)
3226         call get_2d_netcdf("UGDRNOFF",    ncid, soldrain2, units,  ix, jx, .TRUE., ierr)
3227         ierr = nf_close(ncid)
3228 #endif
3229
3230         infxsrt = infxsrt2 - infxsrt
3231         soldrain = soldrain2 - soldrain
3232
3233    end subroutine read_ldasout
3234
3235 !temporary for Noah model
3236
3237   subroutine read_ldasout_seq(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3238       implicit none
3239       logical :: fexist
3240       integer :: ix,jx
3241       character(len=*) :: olddate,hgrid,indir
3242       character(len=19) :: outdate
3243       character(len=256) :: inflnm, inflnm2
3244       real :: dt
3245       real, dimension(ix,jx):: infxsrt,infxsrt2,soldrain,soldrain2
3246       integer :: ncid, ierr
3247       character(len=256) :: units
3248
3249         ! check for file with hours first
3250         inflnm = trim(indir)//"/"//&
3251              olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//&
3252              ".LDASOUT_DOMAIN"//hgrid
3253         inquire (file=trim(inflnm), exist=fexist)
3254
3255         if(.not. fexist) then
3256            ! check for file with minutes
3257              inflnm = trim(indir)//"/"//&
3258                   olddate(1:4)//olddate(6:7)//olddate(9:10)//olddate(12:13)//olddate(15:16)//&
3259                   ".LDASOUT_DOMAIN"//hgrid
3260              inquire (file=trim(inflnm), exist=fexist)
3261         endif
3262         if(.not. fexist) then
3263             write(6,*) "FATAL ERROR: input file does not exist. Check ", trim(olddate)
3264             call hydro_stop( "LDASOUT input Error")
3265         endif
3266
3267         call geth_newdate(outdate,olddate,nint(dt))
3268         ! check file for next date
3269         ! check for file with hours first
3270         inflnm2 = trim(indir)//"/"//&
3271              outdate(1:4)//outdate(6:7)//outdate(9:10)//outdate(12:13)//&
3272              ".LDASOUT_DOMAIN"//hgrid
3273         inquire (file=trim(inflnm2), exist=fexist)
3274
3275         if(.not. fexist) then
3276            ! check for file with minutes
3277              inflnm2 = trim(indir)//"/"//&
3278                   outdate(1:4)//outdate(6:7)//outdate(9:10)//outdate(12:13)//outdate(15:16)//&
3279                   ".LDASOUT_DOMAIN"//hgrid
3280              inquire (file=trim(inflnm2), exist=fexist)
3281         endif
3282         if(.not. fexist) then
3283             write(6,*) "FATAL ERROR: input file does not exist. Check ", trim(outdate)
3284             call hydro_stop( "LDASOUT input Error")
3285         endif
3286 !       read file1
3287         ierr = nf_open(trim(inflnm), NF_NOWRITE, ncid)
3288         call get_2d_netcdf("SFCRNOFF",    ncid, infxsrt, units,  ix, jx, .TRUE., ierr)
3289         call get_2d_netcdf("UGDRNOFF",    ncid, soldrain, units,  ix, jx, .TRUE., ierr)
3290         ierr = nf_close(ncid)
3291 !       read file2
3292         ierr = nf_open(trim(inflnm2), NF_NOWRITE, ncid)
3293         call get_2d_netcdf("SFCRNOFF",    ncid, infxsrt2, units,  ix, jx, .TRUE., ierr)
3294         call get_2d_netcdf("UGDRNOFF",    ncid, soldrain2, units,  ix, jx, .TRUE., ierr)
3295         ierr = nf_close(ncid)
3296
3297         infxsrt = infxsrt2 - infxsrt
3298         soldrain = soldrain2 - soldrain
3299
3300    end subroutine read_ldasout_seq
3301 end module module_lsm_forcing
3302
3303      subroutine read_forc_ldasout(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3304       use module_lsm_forcing, only: read_ldasout
3305       implicit none
3306       integer :: ix,jx
3307       character(len=*) :: olddate,hgrid,indir
3308       real :: dt
3309       real, dimension(ix,jx):: infxsrt,soldrain
3310       call read_ldasout(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3311     end subroutine read_forc_ldasout
3312
3313     subroutine read_forc_ldasout_seq(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3314 ! temporary for Noah model
3315       use module_lsm_forcing, only: read_ldasout_seq
3316       implicit none
3317       integer :: ix,jx
3318       character(len=*) :: olddate,hgrid,indir
3319       real :: dt
3320       real, dimension(ix,jx):: infxsrt,soldrain
3321       call read_ldasout_seq(olddate,hgrid, indir, dt,ix,jx,infxsrt,soldrain)
3322     end subroutine read_forc_ldasout_seq