dyn_em/module_wps_io_arw.F

   1 MODULE module_wps_io_arw
   2
   3    USE module_optional_input
   4
   5    IMPLICIT NONE
   6 !! FROM MODULE_KINDS
   7
   8 !   The numerical data types defined in this module are:
   9 !      i_byte    - specification kind for byte (1-byte) integer variable
  10 !      i_short   - specification kind for short (2-byte) integer variable
  11 !      i_long    - specification kind for long (4-byte) integer variable
  12 !      i_llong   - specification kind for double long (8-byte) integer variable
  13 !      r_single  - specification kind for single precision (4-byte) real variable
  14 !      r_double  - specification kind for double precision (8-byte) real variable
  15 !      r_quad    - specification kind for quad precision (16-byte) real variable
  16 !
  17 !      i_kind    - generic specification kind for default integer
  18 !      r_kind    - generic specification kind for default floating point
  19 !
  20 !
  21 ! Integer type definitions below
  22
  23 ! Integer types
  24   integer, parameter, public  :: i_byte  = selected_int_kind(1)      ! byte  integer
  25   integer, parameter, public  :: i_short = selected_int_kind(4)      ! short integer
  26 ! integer, parameter, public  :: i_long  = selected_int_kind(8)      ! long  integer
  27   integer, parameter, public  :: i_long  = kind(1)                   ! long  integer
  28   integer, parameter, private :: llong_t = selected_int_kind(16)     ! llong integer
  29   integer, parameter, public  :: i_llong = max( llong_t, i_long )
  30
  31 ! Expected 8-bit byte sizes of the integer kinds
  32   integer, parameter, public :: num_bytes_for_i_byte  = 1
  33   integer, parameter, public :: num_bytes_for_i_short = 2
  34   integer, parameter, public :: num_bytes_for_i_long  = 4
  35   integer, parameter, public :: num_bytes_for_i_llong = 8
  36
  37 ! Define arrays for default definition
  38   integer, parameter, private :: num_i_kinds = 4
  39   integer, parameter, dimension( num_i_kinds ), private :: integer_types = (/ &
  40        i_byte, i_short, i_long,  i_llong  /)
  41   integer, parameter, dimension( num_i_kinds ), private :: integer_byte_sizes = (/ &
  42        num_bytes_for_i_byte, num_bytes_for_i_short, &
  43        num_bytes_for_i_long, num_bytes_for_i_llong  /)
  44
  45 ! Default values
  46 ! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT INTEGER TYPE KIND ***
  47   integer, parameter, private :: default_integer = 2  ! 1=byte,
  48                                                       ! 2=short,
  49                                                       ! 3=long,
  50                                                       ! 4=llong
  51   integer, parameter, public  :: i_kind = integer_types( default_integer )
  52   integer, parameter, public  :: num_bytes_for_i_kind = &
  53        integer_byte_sizes( default_integer )
  54
  55
  56 ! Real definitions below
  57
  58 ! Real types
  59   integer, parameter, public  :: r_single = selected_real_kind(6)  ! single precision
  60   integer, parameter, public  :: r_double = selected_real_kind(15) ! double precision
  61   integer, parameter, private :: quad_t   = selected_real_kind(20) ! quad precision
  62   integer, parameter, public  :: r_quad   = max( quad_t, r_double )
  63
  64 ! Expected 8-bit byte sizes of the real kinds
  65   integer, parameter, public :: num_bytes_for_r_single = 4
  66   integer, parameter, public :: num_bytes_for_r_double = 8
  67   integer, parameter, public :: num_bytes_for_r_quad   = 16
  68
  69 ! Define arrays for default definition
  70   integer, parameter, private :: num_r_kinds = 3
  71   integer, parameter, dimension( num_r_kinds ), private :: real_kinds = (/ &
  72        r_single, r_double, r_quad    /)
  73   integer, parameter, dimension( num_r_kinds ), private :: real_byte_sizes = (/ &
  74        num_bytes_for_r_single, num_bytes_for_r_double, &
  75        num_bytes_for_r_quad    /)
  76
  77 ! Default values
  78 ! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT REAL TYPE KIND ***
  79   integer, parameter, private :: default_real = 2  ! 1=single,
  80                                                    ! 2=double,
  81 !! END FROM MODULE_KINDS
  82
  83       !  Input 3D meteorological fields.
  84
  85
  86       REAL , DIMENSION(:,:,:) , ALLOCATABLE :: landuse_frac_input , &
  87                                                soil_top_cat_input , &
  88                                                soil_bot_cat_input
  89
  90
  91       !  Input 2D surface fields.
  92
  93       REAL , DIMENSION(:,:)   , ALLOCATABLE :: soilt010_input , soilt040_input , &
  94                                                soilt100_input , soilt200_input , &
  95                                                soilm010_input , soilm040_input , &
  96                                                soilm100_input , soilm200_input , &
  97                                                psfc_in,pmsl
  98
  99       REAL , DIMENSION(:,:)   , ALLOCATABLE :: lat_wind, lon_wind
 100
 101
 102       !  Local input arrays
 103
 104       REAL,DIMENSION(:,:),ALLOCATABLE :: dum2d
 105       INTEGER,DIMENSION(:,:),ALLOCATABLE :: idum2d
 106       REAL,DIMENSION(:,:,:),ALLOCATABLE :: dum3d
 107
 108       LOGICAL , SAVE :: first_time_in = .TRUE.
 109
 110       INTEGER :: flag_soilt010 , flag_soilt100 , flag_soilt200 , &
 111                  flag_soilm010 , flag_soilm100 , flag_soilm200
 112
 113 !   Some constants to allow simple dimensions in the defined types
 114 !   given below.
 115
 116       CHARACTER (LEN=256) , PRIVATE :: a_message
 117
 118
 119 CONTAINS
 120
 121    SUBROUTINE read_wps ( grid, filename, file_date_string, num_metgrid_levels )
 122
 123       USE module_soil_pre
 124       USE module_domain
 125
 126       IMPLICIT NONE
 127
 128       TYPE(domain) , INTENT(INOUT)  :: grid
 129       CHARACTER (LEN=19) , INTENT(IN) :: file_date_string
 130       CHARACTER (LEN=4)              :: dummychar
 131       CHARACTER (LEN=132)              :: VarName
 132       CHARACTER (LEN=150)             :: chartemp
 133       CHARACTER (*) , INTENT(IN) :: filename
 134
 135       INTEGER :: ids,ide,jds,jde,kds,kde           &
 136                 ,ims,ime,jms,jme,kms,kme           &
 137                 ,its,ite,jts,jte,kts,kte
 138
 139       INTEGER :: i , j , k , loop, IMAX, JMAX
 140       INTEGER :: DATAHANDLE, num_metgrid_levels
 141       INTEGER :: Sysdepinfo, Status
 142       INTEGER :: istatus,ioutcount,iret,index,ierr
 143
 144       integer :: nrecs,iunit, L,hor_size,hor_size_u,hor_size_v
 145
 146 !!
 147       character*132, allocatable :: datestr_all(:)
 148       character*132, allocatable :: varname_all(:)
 149       integer, allocatable       :: domainend_all(:,:)
 150       integer, allocatable       :: start_block(:)
 151       integer, allocatable       :: end_block(:)
 152       integer, allocatable       :: start_byte(:)
 153       integer, allocatable       :: end_byte(:)
 154       integer(kind=i_llong), allocatable           :: file_offset(:)
 155 !!
 156
 157       REAL :: dummy,tmp,garb
 158       REAL, ALLOCATABLE:: dumdata(:,:,:)
 159       REAL, ALLOCATABLE:: dumdata_u(:,:,:)
 160       REAL, ALLOCATABLE:: dumdata_v(:,:,:)
 161
 162       REAL :: lats16(16),lons16(16)
 163
 164       CHARACTER (LEN= 8) :: dummy_char
 165
 166       INTEGER :: ok , map_proj , ok_open, igarb,igarb2, N
 167       REAL :: pt
 168       INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
 169
 170 #if defined(DM_PARALLEL) && !defined(STUBMPI)
 171       INCLUDE "mpif.h"
 172
 173       SELECT CASE ( model_data_order )
 174          CASE ( DATA_ORDER_ZXY )
 175             kds = grid%sd31 ; kde = grid%ed31 ;
 176             ids = grid%sd32 ; ide = grid%ed32 ;
 177             jds = grid%sd33 ; jde = grid%ed33 ;
 178
 179             kms = grid%sm31 ; kme = grid%em31 ;
 180             ims = grid%sm32 ; ime = grid%em32 ;
 181             jms = grid%sm33 ; jme = grid%em33 ;
 182
 183             kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
 184             its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
 185             jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
 186
 187          CASE ( DATA_ORDER_XYZ )
 188             ids = grid%sd31 ; ide = grid%ed31 ;
 189             jds = grid%sd32 ; jde = grid%ed32 ;
 190             kds = grid%sd33 ; kde = grid%ed33 ;
 191
 192             ims = grid%sm31 ; ime = grid%em31 ;
 193             jms = grid%sm32 ; jme = grid%em32 ;
 194             kms = grid%sm33 ; kme = grid%em33 ;
 195
 196             its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
 197             jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
 198             kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
 199
 200          CASE ( DATA_ORDER_XZY )
 201             ids = grid%sd31 ; ide = grid%ed31 ;
 202             kds = grid%sd32 ; kde = grid%ed32 ;
 203             jds = grid%sd33 ; jde = grid%ed33 ;
 204
 205             ims = grid%sm31 ; ime = grid%em31 ;
 206             kms = grid%sm32 ; kme = grid%em32 ;
 207             jms = grid%sm33 ; jme = grid%em33 ;
 208
 209             its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
 210             kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
 211             jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
 212
 213       END SELECT
 214
 215       !  Initialize what soil temperature and moisture is available.
 216
 217
 218       flag_st000010 = 0
 219       flag_st010040 = 0
 220       flag_st040100 = 0
 221       flag_st100200 = 0
 222       flag_sm000010 = 0
 223       flag_sm010040 = 0
 224       flag_sm040100 = 0
 225       flag_sm100200 = 0
 226       flag_st010200 = 0
 227       flag_sm010200 = 0
 228
 229       flag_soilt010 = 0
 230       flag_soilt040 = 0
 231       flag_soilt100 = 0
 232       flag_soilt200 = 0
 233       flag_soilm010 = 0
 234       flag_soilm040 = 0
 235       flag_soilm100 = 0
 236       flag_soilm200 = 0
 237
 238       flag_sst      = 0
 239       flag_toposoil = 0
 240
 241       !  How many soil levels have we found?  Well, right now, none.
 242
 243       num_st_levels_input = 0
 244       num_sm_levels_input = 0
 245       st_levels_input = -1
 246       sm_levels_input = -1
 247
 248       !  Get the space for the data if this is the first time here.
 249
 250
 251         IF (.NOT. ALLOCATED (pmsl)              ) ALLOCATE ( pmsl(its:ite,jts:jte) )
 252         IF (.NOT. ALLOCATED (psfc_in)           ) ALLOCATE ( psfc_in(its:ite,jts:jte) )
 253
 254 !!! MPI IO
 255
 256       iunit=33
 257       call count_recs_wrf_binary_file(iunit, trim(fileName), nrecs)
 258
 259       allocate (datestr_all(nrecs))
 260       allocate (varname_all(nrecs))
 261       allocate (domainend_all(3,nrecs))
 262       allocate (start_block(nrecs))
 263       allocate (end_block(nrecs))
 264       allocate (start_byte(nrecs))
 265       allocate (end_byte(nrecs))
 266       allocate (file_offset(nrecs))
 267
 268       call inventory_wrf_binary_file(iunit, trim(filename), nrecs,  &
 269                       datestr_all,varname_all,domainend_all,   &
 270                       start_block,end_block,start_byte,end_byte,file_offset)
 271
 272       call mpi_file_open(mpi_comm_world, trim(filename),     &
 273                          mpi_mode_rdonly,mpi_info_null, iunit, ierr)
 274       if (ierr /= 0) then
 275        call wrf_error_fatal("Error opening file with mpi io")
 276       end if
 277
 278       VarName='CEN_LAT'
 279       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 280       if (iret /= 0) then
 281         print*,VarName," not found in file"
 282       else
 283
 284         call mpi_file_read_at(iunit,file_offset(index)+5*4,      &
 285                               garb,1,mpi_real4,                  &
 286                               mpi_status_ignore, ierr)
 287
 288         if (ierr /= 0) then
 289           print*,"Error reading ", VarName," using MPIIO"
 290         else
 291           print*,VarName, ' from MPIIO READ= ',garb
 292           CALL nl_set_cen_lat ( grid%id , garb )
 293         end if
 294       end if
 295
 296       VarName='CEN_LON'
 297       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 298         call mpi_file_read_at(iunit,file_offset(index)+5*4,      &
 299                               garb,1,mpi_real4,                  &
 300                               mpi_status_ignore, ierr)
 301           CALL nl_set_cen_lon ( grid%id , garb )
 302           CALL nl_set_stand_lon ( grid%id , garb )
 303
 304       VarName='POLE_LAT'
 305       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 306         call mpi_file_read_at(iunit,file_offset(index)+5*4,      &
 307                               garb,1,mpi_real4,                  &
 308                               mpi_status_ignore, ierr)
 309           CALL nl_set_pole_lat ( grid%id , garb )
 310
 311       VarName='TRUELAT1'
 312       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 313         call mpi_file_read_at(iunit,file_offset(index)+5*4,      &
 314                               garb,1,mpi_real4,                  &
 315                               mpi_status_ignore, ierr)
 316           CALL nl_set_truelat1 ( grid%id , garb )
 317
 318       VarName='TRUELAT2'
 319       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 320         call mpi_file_read_at(iunit,file_offset(index)+5*4,      &
 321                               garb,1,mpi_real4,                  &
 322                               mpi_status_ignore, ierr)
 323           CALL nl_set_truelat2 ( grid%id , garb )
 324
 325       VarName='MAP_PROJ'
 326       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 327
 328         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 329                               igarb,1,mpi_integer4,             &
 330                               mpi_status_ignore, ierr)
 331
 332           CALL  nl_set_map_proj( grid%id, igarb)
 333
 334       VarName='ISURBAN'
 335       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 336         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 337                               igarb,1,mpi_integer4,             &
 338                               mpi_status_ignore, ierr)
 339           CALL  nl_set_isurban ( grid%id, igarb)
 340
 341       VarName='ISWATER'
 342       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 343         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 344                               igarb,1,mpi_integer4,             &
 345                               mpi_status_ignore, ierr)
 346         write(a_message,*) 'setting iswater to be: ', igarb
 347         CALL wrf_message ( a_message )
 348         CALL nl_set_iswater (grid%id, igarb )
 349
 350       VarName='ISICE'
 351       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 352         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 353                               igarb2,1,mpi_integer4,             &
 354                               mpi_status_ignore, ierr)
 355         write(a_message,*) 'setting isice to be: ', igarb2
 356         CALL wrf_message ( a_message )
 357         CALL nl_set_isice (grid%id, igarb2 )
 358
 359       IF ( igarb .eq. 16 .and. igarb2 .eq. 24 ) THEN
 360       CALL nl_set_mminlu ( grid%id, 'USGS')
 361       ENDIF
 362
 363       IF ( igarb .eq. 17 .and. igarb2 .eq. 15 ) THEN
 364 ! ambiguous
 365       CALL nl_set_mminlu ( grid%id, 'MODIFIED_IGBP_MODIS_NOAH')
 366 !      CALL nl_set_mminlu ( grid%id, 'MODIS')
 367       ENDIF
 368
 369       IF ( igarb .eq. 15 .and. igarb2 .eq. 16 ) THEN
 370       CALL nl_set_mminlu ( grid%id, 'SiB')
 371       ENDIF
 372
 373       VarName='FLAG_SNOWH'
 374       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 375         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 376                               igarb,1,mpi_integer4,             &
 377                               mpi_status_ignore, ierr)
 378         flag_snowh=igarb
 379
 380       VarName='FLAG_SNOW'
 381       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 382         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 383                               igarb,1,mpi_integer4,             &
 384                               mpi_status_ignore, ierr)
 385         flag_snow=igarb
 386
 387       VarName='FLAG_METGRID'
 388       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 389         if (iret .eq. 0) then
 390         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 391                               igarb,1,mpi_integer4,             &
 392                               mpi_status_ignore, ierr)
 393         flag_metgrid=igarb
 394         endif
 395
 396       VarName='FLAG_SOILHGT'
 397       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 398         if (iret .eq. 0) then
 399         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 400                               igarb,1,mpi_integer4,             &
 401                               mpi_status_ignore, ierr)
 402         flag_soilhgt=igarb
 403         endif
 404
 405       VarName='FLAG_PSFC'
 406       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 407         if (iret .eq. 0) then
 408         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 409                               igarb,1,mpi_integer4,             &
 410                               mpi_status_ignore, ierr)
 411         flag_psfc=igarb
 412         endif
 413
 414       VarName='FLAG_SLP'
 415       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 416         if (iret .eq. 0) then
 417         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 418                               igarb,1,mpi_integer4,             &
 419                               mpi_status_ignore, ierr)
 420         flag_slp=igarb
 421         endif
 422
 423       VarName='NUM_METGRID_SOIL_LEVELS'
 424       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 425         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 426                               igarb,1,mpi_integer4,             &
 427                               mpi_status_ignore, ierr)
 428
 429         CALL nl_set_num_metgrid_soil_levels(grid%id, igarb)
 430         num_sw_levels_input=igarb
 431
 432
 433       VarName='FLAG_SOIL_LEVELS'
 434       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 435         if (iret .eq. 0) then
 436         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 437                               igarb,1,mpi_integer4,             &
 438                               mpi_status_ignore, ierr)
 439         flag_soil_levels=igarb
 440         endif
 441
 442
 443       VarName='FLAG_SOIL_LAYERS'
 444       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 445         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 446                               igarb,1,mpi_integer4,             &
 447                               mpi_status_ignore, ierr)
 448         flag_soil_layers=igarb
 449
 450
 451
 452       VarName='ISLAKE'
 453       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 454         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 455                               igarb,1,mpi_integer4,             &
 456                               mpi_status_ignore, ierr)
 457           CALL  nl_set_islake ( grid%id, igarb)
 458
 459       VarName='ISOILWATER'
 460       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 461         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 462                               igarb,1,mpi_integer4,             &
 463                               mpi_status_ignore, ierr)
 464           CALL  nl_set_isoilwater ( grid%id, igarb)
 465
 466       VarName='MOAD_CEN_LAT'
 467       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 468         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 469                               garb,1,mpi_real4,             &
 470                               mpi_status_ignore, ierr)
 471           CALL  nl_set_moad_cen_lat ( grid%id,garb)
 472
 473       VarName='corner_lats'
 474       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 475         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 476                               lats16,16,mpi_real4,             &
 477                               mpi_status_ignore, ierr)
 478
 479       VarName='corner_lons'
 480       call retrieve_index(index,VarName,varname_all,nrecs,iret)
 481         call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 482                               lons16,16,mpi_real4,             &
 483                               mpi_status_ignore, ierr)
 484
 485 !      VarName='MMINLU'
 486 !      call retrieve_index(index,VarName,varname_all,nrecs,iret)
 487 !       if (iret .eq. 0) then
 488 !       call mpi_file_read_at(iunit,file_offset(index)+5*4,     &
 489 !                              dummychar,1,mpi_char,             &
 490 !                              mpi_status_ignore, ierr)
 491 !        endif
 492
 493     hor_size=(IDE-IDS)*(JDE-JDS)
 494     hor_size_u=(IDE+1-IDS)*(JDE-JDS)
 495     hor_size_v=(IDE-IDS)*(JDE+1-JDS)
 496
 497     varName='PRES'
 498     allocate(dumdata(IDS:IDE-1,JDS:JDE-1,num_metgrid_levels))
 499     allocate(dumdata_u(IDS:IDE,JDS:JDE-1,num_metgrid_levels))
 500     allocate(dumdata_v(IDS:IDE-1,JDS:JDE,num_metgrid_levels))
 501
 502      CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 503      CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 504                           dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 505                           mpi_status_ignore, ierr)
 506
 507
 508        DO K=1,num_metgrid_levels
 509         DO J=JTS,min(JTE,JDE-1)
 510          DO I=ITS,min(ITE,IDE-1)
 511            grid%p_gc(I,K,J)=dumdata(I,J,K)
 512          ENDDO
 513         ENDDO
 514        ENDDO
 515
 516
 517     varName='GHT'
 518     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 519     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 520                              dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 521                              mpi_status_ignore, ierr)
 522
 523
 524        DO K=1,num_metgrid_levels
 525         DO J=JTS,min(JTE,JDE-1)
 526          DO I=ITS,min(ITE,IDE-1)
 527            grid%ght_gc(I,K,J)=dumdata(I,J,K)
 528          ENDDO
 529         ENDDO
 530        ENDDO
 531
 532     varName='VEGCAT'
 533     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 534     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 535                              dumdata,hor_size,mpi_real4,             &
 536                              mpi_status_ignore, ierr)
 537
 538         DO J=JTS,min(JTE,JDE-1)
 539          DO I=ITS,min(ITE,IDE-1)
 540            grid%vegcat(I,J)=dumdata(I,J,1)
 541          ENDDO
 542         ENDDO
 543
 544 !    varName='SOIL_CAT'
 545 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 546 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 547 !                             dumdata,hor_size,mpi_real4,             &
 548 !                             mpi_status_ignore, ierr)
 549 !        DO J=JTS,min(JTE,JDE-1)
 550 !         DO I=ITS,min(ITE,IDE-1)
 551 !           grid%input_soil_cat(I,J)=dumdata(I,J,1)
 552 !         ENDDO
 553 !        ENDDO
 554
 555     varName='CANWAT'
 556     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 557     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 558                              dumdata,hor_size,mpi_real4,             &
 559                              mpi_status_ignore, ierr)
 560
 561         DO J=JTS,min(JTE,JDE-1)
 562          DO I=ITS,min(ITE,IDE-1)
 563            grid%canwat(I,J)=dumdata(I,J,1)
 564          ENDDO
 565         ENDDO
 566
 567     varName='SNOW'
 568     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 569     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 570                              dumdata,hor_size,mpi_real4,             &
 571                              mpi_status_ignore, ierr)
 572         DO J=JTS,min(JTE,JDE-1)
 573          DO I=ITS,min(ITE,IDE-1)
 574            grid%snow(I,J)=dumdata(I,J,1)
 575          ENDDO
 576         ENDDO
 577
 578     varName='SNOWH'
 579     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 580     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 581                              dumdata,hor_size,mpi_real4,             &
 582                              mpi_status_ignore, ierr)
 583         DO J=JTS,min(JTE,JDE-1)
 584          DO I=ITS,min(ITE,IDE-1)
 585            grid%snowh(I,J)=dumdata(I,J,1)
 586          ENDDO
 587         ENDDO
 588
 589     varName='SKINTEMP'
 590     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 591     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 592                              dumdata,hor_size,mpi_real4,             &
 593                              mpi_status_ignore, ierr)
 594         DO J=JTS,min(JTE,JDE-1)
 595          DO I=ITS,min(ITE,IDE-1)
 596            grid%tsk_gc(I,J)=dumdata(I,J,1)
 597          ENDDO
 598         ENDDO
 599
 600     varName='SOILHGT'
 601     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 602     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 603                              dumdata,hor_size,mpi_real4,             &
 604                              mpi_status_ignore, ierr)
 605         DO J=JTS,min(JTE,JDE-1)
 606          DO I=ITS,min(ITE,IDE-1)
 607            grid%toposoil(I,J)=dumdata(I,J,1)
 608          ENDDO
 609         ENDDO
 610
 611 !    varName='SEAICE'
 612 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 613 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 614 !                             dumdata,hor_size,mpi_real4,             &
 615 !                             mpi_status_ignore, ierr)
 616 !        DO J=JTS,min(JTE,JDE-1)
 617 !         DO I=ITS,min(ITE,IDE-1)
 618 !           grid%xice_gc(I,J)=dumdata(I,J,1)
 619 !         ENDDO
 620 !        ENDDO
 621
 622     varName='ST100200'
 623     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 624     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 625                              dumdata,hor_size,mpi_real4,             &
 626                              mpi_status_ignore, ierr)
 627         DO J=JTS,min(JTE,JDE-1)
 628          DO I=ITS,min(ITE,IDE-1)
 629            grid%st100200(I,J)=dumdata(I,J,1)
 630          ENDDO
 631         ENDDO
 632
 633     varName='ST040100'
 634     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 635     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 636                              dumdata,hor_size,mpi_real4,             &
 637                              mpi_status_ignore, ierr)
 638         DO J=JTS,min(JTE,JDE-1)
 639          DO I=ITS,min(ITE,IDE-1)
 640            grid%st040100(I,J)=dumdata(I,J,1)
 641          ENDDO
 642         ENDDO
 643
 644     varName='ST010040'
 645     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 646     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 647                              dumdata,hor_size,mpi_real4,             &
 648                              mpi_status_ignore, ierr)
 649
 650         DO J=JTS,min(JTE,JDE-1)
 651          DO I=ITS,min(ITE,IDE-1)
 652            grid%st010040(I,J)=dumdata(I,J,1)
 653          ENDDO
 654         ENDDO
 655
 656     varName='ST000010'
 657     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 658     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 659                              dumdata,hor_size,mpi_real4,             &
 660                              mpi_status_ignore, ierr)
 661
 662         DO J=JTS,min(JTE,JDE-1)
 663          DO I=ITS,min(ITE,IDE-1)
 664            grid%st000010(I,J)=dumdata(I,J,1)
 665          ENDDO
 666         ENDDO
 667
 668     varName='SM100200'
 669     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 670     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 671                              dumdata,hor_size,mpi_real4,             &
 672                              mpi_status_ignore, ierr)
 673
 674         DO J=JTS,min(JTE,JDE-1)
 675          DO I=ITS,min(ITE,IDE-1)
 676            grid%sm100200(I,J)=dumdata(I,J,1)
 677          ENDDO
 678         ENDDO
 679
 680     varName='SM040100'
 681     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 682     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 683                              dumdata,hor_size,mpi_real4,             &
 684                              mpi_status_ignore, ierr)
 685
 686         DO J=JTS,min(JTE,JDE-1)
 687          DO I=ITS,min(ITE,IDE-1)
 688            grid%sm040100(I,J)=dumdata(I,J,1)
 689          ENDDO
 690         ENDDO
 691
 692     varName='SM010040'
 693     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 694     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 695                              dumdata,hor_size,mpi_real4,             &
 696                              mpi_status_ignore, ierr)
 697
 698         DO J=JTS,min(JTE,JDE-1)
 699          DO I=ITS,min(ITE,IDE-1)
 700            grid%sm010040(I,J)=dumdata(I,J,1)
 701          ENDDO
 702         ENDDO
 703
 704     varName='SM000010'
 705     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 706     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 707                              dumdata,hor_size,mpi_real4,             &
 708                              mpi_status_ignore, ierr)
 709
 710         DO J=JTS,min(JTE,JDE-1)
 711          DO I=ITS,min(ITE,IDE-1)
 712            grid%sm000010(I,J)=dumdata(I,J,1)
 713          ENDDO
 714         ENDDO
 715
 716     varName='PSFC'
 717     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 718     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 719                              dumdata,hor_size,mpi_real4,             &
 720                              mpi_status_ignore, ierr)
 721
 722         if ( ierr .eq. 0 ) flag_psfc=1
 723
 724         DO J=JTS,min(JTE,JDE-1)
 725          DO I=ITS,min(ITE,IDE-1)
 726            grid%psfc(I,J)=dumdata(I,J,1)
 727          ENDDO
 728         ENDDO
 729
 730     varName='RH'
 731     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 732     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 733                              dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 734                              mpi_status_ignore, ierr)
 735
 736        DO K=1,num_metgrid_levels
 737         DO J=JTS,min(JTE,JDE-1)
 738          DO I=ITS,min(ITE,IDE-1)
 739            grid%rh_gc(I,K,J)=dumdata(I,J,K)
 740          ENDDO
 741         ENDDO
 742        ENDDO
 743
 744     varName='VV'
 745     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 746     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 747                              dumdata_v,hor_size_v*num_metgrid_levels,mpi_real4,             &
 748                              mpi_status_ignore, ierr)
 749
 750        DO K=1,num_metgrid_levels
 751         DO J=JTS,min(JTE,JDE)
 752          DO I=ITS,min(ITE,IDE-1)
 753            grid%v_gc(I,K,J)=dumdata_v(I,J,K)
 754         if (grid%v_gc(I,K,J) .ne. grid%v_gc(I,K,J) .or. abs(grid%v_gc(I,K,J)) .gt. 100.) then
 755         write(a_message,*) 'bad v_gc defined: ', I,K,J,grid%v_gc(I,K,J)
 756         CALL wrf_message ( a_message )
 757         call wrf_error_fatal(" bad v_gc")
 758         endif
 759          ENDDO
 760         ENDDO
 761        ENDDO
 762
 763     varName='UU'
 764     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 765     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 766                              dumdata_u,hor_size_u*num_metgrid_levels,mpi_real4,             &
 767                              mpi_status_ignore, ierr)
 768        DO K=1,num_metgrid_levels
 769         DO J=JTS,min(JTE,JDE-1)
 770          DO I=ITS,min(ITE,IDE)
 771            grid%u_gc(I,K,J)=dumdata_u(I,J,K)
 772         if (grid%u_gc(I,K,J) .ne. grid%u_gc(I,K,J) .or. abs(grid%u_gc(I,K,J)) .gt. 100.) then
 773         write(a_message,*) 'bad u_gc defined: ', I,K,J,grid%u_gc(I,K,J)
 774         CALL wrf_message ( a_message )
 775         call wrf_error_fatal(" bad u_gc")
 776         endif
 777          ENDDO
 778         ENDDO
 779        ENDDO
 780
 781     varName='TT'
 782     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 783     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 784                              dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 785                              mpi_status_ignore, ierr)
 786        DO K=1,num_metgrid_levels
 787         DO J=JTS,min(JTE,JDE-1)
 788          DO I=ITS,min(ITE,IDE-1)
 789            grid%t_gc(I,K,J)=dumdata(I,J,K)
 790         if (grid%t_gc(I,K,J) .ne. grid%t_gc(I,K,J) .or. abs(grid%t_gc(I,K,J)) .gt. 350.) then
 791         write(a_message,*) 'bad t_gc defined: ', I,K,J,grid%t_gc(I,K,J)
 792         CALL wrf_message ( a_message )
 793         call wrf_error_fatal(" bad t_gc")
 794         endif
 795          ENDDO
 796         ENDDO
 797        ENDDO
 798
 799 !    varName='RWMR'
 800 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 801 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 802 !                             dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 803 !                             mpi_status_ignore, ierr)
 804 !       DO K=1,num_metgrid_levels
 805 !        DO J=JTS,min(JTE,JDE-1)
 806 !         DO I=ITS,min(ITE,IDE-1)
 807 !           grid%nmm_rwmr_gc(I,J,K)=dumdata(I,J,K)
 808 !         ENDDO
 809 !        ENDDO
 810 !       ENDDO
 811
 812 !    varName='SNMR'
 813 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 814 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 815 !                             dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 816 !                             mpi_status_ignore, ierr)
 817 !       DO K=1,num_metgrid_levels
 818 !        DO J=JTS,min(JTE,JDE-1)
 819 !         DO I=ITS,min(ITE,IDE-1)
 820 !           grid%nmm_snmr_gc(I,J,K)=dumdata(I,J,K)
 821 !         ENDDO
 822 !        ENDDO
 823 !       ENDDO
 824
 825 !    varName='CLWMR'
 826 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 827 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 828 !                             dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 829 !                             mpi_status_ignore, ierr)
 830 !       DO K=1,num_metgrid_levels
 831 !        DO J=JTS,min(JTE,JDE-1)
 832 !         DO I=ITS,min(ITE,IDE-1)
 833 !           grid%nmm_clwmr_gc(I,J,K)=dumdata(I,J,K)
 834 !         ENDDO
 835 !        ENDDO
 836 !       ENDDO
 837
 838 !    varName='CICE'
 839 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 840 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 841 !                             dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 842 !                             mpi_status_ignore, ierr)
 843 !       DO K=1,num_metgrid_levels
 844 !        DO J=JTS,min(JTE,JDE-1)
 845 !         DO I=ITS,min(ITE,IDE-1)
 846 !           grid%nmm_cice_gc(I,J,K)=dumdata(I,J,K)
 847 !         ENDDO
 848 !        ENDDO
 849 !       ENDDO
 850
 851 !    varName='FRIMEF'
 852 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 853 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 854 !                             dumdata,hor_size*num_metgrid_levels,mpi_real4,             &
 855 !                             mpi_status_ignore, ierr)
 856 !       DO K=1,num_metgrid_levels
 857 !        DO J=JTS,min(JTE,JDE-1)
 858 !         DO I=ITS,min(ITE,IDE-1)
 859 !           grid%nmm_rimef_gc(I,J,K)=dumdata(I,J,K)
 860 !         ENDDO
 861 !        ENDDO
 862 !       ENDDO
 863
 864
 865      varName='PMSL'
 866     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 867     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 868                              dumdata,hor_size,mpi_real4,             &
 869                              mpi_status_ignore, ierr)
 870         DO J=JTS,min(JTE,JDE-1)
 871          DO I=ITS,min(ITE,IDE-1)
 872            grid%pslv_gc(I,J)=dumdata(I,J,1)
 873          ENDDO
 874         ENDDO
 875
 876
 877     varName='SNOALB'
 878     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 879     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 880                              dumdata,hor_size,mpi_real4,             &
 881                              mpi_status_ignore, ierr)
 882
 883         DO J=JTS,min(JTE,JDE-1)
 884          DO I=ITS,min(ITE,IDE-1)
 885            grid%snoalb(I,J)=dumdata(I,J,1)
 886          ENDDO
 887         ENDDO
 888
 889          num_veg_cat      = SIZE ( grid%landusef , DIM=2 )
 890          num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
 891          num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
 892
 893     varName='GREENFRAC'
 894     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 895     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 896                              dumdata,hor_size*12,mpi_real4,             &
 897                              mpi_status_ignore, ierr)
 898
 899        DO K=1,12
 900         DO J=JTS,min(JTE,JDE-1)
 901          DO I=ITS,min(ITE,IDE-1)
 902            grid%greenfrac(I,K,J)=dumdata(I,J,K)
 903          ENDDO
 904         ENDDO
 905        ENDDO
 906
 907     varName='ALBEDO12M'
 908     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 909     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 910                              dumdata,hor_size*12,mpi_real4,             &
 911                              mpi_status_ignore, ierr)
 912
 913        DO K=1,12
 914         DO J=JTS,min(JTE,JDE-1)
 915          DO I=ITS,min(ITE,IDE-1)
 916            grid%albedo12m(I,K,J)=dumdata(I,J,K)
 917          ENDDO
 918         ENDDO
 919        ENDDO
 920
 921     varName='SOILCBOT'
 922     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 923     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 924                              dumdata,hor_size*num_soil_bot_cat,mpi_real4,             &
 925                              mpi_status_ignore, ierr)
 926
 927        DO K=1,num_soil_bot_cat
 928         DO J=JTS,min(JTE,JDE-1)
 929          DO I=ITS,min(ITE,IDE-1)
 930            grid%soilcbot(I,K,J)=dumdata(I,J,K)
 931          ENDDO
 932         ENDDO
 933        ENDDO
 934
 935     varName='SOILCAT'
 936     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 937     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 938                              dumdata,hor_size,mpi_real4,             &
 939                              mpi_status_ignore, ierr)
 940
 941         DO J=JTS,min(JTE,JDE-1)
 942          DO I=ITS,min(ITE,IDE-1)
 943            grid%soilcat(I,J)=dumdata(I,J,1)
 944          ENDDO
 945         ENDDO
 946
 947         write(a_message,*) 'veg_cat and soil_cat sizes:::: ', num_veg_cat , num_soil_top_cat
 948         CALL wrf_message ( a_message )
 949
 950     varName='SOILCTOP'
 951     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 952     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 953                              dumdata,hor_size*num_soil_top_cat,mpi_real4,             &
 954                              mpi_status_ignore, ierr)
 955        DO K=1,num_soil_top_cat
 956         DO J=JTS,min(JTE,JDE-1)
 957          DO I=ITS,min(ITE,IDE-1)
 958            grid%soilctop(I,K,J)=dumdata(I,J,K)
 959          ENDDO
 960         ENDDO
 961        ENDDO
 962
 963     varName='SOILTEMP'
 964     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 965     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 966                              dumdata,hor_size,mpi_real4,             &
 967                              mpi_status_ignore, ierr)
 968
 969         DO J=JTS,min(JTE,JDE-1)
 970          DO I=ITS,min(ITE,IDE-1)
 971            grid%tmn_gc(I,J)=dumdata(I,J,1)
 972          ENDDO
 973         ENDDO
 974
 975 !    varName='HGT_V'
 976 !    CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 977 !    CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 978 !                             dumdata,hor_size,mpi_real4,             &
 979 !                             mpi_status_ignore, ierr)
 980 !
 981 !        DO J=JTS,min(JTE,JDE-1)
 982 !         DO I=ITS,min(ITE,IDE-1)
 983 !           grid%nmm_htv_gc(I,J)=dumdata(I,J,1)
 984 !         ENDDO
 985 !        ENDDO
 986
 987     varName='HGT_M'
 988     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
 989     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
 990                              dumdata,hor_size,mpi_real4,             &
 991                              mpi_status_ignore, ierr)
 992
 993         DO J=JTS,min(JTE,JDE-1)
 994          DO I=ITS,min(ITE,IDE-1)
 995            grid%ht_gc(I,J)=dumdata(I,J,1)
 996          ENDDO
 997         ENDDO
 998
 999     varName='LU_INDEX'
1000     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1001     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1002                              dumdata,hor_size,mpi_real4,             &
1003                              mpi_status_ignore, ierr)
1004
1005         DO J=JTS,min(JTE,JDE-1)
1006          DO I=ITS,min(ITE,IDE-1)
1007            grid%lu_index(I,J)=dumdata(I,J,1)
1008          ENDDO
1009         ENDDO
1010
1011     varName='LANDUSEF'
1012     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1013     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1014                              dumdata,hor_size*num_veg_cat,mpi_real4,             &
1015                              mpi_status_ignore, ierr)
1016
1017        DO K=1,num_veg_cat
1018         DO J=JTS,min(JTE,JDE-1)
1019          DO I=ITS,min(ITE,IDE-1)
1020            grid%landusef(I,K,J)=dumdata(I,J,K)
1021          ENDDO
1022         ENDDO
1023        ENDDO
1024
1025     varName='LANDMASK'
1026     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1027     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1028                              dumdata,hor_size,mpi_real4,             &
1029                              mpi_status_ignore, ierr)
1030
1031         DO J=JTS,min(JTE,JDE-1)
1032          DO I=ITS,min(ITE,IDE-1)
1033            grid%landmask(I,J)=dumdata(I,J,1)
1034          ENDDO
1035         ENDDO
1036
1037     varName='F'
1038     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1039     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1040                              dumdata,hor_size,mpi_real4,             &
1041                              mpi_status_ignore, ierr)
1042         DO J=JTS,min(JTE,JDE-1)
1043          DO I=ITS,min(ITE,IDE-1)
1044            grid%f(I,J)=dumdata(I,J,1)
1045          ENDDO
1046         ENDDO
1047
1048     varName='E'
1049     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1050     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1051                              dumdata,hor_size,mpi_real4,             &
1052                              mpi_status_ignore, ierr)
1053         DO J=JTS,min(JTE,JDE-1)
1054          DO I=ITS,min(ITE,IDE-1)
1055            grid%e(I,J)=dumdata(I,J,1)
1056          ENDDO
1057         ENDDO
1058
1059     varName='COSALPHA'
1060     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1061     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1062                              dumdata,hor_size,mpi_real4,             &
1063                              mpi_status_ignore, ierr)
1064         DO J=JTS,min(JTE,JDE-1)
1065          DO I=ITS,min(ITE,IDE-1)
1066            grid%cosa(I,J)=dumdata(I,J,1)
1067          ENDDO
1068         ENDDO
1069
1070     varName='SINALPHA'
1071     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1072     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1073                              dumdata,hor_size,mpi_real4,             &
1074                              mpi_status_ignore, ierr)
1075         DO J=JTS,min(JTE,JDE-1)
1076          DO I=ITS,min(ITE,IDE-1)
1077            grid%sina(I,J)=dumdata(I,J,1)
1078          ENDDO
1079         ENDDO
1080
1081     varName='MAPFAC_U'
1082     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1083     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1084                              dumdata_u,hor_size_u,mpi_real4,             &
1085                              mpi_status_ignore, ierr)
1086         DO J=JTS,min(JTE,JDE-1)
1087          DO I=ITS,min(ITE,IDE)
1088            grid%msfu(I,J)=dumdata_u(I,J,1)
1089
1090 !! bandaid for newly created WPS geogrid static files
1091         if (grid%msfu(I,J) .lt. 0.7) then
1092         write(a_message,*) 'weird msfu at I,J: ', I,J,grid%msfu(I,J)
1093         CALL wrf_message ( a_message )
1094
1095         if(J .eq. min(JTE,JDE-1)) then
1096         grid%msfu(I,J)=dumdata_u(I,J-1,1)
1097         write(a_message,*) 'changing msfu to: ',I,J,  grid%msfu(I,J)
1098         CALL wrf_message ( a_message )
1099         endif
1100
1101         endif
1102 !! end bandaid
1103
1104          ENDDO
1105         ENDDO
1106
1107     varName='MAPFAC_V'
1108     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1109     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1110                              dumdata_v,hor_size_v,mpi_real4,             &
1111                              mpi_status_ignore, ierr)
1112         DO J=JTS,min(JTE,JDE)
1113          DO I=ITS,min(ITE,IDE-1)
1114            grid%msfv(I,J)=dumdata_v(I,J,1)
1115
1116         if (grid%msfv(I,J) .lt. 0.7 ) then
1117         write(a_message,*) 'weird msfv at I,J: ', I,J,grid%msfv(I,J)
1118         CALL wrf_message ( a_message )
1119         grid%msfv(I,J)=dumdata_v(I,J-1,1)
1120         if( J .eq. min(JTE,JDE)) then
1121         write(a_message,*) 'changing msfv to: ',I,J,  grid%msfv(I,J)
1122         CALL wrf_message ( a_message )
1123         endif
1124
1125         endif
1126
1127          ENDDO
1128         ENDDO
1129
1130     varName='MAPFAC_M'
1131     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1132     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1133                              dumdata,hor_size,mpi_real4,             &
1134                              mpi_status_ignore, ierr)
1135         DO J=JTS,min(JTE,JDE-1)
1136          DO I=ITS,min(ITE,IDE-1)
1137            grid%msft(I,J)=dumdata(I,J,1)
1138          ENDDO
1139         ENDDO
1140
1141     varName='CLAT'
1142     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1143     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1144                              dumdata,hor_size,mpi_real4,             &
1145                              mpi_status_ignore, ierr)
1146         DO J=JTS,min(JTE,JDE-1)
1147          DO I=ITS,min(ITE,IDE-1)
1148            grid%clat(I,J)=dumdata(I,J,1)
1149          ENDDO
1150         ENDDO
1151
1152     varName='XLONG_U'
1153     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1154     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1155                              dumdata_u,hor_size_u,mpi_real4,             &
1156                              mpi_status_ignore, ierr)
1157         DO J=JTS,min(JTE,JDE-1)
1158          DO I=ITS,min(ITE,IDE)
1159            grid%xlong_u(I,J)=dumdata_u(I,J,1)
1160          ENDDO
1161         ENDDO
1162
1163     varName='XLAT_U'
1164     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1165     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1166                              dumdata_u,hor_size_u,mpi_real4,             &
1167                              mpi_status_ignore, ierr)
1168         DO J=JTS,min(JTE,JDE-1)
1169          DO I=ITS,min(ITE,IDE)
1170            grid%xlat_u(I,J)=dumdata_u(I,J,1)
1171          ENDDO
1172         ENDDO
1173
1174     varName='XLONG_V'
1175     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1176     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1177                              dumdata_v,hor_size_v,mpi_real4,             &
1178                              mpi_status_ignore, ierr)
1179         DO J=JTS,min(JTE,JDE)
1180          DO I=ITS,min(ITE,IDE-1)
1181            grid%xlong_v(I,J)=dumdata_v(I,J,1)
1182          ENDDO
1183         ENDDO
1184
1185     varName='XLAT_V'
1186     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1187     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1188                              dumdata_v,hor_size_v,mpi_real4,             &
1189                              mpi_status_ignore, ierr)
1190         DO J=JTS,min(JTE,JDE)
1191          DO I=ITS,min(ITE,IDE-1)
1192            grid%xlat_v(I,J)=dumdata_v(I,J,1)
1193          ENDDO
1194         ENDDO
1195
1196     varName='XLONG_M'
1197     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1198     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1199                              dumdata,hor_size,mpi_real4,             &
1200                              mpi_status_ignore, ierr)
1201
1202         DO J=JTS,min(JTE,JDE-1)
1203          DO I=ITS,min(ITE,IDE-1)
1204            grid%xlong_gc(I,J)=dumdata(I,J,1)
1205          ENDDO
1206         ENDDO
1207
1208         write(a_message,*) 'reading XLAT_M'
1209         CALL wrf_message ( a_message )
1210     varName='XLAT_M'
1211     CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1212     CALL mpi_file_read_at(iunit,file_offset(index+1),     &
1213                              dumdata,hor_size,mpi_real4,             &
1214                              mpi_status_ignore, ierr)
1215
1216         DO J=JTS,min(JTE,JDE-1)
1217          DO I=ITS,min(ITE,IDE-1)
1218            grid%xlat_gc(I,J)=dumdata(I,J,1)
1219          ENDDO
1220         ENDDO
1221         write(a_message,*) 'xlat_gc defined'
1222         CALL wrf_message ( a_message )
1223
1224       call mpi_file_close(mpi_comm_world, ierr)
1225
1226         write(a_message,*) 'to ST000010 def'
1227         CALL wrf_message ( a_message )
1228        varName='ST000010'
1229        flag_st000010 = 1
1230        num_st_levels_input = num_st_levels_input + 1
1231        st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
1232        DO J=JTS,min(JTE,JDE-1)
1233         DO I=ITS,min(ITE,IDE-1)
1234            st_input(I,num_st_levels_input + 1,J) = grid%st000010(i,j)
1235         ENDDO
1236        ENDDO
1237         write(a_message,*) 'past ST000010 def'
1238         CALL wrf_message ( a_message )
1239
1240        varName='ST010040'
1241        flag_st010040 = 1
1242        num_st_levels_input = num_st_levels_input + 1
1243        st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
1244        DO J=JTS,min(JTE,JDE-1)
1245         DO I=ITS,min(ITE,IDE-1)
1246             st_input(I,num_st_levels_input + 1,J) = grid%st010040(i,j)
1247         ENDDO
1248        ENDDO
1249
1250        varName='ST040100'
1251        flag_st040100 = 1
1252        num_st_levels_input = num_st_levels_input + 1
1253        st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
1254        DO J=JTS,min(JTE,JDE-1)
1255         DO I=ITS,min(ITE,IDE-1)
1256               st_input(I,num_st_levels_input + 1,J) = grid%st040100(i,j)
1257         ENDDO
1258        ENDDO
1259
1260        varName='ST100200'
1261        flag_st100200 = 1
1262        num_st_levels_input = num_st_levels_input + 1
1263        st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
1264        DO J=JTS,min(JTE,JDE-1)
1265         DO I=ITS,min(ITE,IDE-1)
1266               st_input(I,num_st_levels_input + 1,J) = grid%st100200(i,j)
1267         ENDDO
1268        ENDDO
1269
1270        varName='SM000010'
1271        flag_sm000010 = 1
1272        num_sm_levels_input = num_sm_levels_input + 1
1273        sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
1274        DO J=JTS,min(JTE,JDE-1)
1275         DO I=ITS,min(ITE,IDE-1)
1276            sm_input(I,num_sm_levels_input + 1,J) = grid%sm000010(i,j)
1277         ENDDO
1278        ENDDO
1279
1280        varName='SM010040'
1281        flag_sm010040 = 1
1282        num_sm_levels_input = num_sm_levels_input + 1
1283        sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
1284        DO J=JTS,min(JTE,JDE-1)
1285         DO I=ITS,min(ITE,IDE-1)
1286             sm_input(I,num_sm_levels_input + 1,J) = grid%sm010040(i,j)
1287         ENDDO
1288        ENDDO
1289
1290        varName='SM040100'
1291        flag_sm040100 = 1
1292        num_sm_levels_input = num_sm_levels_input + 1
1293        sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
1294        DO J=JTS,min(JTE,JDE-1)
1295         DO I=ITS,min(ITE,IDE-1)
1296               sm_input(I,num_sm_levels_input + 1,J) = grid%sm040100(i,j)
1297         ENDDO
1298        ENDDO
1299
1300        varName='SM100200'
1301        flag_sm100200 = 1
1302        num_sm_levels_input = num_sm_levels_input + 1
1303        sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
1304        DO J=JTS,min(JTE,JDE-1)
1305         DO I=ITS,min(ITE,IDE-1)
1306               sm_input(I,num_sm_levels_input + 1,J) = grid%sm100200(i,j)
1307         ENDDO
1308        ENDDO
1309
1310 !           flag_sst = 1
1311
1312         write(a_message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
1313         CALL wrf_message ( a_message )
1314         write(a_message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
1315         CALL wrf_message ( a_message )
1316         write(a_message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
1317         CALL wrf_message ( a_message )
1318         write(a_message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
1319         CALL wrf_message ( a_message )
1320
1321         DEALLOCATE(pmsl)
1322         DEALLOCATE(psfc_in)
1323         DEALLOCATE(dumdata)
1324         DEALLOCATE(dumdata_u)
1325         DEALLOCATE(dumdata_v)
1326
1327 #endif
1328      end subroutine read_wps
1329
1330 ! -------------------------------------------------------------------------
1331 ! -------------------------------------------------------------------------
1332
1333 !!!! MPI-IO pieces
1334
1335 subroutine retrieve_index(index,string,varname_all,nrecs,iret)
1336 !$$$  subprogram documentation block
1337 !                .      .    .                                       .
1338 ! subprogram:    retrieve_index  get record number of desired variable
1339 !   prgmmr: parrish          org: np22                date: 2004-11-29
1340 !
1341 ! abstract: by examining previously generated inventory of wrf binary restart file,
1342 !             find record number that contains the header record for variable
1343 !             identified by input character variable "string".
1344 !
1345 ! program history log:
1346 !   2004-11-29  parrish
1347 !
1348 !   input argument list:
1349 !     string           - mnemonic for variable desired
1350 !     varname_all      - list of all mnemonics obtained from inventory of file
1351 !     nrecs            - total number of sequential records counted in wrf
1352 !                        binary restart file
1353 !
1354 !   output argument list:
1355 !     index            - desired record number
1356 !     iret             - return status, set to 0 if variable was found,
1357 !                        non-zero if not.
1358 !
1359 ! attributes:
1360 !   language: f90
1361 !   machine:  ibm RS/6000 SP
1362 !
1363 !$$$
1364
1365   implicit none
1366
1367   integer,intent(out)::iret
1368   integer,intent(in)::nrecs
1369   integer,intent(out):: index
1370   character(*), intent(in):: string
1371   character(132),intent(in)::varname_all(nrecs)
1372
1373   integer i
1374
1375   iret=0
1376
1377   do i=1,nrecs
1378    if(trim(string) == trim(varname_all(i))) then
1379       index=i
1380       return
1381    end if
1382   end do
1383
1384   write(6,*)' problem reading wrf nmm binary file, rec id "',trim(string),'" not found'
1385
1386   iret=-1
1387
1388 end subroutine retrieve_index
1389 subroutine next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1390 !$$$  subprogram documentation block
1391 !                .      .    .                                       .
1392 ! subprogram:    next_buf    bring in next direct access block
1393 !   prgmmr: parrish          org: np22                date: 2004-11-29
1394 !
1395 ! abstract: bring in next direct access block when needed, as the file is scanned
1396 !             from beginning to end during counting and inventory of records.
1397 !             (subroutines count_recs_wrf_binary_file and inventory_wrf_binary_file)
1398 !
1399 ! program history log:
1400 !   2004-11-29  parrish
1401 !
1402 !   input argument list:
1403 !     in_unit          - fortran unit number where input file is opened through.
1404 !     nextbyte         - byte number from beginning of file that is desired
1405 !     locbyte          - byte number from beginning of last block read for desired byt
1406 !     lrecl            - direct access block length
1407 !     nreads           - number of blocks read before now (for diagnostic information
1408 !     lastbuf          - logical, if true, then no more blocks, so return
1409 !
1410 !   output argument list:
1411 !     buf              - output array containing contents of next block
1412 !     locbyte          - byte number from beginning of new block read for desired byte
1413 !     thisblock        - number of new block being read by this routine
1414 !     nreads           - number of blocks read now (for diagnostic information only)
1415 !     lastbuf          - logical, if true, then at end of file.
1416 !
1417 ! attributes:
1418 !   language: f90
1419 !   machine:  ibm RS/6000 SP
1420 !
1421 !$$$
1422 !  use kinds, only: i_byte,i_llong
1423   implicit none
1424
1425   integer(i_llong) lrecl
1426   integer in_unit,nreads
1427   integer(i_byte) buf(lrecl)
1428   integer(i_llong) nextbyte,locbyte,thisblock
1429   logical lastbuf
1430
1431   integer ierr
1432
1433   if(lastbuf) return
1434
1435   ierr=0
1436   nreads=nreads+1
1437
1438 !  compute thisblock:
1439
1440   thisblock = 1_i_llong + (nextbyte-1_i_llong)/lrecl
1441
1442   locbyte = 1_i_llong+mod(locbyte-1_i_llong,lrecl)
1443
1444   read(in_unit,rec=thisblock,iostat=ierr)buf
1445   lastbuf = ierr /= 0
1446
1447 end subroutine next_buf
1448
1449 subroutine inventory_wrf_binary_file(in_unit,wrf_ges_filename,nrecs, &
1450                                      datestr_all,varname_all,domainend_all, &
1451                                      start_block,end_block,start_byte,end_byte,file_offset)
1452 !$$$  subprogram documentation block
1453 !                .      .    .                                       .
1454 ! subprogram:    inventory_wrf_binary_file  get contents of wrf binary file
1455 !   prgmmr: parrish          org: np22                date: 2004-11-29
1456 !
1457 ! abstract: generate list of contents and map of wrf binary file which can be
1458 !             used for reading and writing with mpi io routines.
1459 !             same basic routine as count_recs_wrf_binary_file, except
1460 !             now wrf unpacking routines are used to decode wrf header
1461 !             records, and send back lists of variable mnemonics, dates,
1462 !             grid dimensions, and byte addresses relative to start of
1463 !             file for each field (this is used by mpi io routines).
1464 !
1465 ! program history log:
1466 !   2004-11-29  parrish
1467 !
1468 !
1469 !   input argument list:
1470 !     in_unit          - fortran unit number where input file is opened through.
1471 !     wrf_ges_filename - filename of input wrf binary restart file
1472 !     nrecs            - number of sequential records found on input wrf binary restart file.
1473 !                          (obtained by a previous call to count_recs_wrf_binary_file)
1474 !
1475 !   output argument list:  (all following dimensioned nrecs)
1476 !     datestr_all      - date character string for each field, where applicable (or else blanks)
1477 !     varname_all      - wrf mnemonic for each variable, where applicable (or blank)
1478 !     domainend_all    - dimensions of each field, where applicable (up to 3 dimensions)
1479 !     start_block      - direct access block number containing 1st byte of record
1480 !                            (after 4 byte record mark)
1481 !     end_block        - direct access block number containing last byte of record
1482 !                            (before 4 byte record mark)
1483 !     start_byte       - relative byte address in direct access block of 1st byte of record
1484 !     end_byte         - relative byte address in direct access block of last byte of record
1485 !     file_offset      - absolute address of byte before 1st byte of record (used by mpi io)
1486 !
1487 ! attributes:
1488 !   language: f90
1489 !   machine:  ibm RS/6000 SP
1490 !
1491 !$$$
1492 !   use kinds, only: r_single,i_byte,i_long,i_llong
1493   use module_internal_header_util
1494   implicit none
1495
1496   integer,intent(in)::in_unit,nrecs
1497   character(*),intent(in)::wrf_ges_filename
1498   character(132),intent(out)::datestr_all(nrecs),varname_all(nrecs)
1499   integer,intent(out)::domainend_all(3,nrecs)
1500   integer,intent(out)::start_block(nrecs),end_block(nrecs)
1501   integer,intent(out)::start_byte(nrecs),end_byte(nrecs)
1502   integer(i_llong),intent(out)::file_offset(nrecs)
1503
1504   integer irecs
1505   integer(i_llong) nextbyte,locbyte,thisblock
1506   integer(i_byte) lenrec4(4)
1507   integer(i_long) lenrec,lensave
1508   equivalence (lenrec4(1),lenrec)
1509   integer(i_byte) missing4(4)
1510   integer(i_long) missing
1511   equivalence (missing,missing4(1))
1512   integer(i_llong),parameter:: lrecl=2**20
1513   integer(i_byte) buf(lrecl)
1514   integer i,loc_count,nreads
1515   logical lastbuf
1516   integer(i_byte) hdrbuf4(2048)
1517   integer(i_long) hdrbuf(512)
1518   equivalence(hdrbuf(1),hdrbuf4(1))
1519   integer,parameter:: int_field       =       530
1520   integer,parameter:: int_dom_ti_char =       220
1521   integer,parameter:: int_dom_ti_real =       140
1522   integer,parameter:: int_dom_ti_integer =       180
1523   integer hdrbufsize
1524   integer inttypesize
1525   integer datahandle,count
1526   character(128) element,dumstr,strdata
1527   integer loccode
1528   character(132) blanks
1529   integer typesize
1530   integer fieldtype,comm,iocomm
1531   integer domaindesc
1532   character(132) memoryorder,stagger,dimnames(3)
1533   integer domainstart(3),domainend(3)
1534   integer memorystart(3),memoryend(3)
1535   integer patchstart(3),patchend(3)
1536   character(132) datestr,varname
1537 #if (WRFPLUS == 1)
1538   real(r_double) dummy_field(1)
1539 #else
1540   real dummy_field(1)
1541 #endif
1542 !  real(r_single) dummy_field(1)
1543 !  integer dummy_field
1544 !  real dummy_field
1545   integer itypesize
1546   integer idata(1)
1547   real rdata(1)
1548
1549   call wrf_sizeof_integer(itypesize)
1550   inttypesize=itypesize
1551
1552   blanks=trim(' ')
1553
1554         write(a_message,*) 'inventory subroutine'
1555         CALL wrf_message ( a_message )
1556
1557         write(a_message,*) 'opening file : ', trim(wrf_ges_filename)
1558         CALL wrf_message ( a_message )
1559
1560   open(in_unit,file=trim(wrf_ges_filename),access='direct',recl=lrecl)
1561   irecs=0
1562   missing=-9999
1563   nextbyte=0_i_llong
1564   locbyte=lrecl
1565   nreads=0
1566   lastbuf=.false.
1567   do
1568
1569 !   get length of next record
1570
1571     do i=1,4
1572      nextbyte=nextbyte+1_i_llong
1573      locbyte=locbyte+1_i_llong
1574      if(locbyte > lrecl .and. lastbuf) go to 900
1575      if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1576      lenrec4(i)=buf(locbyte)
1577     end do
1578     if(lenrec <= 0 .and. lastbuf) go to 900
1579     if(lenrec <= 0 .and. .not. lastbuf) go to 885
1580     nextbyte=nextbyte+1_i_llong
1581     locbyte=locbyte+1_i_llong
1582     if(locbyte > lrecl .and. lastbuf) go to 900
1583     if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1584
1585     irecs=irecs+1
1586     start_block(irecs)=thisblock
1587     start_byte(irecs)=locbyte
1588     file_offset(irecs)=nextbyte-1_i_llong
1589     hdrbuf4(1)=buf(locbyte)
1590     hdrbuf4(2:4)=missing4(2:4)
1591     hdrbuf4(5:8)=missing4(1:4)
1592     datestr_all(irecs)=blanks
1593     varname_all(irecs)=blanks
1594     domainend_all(1:3,irecs)=0
1595
1596     loc_count=1
1597 !       write(a_message,*) '1, hdrbuf4(1): 1', hdrbuf4(1)
1598 !       CALL wrf_message ( a_message )
1599     do i=2,8
1600        if(loc_count.ge.lenrec) exit
1601        loc_count=loc_count+1
1602        nextbyte=nextbyte+1_i_llong
1603        locbyte=locbyte+1_i_llong
1604        if(locbyte > lrecl .and. lastbuf) go to 900
1605        if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1606        hdrbuf4(i)=buf(locbyte)
1607 !       write(0,*) 'I, hdrbuf4(I): ', I, hdrbuf4(I)
1608     end do
1609
1610 !         if(lenrec==2048) then write(0,*)' irecs,hdrbuf(2),int_dom_ti_char,int_field=', &
1611 !                                      irecs,hdrbuf(2),int_dom_ti_char,int_field, int_dom_ti_real, int_dom_ti_integer
1612
1613     if(lenrec==2048.and.(hdrbuf(2) == int_dom_ti_char .or. hdrbuf(2) == int_field &
1614     .or. hdrbuf(2) == int_dom_ti_real .or. hdrbuf(2) == int_dom_ti_integer)) then
1615
1616 !    bring in next full record, so we can unpack datestr, varname, and domainend
1617        do i=9,lenrec
1618           loc_count=loc_count+1
1619           nextbyte=nextbyte+1_i_llong
1620           locbyte=locbyte+1_i_llong
1621           if(locbyte > lrecl .and. lastbuf) go to 900
1622           if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1623           hdrbuf4(i)=buf(locbyte)
1624 !       write(0,*) 'I, hdrbuf4(I): ', I, hdrbuf4(I)
1625        end do
1626 !       write(0,*) 'past 9,lenrec loop'
1627 !       write(0,*) 'hdrbuf(2): ', hdrbuf(2)
1628
1629        if(hdrbuf(2) == int_dom_ti_char) then
1630
1631           call int_get_ti_header_char(hdrbuf,hdrbufsize,inttypesize, &
1632                    datahandle,element,dumstr,strdata,loccode)
1633           varname_all(irecs)=trim(element)
1634           datestr_all(irecs)=trim(strdata)
1635               write(a_message,*)'(1) irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),trim(datestr_all(irecs))
1636               CALL wrf_message ( a_message )
1637
1638        else if(hdrbuf(2) == int_dom_ti_real) then
1639
1640 !       write(0,*) 'call int_get_ti_header_real'
1641
1642 !       write(0,*) 'hdrbuf(1:10): ', hdrbuf(1:10)
1643
1644 !       write(0,*) 'call with inttypesize typesize: ', inttypesize,typesize
1645
1646           call int_get_ti_header_real(hdrbuf,hdrbufsize,inttypesize,typesize, &
1647                    datahandle,element,rdata,count,loccode)
1648
1649 !       write(0,*) 'return int_get_ti_header_real'
1650
1651           varname_all(irecs)=trim(element)
1652
1653 !       write(0,*) 'varname_all(irecs): ', varname_all(irecs)
1654
1655 !          datestr_all(irecs)=trim(strdata)
1656
1657               write(a_message,*)'(2) irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),rdata(1)
1658               CALL wrf_message ( a_message )
1659 !       write(0,*) ' --------------------------- '
1660
1661        else if(hdrbuf(2) == int_dom_ti_integer) then
1662
1663           call int_get_ti_header_integer(hdrbuf,hdrbufsize,inttypesize,typesize, &
1664                    datahandle,element,idata,count,loccode)
1665           varname_all(irecs)=trim(element)
1666 !          datestr_all(irecs)=trim(strdata)
1667               write(0,*)'(3) irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),idata(1:count)
1668
1669        else
1670
1671 !       write(0,*) 'call int_get_write_field_header'
1672           call int_get_write_field_header(hdrbuf,hdrbufsize,inttypesize,typesize, &
1673              datahandle,datestr,varname,dummy_field,fieldtype,comm,iocomm, &
1674              domaindesc,memoryorder,stagger,dimnames, &
1675              domainstart,domainend,memorystart,memoryend,patchstart,patchend)
1676           varname_all(irecs)=trim(varname)
1677           datestr_all(irecs)=trim(datestr)
1678           domainend_all(1:3,irecs)=domainend(1:3)
1679 !              write(0,*)'(4) irecs,datestr,domend,varname = ', &
1680 !                  irecs,trim(datestr_all(irecs)),domainend_all(1:3,irecs),trim(varname_all(irecs))
1681
1682        end if
1683     end if
1684
1685     nextbyte=nextbyte-loc_count+lenrec
1686     locbyte=locbyte-loc_count+lenrec
1687     if(locbyte > lrecl .and. lastbuf) go to 900
1688     if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1689     end_block(irecs)=thisblock
1690     end_byte(irecs)=locbyte
1691     lensave=lenrec
1692     do i=1,4
1693      nextbyte=nextbyte+1_i_llong
1694      locbyte=locbyte+1_i_llong
1695      if(locbyte > lrecl .and. lastbuf) go to 900
1696      if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1697      lenrec4(i)=buf(locbyte)
1698     end do
1699     if(lenrec /= lensave) go to 890
1700
1701   end do
1702
1703 880  continue
1704      write(6,*)' reached impossible place in inventory_wrf_binary_file'
1705      close(in_unit)
1706      return
1707
1708 885  continue
1709      write(a_message,*)' problem in inventory_wrf_binary_file, lenrec has bad value before end of file'
1710      CALL wrf_message ( a_message )
1711      write(a_message,*)'     lenrec =',lenrec
1712      CALL wrf_message ( a_message )
1713      close(in_unit)
1714      return
1715
1716 890  continue
1717      write(a_message,*)' problem in inventory_wrf_binary_file, beginning and ending rec len words unequal'
1718      CALL wrf_message ( a_message )
1719      write(a_message,*)'     begining reclen =',lensave
1720      CALL wrf_message ( a_message )
1721      write(a_message,*)'       ending reclen =',lenrec
1722      CALL wrf_message ( a_message )
1723      write(a_message,*)'               irecs =',irecs
1724      CALL wrf_message ( a_message )
1725      write(a_message,*)'               nrecs =',nrecs
1726      CALL wrf_message ( a_message )
1727      call wrf_error_fatal("curious reclen discrepancy")
1728      close(in_unit)
1729      return
1730
1731 900  continue
1732      write(a_message,*)' normal end of file reached in inventory_wrf_binary_file'
1733      CALL wrf_message ( a_message )
1734      write(a_message,*)'        nblocks=',thisblock
1735      CALL wrf_message ( a_message )
1736      write(a_message,*)'          irecs,nrecs=',irecs,nrecs
1737      CALL wrf_message ( a_message )
1738      write(a_message,*)'         nreads=',nreads
1739      CALL wrf_message ( a_message )
1740      close(in_unit)
1741
1742 end subroutine inventory_wrf_binary_file
1743
1744 SUBROUTINE wrf_sizeof_integer( retval )
1745   IMPLICIT NONE
1746   INTEGER retval
1747 ! 4 is defined by CPP
1748   retval = 4
1749   RETURN
1750 END SUBROUTINE wrf_sizeof_integer
1751
1752 SUBROUTINE wrf_sizeof_real( retval )
1753   IMPLICIT NONE
1754   INTEGER retval
1755 ! 4 is defined by CPP
1756   retval = 4
1757   RETURN
1758 END SUBROUTINE wrf_sizeof_real
1759
1760 !!!
1761 !!!
1762 !!!
1763
1764 subroutine count_recs_wrf_binary_file(in_unit,wrf_ges_filename,nrecs)
1765 !$$$  subprogram documentation block
1766 !                .      .    .                                       .
1767 ! subprogram:    count_recs_binary_file  count # recs on wrf binary file
1768 !   prgmmr: parrish          org: np22                date: 2004-11-29
1769 !
1770 ! abstract: count number of sequential records contained in wrf binary
1771 !             file.  this is done by opening the file in direct access
1772 !             mode with block length of 2**20, the size of the physical
1773 !             blocks on ibm "blue" and "white" machines.  for optimal
1774 !             performance, change block length to correspond to the
1775 !             physical block length of host machine disk space.
1776 !             records are counted by looking for the 4 byte starting
1777 !             and ending sequential record markers, which contain the
1778 !             record size in bytes.  only blocks are read which are known
1779 !             by simple calculation to contain these record markers.
1780 !             even though this is done on one processor, it is still
1781 !             very fast, and the time will always scale by the number of
1782 !             sequential records, not their size.  this step and the
1783 !             following inventory step consistently take less than 0.1 seconds
1784 !             to complete.
1785 !
1786 ! program history log:
1787 !   2004-11-29  parrish
1788 !
1789 !   input argument list:
1790 !     in_unit          - fortran unit number where input file is opened through.
1791 !     wrf_ges_filename - filename of input wrf binary restart file
1792 !
1793 !   output argument list:
1794 !     nrecs            - number of sequential records found on input wrf binary restart fil
1795 !
1796 ! attributes:
1797 !   language: f90
1798 !   machine:  ibm RS/6000 SP
1799 !
1800 !$$$
1801
1802 !   do an initial read through of a wrf binary file, and get total number of sequential fil
1803
1804 !   use kinds, only: r_single,i_byte,i_long,i_llong
1805   implicit none
1806
1807   integer,intent(in)::in_unit
1808   character(*),intent(in)::wrf_ges_filename
1809   integer,intent(out)::nrecs
1810
1811   integer(i_llong) nextbyte,locbyte,thisblock
1812   integer(i_byte) lenrec4(4)
1813   integer(i_long) lenrec,lensave
1814   equivalence (lenrec4(1),lenrec)
1815   integer(i_byte) missing4(4)
1816   integer(i_long) missing
1817   equivalence (missing,missing4(1))
1818   integer(i_llong),parameter:: lrecl=2**20
1819   integer(i_byte) buf(lrecl)
1820   integer i,loc_count,nreads
1821   logical lastbuf
1822
1823   open(in_unit,file=trim(wrf_ges_filename),access='direct',recl=lrecl)
1824   nrecs=0
1825   missing=-9999
1826   nextbyte=0_i_llong
1827   locbyte=lrecl
1828   nreads=0
1829   lastbuf=.false.
1830   do
1831
1832 !   get length of next record
1833
1834     do i=1,4
1835      nextbyte=nextbyte+1_i_llong
1836      locbyte=locbyte+1_i_llong
1837      if(locbyte > lrecl .and. lastbuf) go to 900
1838      if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1839      lenrec4(i)=buf(locbyte)
1840     end do
1841     if(lenrec <= 0 .and. lastbuf) go to 900
1842     if(lenrec <= 0 .and. .not.lastbuf) go to 885
1843     nextbyte=nextbyte+1_i_llong
1844     locbyte=locbyte+1_i_llong
1845     if(locbyte > lrecl .and. lastbuf) go to 900
1846     if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1847
1848     nrecs=nrecs+1
1849
1850     loc_count=1
1851     do i=2,4
1852        if(loc_count.ge.lenrec) exit
1853        loc_count=loc_count+1
1854        nextbyte=nextbyte+1_i_llong
1855        locbyte=locbyte+1_i_llong
1856        if(locbyte > lrecl .and. lastbuf) go to 900
1857        if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1858     end do
1859     do i=1,4
1860        if(loc_count.ge.lenrec) exit
1861        loc_count=loc_count+1
1862        nextbyte=nextbyte+1_i_llong
1863        locbyte=locbyte+1_i_llong
1864        if(locbyte > lrecl .and. lastbuf) go to 900
1865        if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1866     end do
1867     nextbyte=nextbyte-loc_count+lenrec
1868     locbyte=locbyte-loc_count+lenrec
1869     if(locbyte > lrecl .and. lastbuf) go to 900
1870     if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1871     lensave=lenrec
1872     do i=1,4
1873      nextbyte=nextbyte+1_i_llong
1874      locbyte=locbyte+1_i_llong
1875      if(locbyte > lrecl .and. lastbuf) go to 900
1876      if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
1877      lenrec4(i)=buf(locbyte)
1878     end do
1879     if(lenrec /= lensave) go to 890
1880
1881   end do
1882
1883 880  continue
1884      write(6,*)' reached impossible place in count_recs_wrf_binary_file'
1885      close(in_unit)
1886      return
1887
1888 885  continue
1889      write(6,*)' problem in count_recs_wrf_binary_file, lenrec has bad value before end of file'
1890      write(6,*)'     lenrec =',lenrec
1891      close(in_unit)
1892      return
1893
1894 890  continue
1895      write(6,*)' problem in count_recs_wrf_binary_file, beginning and ending rec len words unequal'
1896      write(6,*)'     begining reclen =',lensave
1897      write(6,*)'       ending reclen =',lenrec
1898      close(in_unit)
1899      call wrf_error_fatal("bad reclen stuff")
1900      return
1901
1902 900  continue
1903      write(6,*)' normal end of file reached in count_recs_wrf_binary_file'
1904      write(6,*)'        nblocks=',thisblock
1905      write(6,*)'          nrecs=',nrecs
1906      write(6,*)'         nreads=',nreads
1907      close(in_unit)
1908
1909 end subroutine count_recs_wrf_binary_file
1910
1911 subroutine retrieve_field(in_unit,wrfges,out,start_block,end_block,start_byte,end_byte)
1912 !$$$  subprogram documentation block
1913 !                .      .    .                                       .
1914 ! subprogram:    retrieve_field  retrieve field from wrf binary file
1915 !   prgmmr: parrish          org: np22                date: 2004-11-29
1916 !
1917 ! abstract: still using direct access, retrieve a field from the wrf binary restart file.
1918 !
1919 ! program history log:
1920 !   2004-11-29  parrish
1921 !
1922 !   input argument list:
1923 !     in_unit          - fortran unit number where input file is opened through.
1924 !     wrfges - filename of input wrf binary restart file
1925 !     start_block      - direct access block number containing 1st byte of record
1926 !                            (after 4 byte record mark)
1927 !     end_block        - direct access block number containing last byte of record
1928 !                            (before 4 byte record mark)
1929 !     start_byte       - relative byte address in direct access block of 1st byte of record
1930 !     end_byte         - relative byte address in direct access block of last byte of record
1931 !
1932 !   output argument list:
1933 !     out              - output buffer where desired field is deposited
1934 !
1935 ! attributes:
1936 !   language: f90
1937 !   machine:  ibm RS/6000 SP
1938 !
1939 !$$$
1940
1941  ! use kinds, only: i_byte,i_kind
1942   implicit none
1943
1944   integer(i_kind),intent(in)::in_unit
1945   character(50),intent(in)::wrfges
1946   integer(i_kind),intent(in)::start_block,end_block,start_byte,end_byte
1947   integer(i_byte),intent(out)::out(*)
1948
1949   integer(i_llong),parameter:: lrecl=2**20
1950   integer(i_byte) buf(lrecl)
1951   integer(i_kind) i,ii,k
1952   integer(i_llong) ibegin,iend
1953   integer :: ierr
1954
1955   open(in_unit,file=trim(wrfges),access='direct',recl=lrecl)
1956
1957      write(6,*)' in retrieve_field, start_block,end_block=',start_block,end_block
1958      write(6,*)' in retrieve_field, start_byte,end_byte=',start_byte,end_byte
1959   ii=0
1960   do k=start_block,end_block
1961      read(in_unit,rec=k,iostat=ierr)buf
1962      ibegin=1 ; iend=lrecl
1963      if(k == start_block) ibegin=start_byte
1964      if(k == end_block) iend=end_byte
1965      do i=ibegin,iend
1966         ii=ii+1
1967         out(ii)=buf(i)
1968      end do
1969   end do
1970   close(in_unit)
1971
1972 end subroutine retrieve_field
1973
1974
1975 END MODULE module_wps_io_arw