main/ideal_nmm.F

   1 !  Create an initial data set for the WRF model based on an ideal condition
   2 !  This program is specifically set up for the NMM core.
   3
   4 PROGRAM ideal_nmm
   5
   6    USE module_machine
   7    USE module_domain
   8    USE module_initialize_ideal
   9    USE module_io_domain
  10    USE module_driver_constants
  11    USE module_configure
  12    USE module_timing
  13    USE module_check_a_mundo
  14 #if ( WRF_CHEM == 1 )
  15    USE module_input_chem_data
  16    USE module_input_chem_bioemiss
  17 #endif
  18    USE module_utility
  19 #ifdef DM_PARALLEL
  20    USE module_dm
  21 #endif
  22
  23    IMPLICIT NONE
  24
  25    REAL    :: time , bdyfrq
  26
  27    INTEGER :: loop , levels_to_process , debug_level
  28
  29
  30    TYPE(domain) , POINTER :: null_domain
  31    TYPE(domain) , POINTER :: grid
  32    TYPE (grid_config_rec_type)              :: config_flags
  33    INTEGER                :: number_at_same_level
  34
  35    INTEGER :: max_dom, domain_id
  36    INTEGER :: idum1, idum2
  37 #ifdef DM_PARALLEL
  38    INTEGER                 :: nbytes
  39 !   INTEGER, PARAMETER      :: configbuflen = 2*1024
  40    INTEGER, PARAMETER      :: configbuflen = 4*CONFIG_BUF_LEN
  41    INTEGER                 :: configbuf( configbuflen )
  42    LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
  43 #endif
  44
  45    INTEGER :: ids , ide , jds , jde , kds , kde
  46    INTEGER :: ims , ime , jms , jme , kms , kme
  47    INTEGER :: ips , ipe , jps , jpe , kps , kpe
  48    INTEGER :: ijds , ijde , spec_bdy_width
  49    INTEGER :: i , j , k , idts
  50
  51 #ifdef DEREF_KLUDGE
  52 !  see http://www.mmm.ucar.edu/wrf/WG2/topics/deref_kludge.htm
  53    INTEGER     :: sm31 , em31 , sm32 , em32 , sm33 , em33
  54    INTEGER     :: sm31x, em31x, sm32x, em32x, sm33x, em33x
  55    INTEGER     :: sm31y, em31y, sm32y, em32y, sm33y, em33y
  56 #endif
  57
  58    CHARACTER (LEN=80)     :: message
  59
  60    INTEGER :: start_year , start_month , start_day
  61    INTEGER :: start_hour , start_minute , start_second
  62    INTEGER :: end_year ,   end_month ,   end_day ,   &
  63               end_hour ,   end_minute ,   end_second
  64    INTEGER :: interval_seconds , real_data_init_type
  65    INTEGER :: time_loop_max , time_loop, rc
  66    REAL    :: t1,t2
  67
  68 #include "version_decl"
  69 #include "commit_decl"
  70
  71    INTERFACE
  72      SUBROUTINE Setup_Timekeeping( grid )
  73       USE module_domain
  74       TYPE(domain), POINTER :: grid
  75      END SUBROUTINE Setup_Timekeeping
  76    END INTERFACE
  77
  78    !  Define the name of this program (program_name defined in module_domain)
  79
  80    program_name = "IDEAL_HWRF " // TRIM(release_version) // " PREPROCESSOR"
  81
  82 #ifdef DM_PARALLEL
  83    CALL disable_quilting
  84 #endif
  85
  86 !       CALL start()
  87
  88    !  Initialize the modules used by the WRF system.
  89    !  Many of the CALLs made from the
  90    !  init_modules routine are NO-OPs.  Typical initializations
  91    !  are: the size of a
  92    !  REAL, setting the file handles to a pre-use value, defining moisture and
  93    !  chemistry indices, etc.
  94
  95    CALL       wrf_debug ( 100 , 'ideal_hwrf: calling init_modules ' )
  96
  97 !!!!   CALL init_modules
  98    CALL init_modules(1)   ! Phase 1 returns after MPI_INIT() (if it is called)
  99    CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN, rc=rc )
 100    CALL init_modules(2)   ! Phase 2 resumes after MPI_INIT() (if it is called)
 101
 102    !  The configuration switches mostly come from the NAMELIST input.
 103
 104 #ifdef DM_PARALLEL
 105    IF ( wrf_dm_on_monitor() ) THEN
 106       write(message,*) 'call initial_config'
 107       CALL wrf_message ( message )
 108       CALL initial_config
 109    ENDIF
 110    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 111    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 112    CALL set_config_as_buffer( configbuf, configbuflen )
 113    CALL wrf_dm_initialize
 114 #else
 115    CALL initial_config
 116 #endif
 117
 118    CALL check_nml_consistency
 119    CALL set_physics_rconfigs
 120
 121    CALL nl_get_debug_level ( 1, debug_level )
 122    CALL set_wrf_debug_level ( debug_level )
 123
 124    CALL  wrf_message ( program_name )
 125    CALL  wrf_message ( commit_version )
 126
 127    !  Allocate the space for the mother of all domains.
 128
 129    NULLIFY( null_domain )
 130    CALL  wrf_debug ( 100 , 'ideal_hwrf: calling alloc_and_configure_domain ' )
 131    CALL alloc_and_configure_domain ( domain_id  = 1           , &
 132                                      grid       = head_grid   , &
 133                                      parent     = null_domain , &
 134                                      kid        = -1            )
 135
 136    grid => head_grid
 137
 138 #include "deref_kludge.h"
 139    CALL Setup_Timekeeping ( grid )
 140    CALL domain_clock_set( grid, &
 141                           time_step_seconds=model_config_rec%interval_seconds )
 142    CALL wrf_debug ( 100 , 'ideal_hwrf: calling set_scalar_indices_from_config ' )
 143    CALL set_scalar_indices_from_config ( grid%id , idum1, idum2 )
 144
 145    CALL     wrf_debug ( 100 , 'ideal_hwrf: calling model_to_grid_config_rec ' )
 146
 147    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 148
 149    write(message,*) 'after model_to_grid_config_rec, e_we, e_sn are: ', &
 150                     config_flags%e_we, config_flags%e_sn
 151    CALL wrf_message(message)
 152
 153    !  Initialize the WRF IO: open files, init file handles, etc.
 154
 155    CALL       wrf_debug ( 100 , 'ideal_hwrf: calling init_wrfio' )
 156    CALL init_wrfio
 157
 158 !  Some of the configuration values may have been modified from the initial READ
 159 !  of the NAMELIST, so we re-broadcast the configuration records.
 160
 161 #ifdef DM_PARALLEL
 162    CALL wrf_debug ( 100 , 'ideal_hwrf: re-broadcast the configuration records' )
 163    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
 164    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
 165    CALL set_config_as_buffer( configbuf, configbuflen )
 166 #endif
 167
 168    !   No looping in this layer.
 169
 170    CALL med_sidata_input ( grid , config_flags )
 171
 172    !  We are done.
 173
 174    CALL       wrf_debug (   0 , 'ideal_hwrf: SUCCESS COMPLETE IDEAL_HWRF INIT' )
 175
 176 #ifdef DM_PARALLEL
 177     CALL wrf_dm_shutdown
 178 #endif
 179
 180    CALL WRFU_Finalize( rc=rc )
 181
 182 END PROGRAM ideal_nmm
 183
 184 SUBROUTINE med_sidata_input ( grid , config_flags )
 185   ! Driver layer
 186    USE module_domain
 187    USE module_io_domain
 188   ! Model layer
 189    USE module_configure
 190    USE module_bc_time_utilities
 191    USE module_initialize_ideal
 192    USE module_optional_input
 193 #if ( WRF_CHEM == 1 )
 194    USE module_input_chem_data
 195    USE module_input_chem_bioemiss
 196 #endif
 197
 198    USE module_si_io_nmm
 199
 200    USE module_date_time
 201
 202    IMPLICIT NONE
 203
 204
 205   ! Interface
 206    INTERFACE
 207      SUBROUTINE start_domain ( grid , allowed_to_read )
 208        USE module_domain
 209        TYPE (domain) grid
 210        LOGICAL, INTENT(IN) :: allowed_to_read
 211      END SUBROUTINE start_domain
 212    END INTERFACE
 213
 214   ! Arguments
 215    TYPE(domain)                :: grid
 216    TYPE (grid_config_rec_type) :: config_flags
 217   ! Local
 218    INTEGER                :: time_step_begin_restart
 219    INTEGER                :: idsi , ierr , myproc
 220    CHARACTER (LEN=256)     :: si_inpname
 221    CHARACTER (LEN=132)     :: message
 222
 223    CHARACTER(LEN=19) :: start_date_char , end_date_char , &
 224                         current_date_char , next_date_char
 225
 226    INTEGER :: time_loop_max , loop
 227    INTEGER :: julyr , julday , LEN
 228
 229    INTEGER :: io_form_auxinput1
 230    INTEGER, EXTERNAL :: use_package
 231
 232    LOGICAL :: using_binary_wrfsi
 233
 234    REAL :: gmt
 235    REAL :: t1,t2
 236
 237    INTEGER :: numx_sm_levels_input,numx_st_levels_input
 238    REAL,DIMENSION(100) :: smx_levels_input,stx_levels_input
 239
 240
 241 #ifdef DEREF_KLUDGE
 242 !  see http://www.mmm.ucar.edu/wrf/WG2/topics/deref_kludge.htm
 243    INTEGER     :: sm31 , em31 , sm32 , em32 , sm33 , em33
 244    INTEGER     :: sm31x, em31x, sm32x, em32x, sm33x, em33x
 245    INTEGER     :: sm31y, em31y, sm32y, em32y, sm33y, em33y
 246 #endif
 247
 248 #if ( HWRF == 1 )
 249   ! Sam Says:
 250
 251   ! The *INIT arrays are used to read init data written out by hwrf_prep_hybrid
 252   REAL,ALLOCATABLE,DIMENSION(:,:,:)::TINIT,UINIT,VINIT,QINIT,CWMINIT
 253   REAL,ALLOCATABLE,DIMENSION(:,:,:)::PINIT
 254   REAL,ALLOCATABLE,DIMENSION(:,:)::PDINIT
 255
 256   ! The *B arrays are used to read boundary data written out by hwrf_prep_hybrid
 257   REAL,ALLOCATABLE,DIMENSION(:,:,:)::TB,UB,VB,QB,CWMB
 258   REAL,ALLOCATABLE,DIMENSION(:,:)::PDB
 259
 260   INTEGER :: KB, LM, IM, JM, iunit_gfs, N
 261
 262   integer :: i,j,k
 263   LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
 264   integer :: ids,ide, jds,jde, kds,kde
 265   integer :: ims,ime, jms,jme, kms,kme
 266   integer :: its,ite, jts,jte, kts,kte
 267
 268   integer :: ioerror
 269 #endif
 270
 271 #include "deref_kludge.h"
 272
 273
 274    grid%input_from_file = .true.
 275    grid%input_from_file = .false.
 276
 277    CALL compute_si_start_and_end ( model_config_rec%start_year  (grid%id) , &
 278                                    model_config_rec%start_month (grid%id) , &
 279                                    model_config_rec%start_day   (grid%id) , &
 280                                    model_config_rec%start_hour  (grid%id) , &
 281                                    model_config_rec%start_minute(grid%id) , &
 282                                    model_config_rec%start_second(grid%id) , &
 283                                    model_config_rec%  end_year  (grid%id) , &
 284                                    model_config_rec%  end_month (grid%id) , &
 285                                    model_config_rec%  end_day   (grid%id) , &
 286                                    model_config_rec%  end_hour  (grid%id) , &
 287                                    model_config_rec%  end_minute(grid%id) , &
 288                                    model_config_rec%  end_second(grid%id) , &
 289                                    model_config_rec%interval_seconds      , &
 290                                    model_config_rec%real_data_init_type   , &
 291                                    start_date_char , end_date_char , time_loop_max )
 292
 293    !  Here we define the initial time to process, for later use by the code.
 294
 295    current_date_char = start_date_char
 296 !   start_date = start_date_char // '.0000'
 297    start_date = start_date_char
 298    current_date = start_date
 299
 300    CALL nl_set_bdyfrq ( grid%id , REAL(model_config_rec%interval_seconds) )
 301
 302    !  Loop over each time period to process.
 303
 304    write(message,*) 'time_loop_max: ', time_loop_max
 305    CALL wrf_message(message)
 306    DO loop = 1 , time_loop_max
 307
 308      internal_time_loop=loop
 309
 310       write(message,*) 'loop=', loop
 311       CALL wrf_message(message)
 312
 313       write(message,*) '-----------------------------------------------------------'
 314       CALL wrf_message(message)
 315
 316       write(message,*) ' '
 317       CALL wrf_message(message)
 318       write(message,'(A,A,A,I2,A,I2)') ' Current date being processed: ', &
 319         current_date, ', which is loop #',loop,' out of ',time_loop_max
 320       CALL wrf_message(message)
 321
 322       !  After current_date has been set, fill in the julgmt stuff.
 323
 324       CALL geth_julgmt ( config_flags%julyr , config_flags%julday , &
 325                                               config_flags%gmt )
 326
 327       !  Now that the specific Julian info is available,
 328       !  save these in the model config record.
 329
 330       CALL nl_set_gmt (grid%id, config_flags%gmt)
 331       CALL nl_set_julyr (grid%id, config_flags%julyr)
 332       CALL nl_set_julday (grid%id, config_flags%julday)
 333
 334       CALL nl_get_io_form_auxinput1( 1, io_form_auxinput1 )
 335       using_binary_wrfsi=.false.
 336
 337
 338       write(message,*) 'TRIM(config_flags%auxinput1_inname): ', TRIM(config_flags%auxinput1_inname)
 339       CALL wrf_message(message)
 340
 341 #if ( HWRF == 1 )
 342      ifph_onlyfirst: if(.not.grid%use_prep_hybrid .or. loop==1) then
 343 #endif
 344       IF (config_flags%auxinput1_inname(1:10) .eq. 'real_input') THEN
 345          using_binary_wrfsi=.true.
 346       ENDIF
 347
 348       SELECT CASE ( use_package(io_form_auxinput1) )
 349 #if defined(NETCDF)
 350       CASE ( IO_NETCDF )
 351
 352       !  Open the wrfinput file.
 353
 354         current_date_char(11:11)='_'
 355
 356        WRITE ( wrf_err_message , FMT='(A,A)' )'med_sidata_input: calling open_r_dataset for ',TRIM(config_flags%auxinput1_inname)
 357        CALL wrf_debug ( 100 , wrf_err_message )
 358        IF ( config_flags%auxinput1_inname(1:8) .NE. 'wrf_real' ) THEN
 359           CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , current_date_char , &
 360                                      config_flags%io_form_auxinput1 )
 361        ELSE
 362           CALL construct_filename2a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , current_date_char )
 363        END IF
 364        CALL open_r_dataset ( idsi, TRIM(si_inpname) , grid , config_flags , "DATASET=AUXINPUT1", ierr )
 365
 366        IF ( ierr .NE. 0 ) THEN
 367           CALL wrf_error_fatal( 'error opening ' // TRIM(si_inpname) // ' for input; bad date in namelist or file not in directory' )
 368        ENDIF
 369
 370       !  Input data.
 371
 372       CALL wrf_debug (100, 'med_sidata_input: call input_auxinput1_wrf')
 373
 374       CALL input_auxinput1 ( idsi, grid, config_flags, ierr )
 375
 376       !  Possible optional SI input.  This sets flags used by init_domain.
 377
 378       IF ( loop .EQ. 1 ) THEN
 379          CALL  wrf_debug (100, 'med_sidata_input: call init_module_optional_input' )
 380          CALL init_module_optional_input ( grid , config_flags )
 381       CALL wrf_debug ( 100 , 'med_sidata_input: calling optional_input' )
 382 !
 383       CALL optional_input ( grid , idsi , config_flags )
 384         write(0,*) 'maxval st_input(1) within ideal_hwrf: ', maxval(st_input(:,1,:))
 385       END IF
 386 !
 387       CALL close_dataset ( idsi , config_flags , "DATASET=AUXINPUT1" )
 388
 389 #endif
 390 #ifdef INTIO
 391       CASE ( IO_INTIO )
 392
 393       !  Possible optional SI input.  This sets flags used by init_domain.
 394
 395       IF ( loop .EQ. 1 ) THEN
 396          CALL  wrf_debug (100, 'med_sidata_input: call init_module_optional_input' )
 397          CALL init_module_optional_input ( grid , config_flags )
 398       END IF
 399
 400       IF (using_binary_wrfsi) THEN
 401
 402         current_date_char(11:11)='_'
 403         CALL read_si ( grid, current_date_char )
 404         current_date_char(11:11)='T'
 405
 406       ELSE
 407
 408         write(message,*) 'binary WPS branch'
 409         CALL wrf_message(message)
 410         current_date_char(11:11)='_'
 411         CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , current_date_char , &
 412                                      config_flags%io_form_auxinput1 )
 413         CALL read_wps ( grid, trim(si_inpname), current_date_char, config_flags%num_metgrid_levels )
 414 !!! bogus set some flags??
 415       flag_metgrid=1
 416       flag_soilhgt=1
 417
 418
 419           ENDIF
 420
 421 #endif
 422       CASE DEFAULT
 423         CALL wrf_error_fatal('ideal_hwrf: not valid io_form_auxinput1')
 424       END SELECT
 425 #if ( HWRF == 1 )
 426      endif ifph_onlyfirst
 427 #endif
 428
 429       grid%islope=1
 430       grid%vegfra=grid%vegfrc
 431       grid%dfrlg=grid%dfl/9.81
 432
 433       grid%isurban=1
 434       grid%isoilwater=14
 435
 436       !  Initialize the mother domain for this time period with input data.
 437
 438       CALL wrf_debug ( 100 , 'med_sidata_input: calling init_domain' )
 439       grid%input_from_file = .true.
 440
 441       CALL init_domain ( grid )
 442
 443 #if ( HWRF == 1 )
 444      read_phinit: if(grid%use_prep_hybrid) then
 445 #if defined(DM_PARALLEL)
 446         if(.not. wrf_dm_on_monitor()) then
 447            call wrf_error_fatal('ideal_hwrf: in use_prep_hybrid mode, threading and mpi are forbidden.')
 448         endif
 449 #endif
 450
 451         ph_loop1: if(loop==1) then
 452
 453            ! determine kds, ids, jds
 454            SELECT CASE ( model_data_order )
 455            CASE ( DATA_ORDER_ZXY )
 456               kds = grid%sd31 ; kde = grid%ed31 ;
 457               ids = grid%sd32 ; ide = grid%ed32 ;
 458               jds = grid%sd33 ; jde = grid%ed33 ;
 459
 460               kms = grid%sm31 ; kme = grid%em31 ;
 461               ims = grid%sm32 ; ime = grid%em32 ;
 462               jms = grid%sm33 ; jme = grid%em33 ;
 463
 464               kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
 465               its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
 466               jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
 467
 468            CASE ( DATA_ORDER_XYZ )
 469               ids = grid%sd31 ; ide = grid%ed31 ;
 470               jds = grid%sd32 ; jde = grid%ed32 ;
 471               kds = grid%sd33 ; kde = grid%ed33 ;
 472
 473               ims = grid%sm31 ; ime = grid%em31 ;
 474               jms = grid%sm32 ; jme = grid%em32 ;
 475               kms = grid%sm33 ; kme = grid%em33 ;
 476
 477               its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
 478               jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
 479               kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
 480
 481            CASE ( DATA_ORDER_XZY )
 482               ids = grid%sd31 ; ide = grid%ed31 ;
 483               kds = grid%sd32 ; kde = grid%ed32 ;
 484               jds = grid%sd33 ; jde = grid%ed33 ;
 485
 486               ims = grid%sm31 ; ime = grid%em31 ;
 487               kms = grid%sm32 ; kme = grid%em32 ;
 488               jms = grid%sm33 ; jme = grid%em33 ;
 489
 490               its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
 491               kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
 492               jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
 493
 494            END SELECT
 495            ! Allocate 3D initialization arrays:
 496            call wrf_message('ALLOCATE PREP_HYBRID INIT ARRAYS')
 497            ALLOCATE ( TINIT  (ids:(ide-1),kds:(kde-1)  ,jds:(jde-1)) )
 498            ALLOCATE ( PINIT  (ids:(ide-1),kds:kde      ,jds:(jde-1)) )
 499            ALLOCATE ( UINIT  (ids:(ide-1),kds:(kde-1)  ,jds:(jde-1)) )
 500            ALLOCATE ( VINIT  (ids:(ide-1),kds:(kde-1)  ,jds:(jde-1)) )
 501            ALLOCATE ( QINIT  (ids:(ide-1),kds:(kde-1)  ,jds:(jde-1)) )
 502            ALLOCATE ( CWMINIT(ids:(ide-1),kds:(kde-1)  ,jds:(jde-1)) )
 503            ALLOCATE ( PDINIT (ids:(ide-1),              jds:(jde-1)) )
 504
 505            REWIND 900
 506            READ(900,iostat=ioerror) PDINIT,TINIT,QINIT,CWMINIT,UINIT,VINIT,PINIT
 507            if(ioerror/=0) then
 508               call wrf_error_fatal('Unable to read MAKBND output from unit 900.')
 509            endif
 510            WRITE(0,*) 'U V T AT 10 10 10 ',UINIT(10,10,10),VINIT(10,10,10),TINIT(10,10,10)
 511            ! Switch from IKJ to IJK ordering
 512            DO I = ids,ide-1
 513               DO J = jds,jde-1
 514                  grid%pd(I,J) = PDINIT(I,J)
 515                  DO K = kds,kde-1
 516                     grid%q2(I,J,K) = 0
 517                     grid%u(I,J,K) = UINIT(I,K,J)
 518                     grid%v(I,J,K) = VINIT(I,K,J)
 519                     grid%t(I,J,K) = TINIT(I,K,J)
 520                     grid%q(I,J,K) = QINIT(I,K,J)
 521                     grid%cwm(I,J,K) = CWMINIT(I,K,J)
 522                  ENDDO
 523                  !  Was commented out in original V2 HWRF too:
 524                  !      DO K = kds,kde
 525                  !         grid%nmm_pint(I,J,K) = pinit(I,K,J)
 526                  !      ENDDO
 527               ENDDO
 528            ENDDO
 529
 530            call wrf_message('DEALLOCATE PREP_HYBRID INIT ARRAYS')
 531            deallocate(TINIT,PINIT,UINIT,VINIT,QINIT,CWMINIT,PDINIT)
 532         end if ph_loop1
 533      end if read_phinit
 534 #endif
 535
 536       CALL model_to_grid_config_rec ( grid%id, model_config_rec, config_flags )
 537
 538       !  Close this file that is output from the SI and input to this pre-proc.
 539
 540       CALL wrf_debug ( 100 , 'med_sidata_input: back from init_domain' )
 541
 542
 543 !!! not sure about this, but doesnt seem like needs to be called each time
 544       IF ( loop .EQ. 1 ) THEN
 545         CALL start_domain ( grid , .TRUE.)
 546       END IF
 547
 548 #if ( WRF_CHEM == 1 )
 549       IF ( loop == 1 ) THEN
 550 !        IF ( ( grid%chem_opt .EQ. RADM2     ) .OR. &
 551 !             ( grid%chem_opt .EQ. RADM2SORG ) .OR. &
 552 !             ( grid%chem_opt .EQ. RACM      ) .OR. &
 553 !             ( grid%chem_opt .EQ. RACMSORG  ) ) THEN
 554          IF( grid%chem_opt > 0 ) then
 555            ! Read the chemistry data from a previous wrf forecast (wrfout file)
 556            IF(grid%chem_in_opt == 1 ) THEN
 557               message = 'INITIALIZING CHEMISTRY WITH OLD SIMULATION'
 558               CALL  wrf_message ( message )
 559
 560               CALL input_ext_chem_file( grid )
 561
 562               IF(grid%bio_emiss_opt == BEIS311 ) THEN
 563                  message = 'READING BEIS3.11 EMISSIONS DATA'
 564                  CALL  wrf_message ( message )
 565                  CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
 566               else IF(grid%bio_emiss_opt == 3 ) THEN !shc
 567                  message = 'READING MEGAN 2 EMISSIONS DATA'
 568                  CALL  wrf_message ( message )
 569                  CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
 570               END IF
 571
 572            ELSEIF(grid%chem_in_opt == 0)then
 573               ! Generate chemistry data from a idealized vertical profile
 574               message = 'STARTING WITH BACKGROUND CHEMISTRY '
 575               CALL  wrf_message ( message )
 576
 577               write(message,*)' ETA1 '
 578               CALL  wrf_message ( message )
 579 !             write(message,*) grid%eta1
 580 !             CALL  wrf_message ( message )
 581
 582               CALL input_chem_profile ( grid )
 583
 584               IF(grid%bio_emiss_opt == BEIS311 ) THEN
 585                  message = 'READING BEIS3.11 EMISSIONS DATA'
 586                  CALL  wrf_message ( message )
 587                  CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
 588               else IF(grid%bio_emiss_opt == 3 ) THEN !shc
 589                  message = 'READING MEGAN 2 EMISSIONS DATA'
 590                  CALL  wrf_message ( message )
 591                  CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
 592               END IF
 593
 594            ELSE
 595              message = 'RUNNING WITHOUT CHEMISTRY INITIALIZATION'
 596              CALL  wrf_message ( message )
 597            ENDIF
 598          ENDIF
 599       ENDIF
 600 #endif
 601
 602       config_flags%isurban=1
 603       config_flags%isoilwater=14
 604
 605       CALL assemble_output ( grid , config_flags , loop , time_loop_max )
 606
 607       !  Here we define the next time that we are going to process.
 608
 609       CALL geth_newdate ( current_date_char , start_date_char , &
 610                           loop * model_config_rec%interval_seconds )
 611       current_date =  current_date_char // '.0000'
 612
 613       CALL domain_clock_set( grid, current_date(1:19) )
 614
 615       write(message,*) 'current_date= ', current_date
 616       CALL wrf_message(message)
 617
 618    END DO
 619 END SUBROUTINE med_sidata_input
 620
 621 SUBROUTINE compute_si_start_and_end (  &
 622           start_year, start_month, start_day, start_hour, &
 623           start_minute, start_second, &
 624           end_year ,   end_month ,   end_day ,   end_hour , &
 625           end_minute ,   end_second , &
 626           interval_seconds , real_data_init_type , &
 627           start_date_char , end_date_char , time_loop_max )
 628
 629    USE module_date_time
 630
 631    IMPLICIT NONE
 632
 633    INTEGER :: start_year , start_month , start_day , &
 634               start_hour , start_minute , start_second
 635    INTEGER ::   end_year ,   end_month ,   end_day , &
 636                 end_hour ,   end_minute ,   end_second
 637    INTEGER :: interval_seconds , real_data_init_type
 638    INTEGER :: time_loop_max , time_loop
 639
 640    CHARACTER(LEN=132) :: message
 641    CHARACTER(LEN=19)  :: current_date_char , start_date_char , &
 642                         end_date_char , next_date_char
 643
 644 !   WRITE ( start_date_char , FMT = &
 645 !         '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
 646 !         start_year,start_month,start_day,start_hour,start_minute,start_second
 647 !   WRITE (   end_date_char , FMT = &
 648 !         '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
 649 !          end_year,  end_month,  end_day,  end_hour,  end_minute,  end_second
 650
 651    WRITE ( start_date_char , FMT = &
 652          '(I4.4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2)' ) &
 653          start_year,start_month,start_day,start_hour,start_minute,start_second
 654    WRITE (   end_date_char , FMT = &
 655          '(I4.4,"-",I2.2,"-",I2.2,"T",I2.2,":",I2.2,":",I2.2)' ) &
 656           end_year,  end_month,  end_day,  end_hour,  end_minute,  end_second
 657
 658 !  start_date = start_date_char // '.0000'
 659
 660    !  Figure out our loop count for the processing times.
 661
 662    time_loop = 1
 663    PRINT '(A,I4,A,A,A)','Time period #',time_loop, &
 664                         ' to process = ',start_date_char,'.'
 665    current_date_char = start_date_char
 666    loop_count : DO
 667       CALL geth_newdate (next_date_char, current_date_char, interval_seconds )
 668       IF      ( next_date_char .LT. end_date_char ) THEN
 669          time_loop = time_loop + 1
 670          PRINT '(A,I4,A,A,A)','Time period #',time_loop,&
 671                               ' to process = ',next_date_char,'.'
 672          current_date_char = next_date_char
 673       ELSE IF ( next_date_char .EQ. end_date_char ) THEN
 674          time_loop = time_loop + 1
 675          PRINT '(A,I4,A,A,A)','Time period #',time_loop,&
 676                               ' to process = ',next_date_char,'.'
 677          PRINT '(A,I4,A)','Total analysis times to input = ',time_loop,'.'
 678          time_loop_max = time_loop
 679          EXIT loop_count
 680       ELSE IF ( next_date_char .GT. end_date_char ) THEN
 681          PRINT '(A,I4,A)','Total analysis times to input = ',time_loop,'.'
 682          time_loop_max = time_loop
 683          EXIT loop_count
 684       END IF
 685    END DO loop_count
 686         write(message,*) 'done in si_start_and_end'
 687         CALL wrf_message(message)
 688 END SUBROUTINE compute_si_start_and_end
 689
 690 SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max )
 691
 692 !!! replace with something?   USE module_big_step_utilities_em
 693
 694    USE module_domain
 695    USE module_io_domain
 696    USE module_configure
 697    USE module_date_time
 698    USE module_bc
 699    IMPLICIT NONE
 700
 701 #if ( HWRF == 1 )
 702   external get_wrf_debug_level
 703   integer :: debug
 704 #endif
 705
 706    TYPE(domain)                 :: grid
 707    TYPE (grid_config_rec_type)  :: config_flags
 708    INTEGER , INTENT(IN)         :: loop , time_loop_max
 709
 710    INTEGER :: ids , ide , jds , jde , kds , kde
 711    INTEGER :: ims , ime , jms , jme , kms , kme
 712    INTEGER :: ips , ipe , jps , jpe , kps , kpe
 713    INTEGER :: ijds , ijde , spec_bdy_width
 714    INTEGER :: inc_h,inc_v
 715    INTEGER :: i , j , k , idts
 716
 717    INTEGER :: id1 , interval_seconds , ierr, rc, sst_update
 718    INTEGER , SAVE :: id ,id4
 719    CHARACTER (LEN=80) :: inpname , bdyname
 720    CHARACTER(LEN= 4) :: loop_char
 721    CHARACTER(LEN=132) :: message
 722 character *19 :: temp19
 723 character *24 :: temp24 , temp24b
 724
 725    REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ubdy3dtemp1 , vbdy3dtemp1 ,&
 726                                                 tbdy3dtemp1 , &
 727                                                 cwmbdy3dtemp1 , qbdy3dtemp1,&
 728                                                 q2bdy3dtemp1 , pdbdy2dtemp1
 729    REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ubdy3dtemp2 , vbdy3dtemp2 , &
 730                                                 tbdy3dtemp2 , &
 731                                                 cwmbdy3dtemp2 , qbdy3dtemp2, &
 732                                                 q2bdy3dtemp2, pdbdy2dtemp2
 733    REAL :: t1,t2
 734
 735 #ifdef DEREF_KLUDGE
 736 !  see http://www.mmm.ucar.edu/wrf/WG2/topics/deref_kludge.htm
 737    INTEGER     :: sm31 , em31 , sm32 , em32 , sm33 , em33
 738    INTEGER     :: sm31x, em31x, sm32x, em32x, sm33x, em33x
 739    INTEGER     :: sm31y, em31y, sm32y, em32y, sm33y, em33y
 740 #endif
 741
 742 #if ( HWRF == 1 )
 743   ! Sam says:
 744
 745   ! The *B arrays are used to read boundary data written out by hwrf_prep_hybrid
 746   REAL,ALLOCATABLE,DIMENSION(:,:,:)::TB,UB,VB,QB,CWMB
 747   REAL,ALLOCATABLE,DIMENSION(:,:)::PDB
 748
 749   ! Dimensions and looping variables:
 750   INTEGER :: KB, LM, IM, JM, N
 751
 752   ! Unit number to read boundary data from (changes each time)
 753   INTEGER :: iunit_gfs
 754
 755   ! Did we allocate the prep_hybrid input arrays?
 756   LOGICAL :: alloc_ph_arrays
 757
 758   integer :: ioerror
 759 #endif
 760
 761 #include "deref_kludge.h"
 762
 763 #if ( HWRF == 1 )
 764   alloc_ph_arrays=.false.
 765   call get_wrf_debug_level(debug)
 766 #endif
 767
 768    !  Various sizes that we need to be concerned about.
 769
 770    ids = grid%sd31
 771    ide = grid%ed31-1 ! 030730tst
 772    jds = grid%sd32
 773    jde = grid%ed32-1 ! 030730tst
 774    kds = grid%sd33
 775    kde = grid%ed33-1 ! 030730tst
 776
 777    ims = grid%sm31
 778    ime = grid%em31
 779    jms = grid%sm32
 780    jme = grid%em32
 781    kms = grid%sm33
 782    kme = grid%em33
 783
 784    ips = grid%sp31
 785    ipe = grid%ep31-1 ! 030730tst
 786    jps = grid%sp32
 787    jpe = grid%ep32-1 ! 030730tst
 788    kps = grid%sp33
 789    kpe = grid%ep33-1 ! 030730tst
 790
 791         if (IPE .ne. IDE) IPE=IPE+1
 792         if (JPE .ne. JDE) JPE=JPE+1
 793
 794         write(message,*) 'assemble output (ids,ide): ', ids,ide
 795         CALL wrf_message(message)
 796         write(message,*) 'assemble output (ims,ime): ', ims,ime
 797         CALL wrf_message(message)
 798         write(message,*) 'assemble output (ips,ipe): ', ips,ipe
 799         CALL wrf_message(message)
 800
 801         write(message,*) 'assemble output (jds,jde): ', jds,jde
 802         CALL wrf_message(message)
 803         write(message,*) 'assemble output (jms,jme): ', jms,jme
 804         CALL wrf_message(message)
 805         write(message,*) 'assemble output (jps,jpe): ', jps,jpe
 806         CALL wrf_message(message)
 807
 808         write(message,*) 'assemble output (kds,kde): ', kds,kde
 809         CALL wrf_message(message)
 810         write(message,*) 'assemble output (kms,kme): ', kms,kme
 811         CALL wrf_message(message)
 812         write(message,*) 'assemble output (kps,kpe): ', kps,kpe
 813         CALL wrf_message(message)
 814
 815    ijds = MIN ( ids , jds )
 816 !mptest030805   ijde = MAX ( ide , jde )
 817    ijde = MAX ( ide , jde ) + 1   ! to make stuff_bdy dimensions consistent with alloc
 818
 819    !  Boundary width, scalar value.
 820
 821    spec_bdy_width = model_config_rec%spec_bdy_width
 822    interval_seconds = model_config_rec%interval_seconds
 823    sst_update = model_config_rec%sst_update
 824
 825 !-----------------------------------------------------------------------
 826 !
 827    main_loop_test: IF ( loop .EQ. 1 ) THEN
 828 !
 829 !-----------------------------------------------------------------------
 830
 831       IF ( time_loop_max .NE. 1 ) THEN
 832          IF(sst_update .EQ. 1)THEN
 833            CALL construct_filename1( inpname , 'wrflowinp' , grid%id , 2 )
 834            CALL open_w_dataset ( id4, TRIM(inpname) , grid , config_flags , output_auxinput4 , "DATASET=AUXINPUT4", ierr )
 835            IF ( ierr .NE. 0 ) THEN
 836               CALL wrf_error_fatal( 'ideal_hwrf: error opening wrflowinp for writing' )
 837            END IF
 838            CALL output_auxinput4 ( id4, grid , config_flags , ierr )
 839          END IF
 840       END IF
 841
 842
 843    !  This is the space needed to save the current 3d data for use in computing
 844    !  the lateral boundary tendencies.
 845
 846       ALLOCATE ( ubdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 847       ALLOCATE ( vbdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 848       ALLOCATE ( tbdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 849       ALLOCATE ( qbdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 850       ALLOCATE ( cwmbdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 851       ALLOCATE ( q2bdy3dtemp1(ims:ime,jms:jme,kms:kme) )
 852       ALLOCATE ( pdbdy2dtemp1(ims:ime,jms:jme,1:1) )
 853
 854         ubdy3dtemp1=0.
 855         vbdy3dtemp1=0.
 856         tbdy3dtemp1=0.
 857         qbdy3dtemp1=0.
 858         cwmbdy3dtemp1=0.
 859         q2bdy3dtemp1=0.
 860         pdbdy2dtemp1=0.
 861
 862       ALLOCATE ( ubdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 863       ALLOCATE ( vbdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 864       ALLOCATE ( tbdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 865       ALLOCATE ( qbdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 866       ALLOCATE ( cwmbdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 867       ALLOCATE ( q2bdy3dtemp2(ims:ime,jms:jme,kms:kme) )
 868       ALLOCATE ( pdbdy2dtemp2(ims:ime,jms:jme,1:1) )
 869
 870         ubdy3dtemp2=0.
 871         vbdy3dtemp2=0.
 872         tbdy3dtemp2=0.
 873         qbdy3dtemp2=0.
 874         cwmbdy3dtemp2=0.
 875         q2bdy3dtemp2=0.
 876         pdbdy2dtemp2=0.
 877
 878       !  Open the wrfinput file.  From this program, this is an *output* file.
 879
 880       CALL construct_filename1( inpname , 'wrfinput' , grid%id , 2 )
 881
 882       CALL open_w_dataset ( id1, TRIM(inpname) , grid , config_flags , &
 883                             output_input , "DATASET=INPUT", ierr )
 884
 885       IF ( ierr .NE. 0 ) THEN
 886       CALL wrf_error_fatal( 'ideal_hwrf: error opening wrfinput for writing' )
 887       ENDIF
 888
 889 !     CALL calc_current_date ( grid%id , 0. )
 890 !      grid%write_metadata = .true.
 891
 892         write(message,*) 'making call to output_input'
 893         CALL wrf_message(message)
 894
 895         CALL output_input ( id1, grid , config_flags , ierr )
 896
 897 !***
 898 !***  CLOSE THE WRFINPUT DATASET
 899 !***
 900       CALL close_dataset ( id1 , config_flags , "DATASET=INPUT" )
 901
 902       !  We need to save the 3d data to compute a
 903       !  difference during the next loop.
 904
 905 !
 906 !-----------------------------------------------------------------------
 907 !***  SOUTHERN BOUNDARY
 908 !-----------------------------------------------------------------------
 909 !
 910
 911         IF(JPS==JDS)THEN
 912           J=1
 913           DO k = kps , MIN(kde,kpe)
 914           DO i = ips , MIN(ide,ipe)
 915             ubdy3dtemp1(i,j,k) = grid%u(i,j,k)
 916             vbdy3dtemp1(i,j,k) = grid%v(i,j,k)
 917             tbdy3dtemp1(i,j,k) = grid%t(i,j,k)
 918             qbdy3dtemp1(i,j,k) = grid%q(i,j,k)
 919             cwmbdy3dtemp1(i,j,k) = grid%cwm(i,j,k)
 920             q2bdy3dtemp1(i,j,k) = grid%q2(i,j,k)
 921           END DO
 922           END DO
 923
 924           DO i = ips , MIN(ide,ipe)
 925             pdbdy2dtemp1(i,j,1) = grid%pd(i,j)
 926           END DO
 927         ENDIF
 928
 929 !
 930 !-----------------------------------------------------------------------
 931 !***  NORTHERN BOUNDARY
 932 !-----------------------------------------------------------------------
 933 !
 934         IF(JPE==JDE)THEN
 935           J=MIN(JDE,JPE)
 936           DO k = kps , MIN(kde,kpe)
 937           DO i = ips , MIN(ide,ipe)
 938             ubdy3dtemp1(i,j,k) = grid%u(i,j,k)
 939             vbdy3dtemp1(i,j,k) = grid%v(i,j,k)
 940             tbdy3dtemp1(i,j,k) = grid%t(i,j,k)
 941             qbdy3dtemp1(i,j,k) = grid%q(i,j,k)
 942             cwmbdy3dtemp1(i,j,k) = grid%cwm(i,j,k)
 943             q2bdy3dtemp1(i,j,k) = grid%q2(i,j,k)
 944           END DO
 945           END DO
 946
 947           DO i = ips , MIN(ide,ipe)
 948             pdbdy2dtemp1(i,j,1) = grid%pd(i,j)
 949           END DO
 950         ENDIF
 951
 952 !
 953 !-----------------------------------------------------------------------
 954 !***  WESTERN BOUNDARY
 955 !-----------------------------------------------------------------------
 956 !
 957         write(message,*) 'western boundary, store winds over J: ', jps, min(jpe,jde)
 958         CALL wrf_message(message)
 959
 960         IF(IPS==IDS)THEN
 961           I=1
 962           DO k = kps , MIN(kde,kpe)
 963           inc_h=mod(jps+1,2)
 964           DO j = jps+inc_h, min(jde,jpe),2
 965
 966         if (J .ge. 3 .and. J .le. JDE-2 .and. mod(J,2) .eq. 1) then
 967             tbdy3dtemp1(i,j,k) = grid%t(i,j,k)
 968             qbdy3dtemp1(i,j,k) = grid%q(i,j,k)
 969             cwmbdy3dtemp1(i,j,k) = grid%cwm(i,j,k)
 970             q2bdy3dtemp1(i,j,k) = grid%q2(i,j,k)
 971       if(k==1)then
 972         write(message,*)' loop=',loop,' i=',i,' j=',j,' tbdy3dtemp1(i,j,k)=',tbdy3dtemp1(i,j,k)
 973         CALL wrf_debug(10,message)
 974       endif
 975         endif
 976           END DO
 977           END DO
 978
 979           DO k = kps , MIN(kde,kpe)
 980           inc_v=mod(jps,2)
 981           DO j = jps+inc_v, min(jde,jpe),2
 982         if (J .ge. 2 .and. J .le. JDE-1 .and. mod(J,2) .eq. 0) then
 983             ubdy3dtemp1(i,j,k) = grid%u(i,j,k)
 984             vbdy3dtemp1(i,j,k) = grid%v(i,j,k)
 985         endif
 986           END DO
 987           END DO
 988 !
 989           inc_h=mod(jps+1,2)
 990         DO j = jps+inc_h, min(jde,jpe),2
 991         if (J .ge. 3 .and. J .le. JDE-2 .and. mod(J,2) .eq. 1) then
 992             pdbdy2dtemp1(i,j,1) = grid%pd(i,j)
 993           write(message,*)' loop=',loop,' i=',i,' j=',j,' pdbdy2dtemp1(i,j)=',pdbdy2dtemp1(i,j,1)
 994           CALL wrf_debug(10,message)
 995         endif
 996           END DO
 997         ENDIF
 998 !
 999 !-----------------------------------------------------------------------
1000 !***  EASTERN BOUNDARY
1001 !-----------------------------------------------------------------------
1002 !
1003         IF(IPE==IDE)THEN
1004           I=MIN(IDE,IPE)
1005 !
1006           DO k = kps , MIN(kde,kpe)
1007 !
1008 !***   Make sure the J loop is on the global boundary
1009 !
1010           inc_h=mod(jps+1,2)
1011           DO j = jps+inc_h, min(jde,jpe),2
1012         if (J .ge. 3 .and. J .le. JDE-2 .and. mod(J,2) .eq. 1) then
1013             tbdy3dtemp1(i,j,k) = grid%t(i,j,k)
1014             qbdy3dtemp1(i,j,k) = grid%q(i,j,k)
1015             cwmbdy3dtemp1(i,j,k) = grid%cwm(i,j,k)
1016             q2bdy3dtemp1(i,j,k) = grid%q2(i,j,k)
1017         endif
1018           END DO
1019           END DO
1020
1021           DO k = kps , MIN(kde,kpe)
1022           inc_v=mod(jps,2)
1023           DO j = jps+inc_v, min(jde,jpe),2
1024         if (J .ge. 2 .and. J .le. JDE-1 .and. mod(J,2) .eq. 0) then
1025             ubdy3dtemp1(i,j,k) = grid%u(i,j,k)
1026             vbdy3dtemp1(i,j,k) = grid%v(i,j,k)
1027         endif
1028           END DO
1029           END DO
1030 !
1031           inc_h=mod(jps+1,2)
1032           DO j = jps+inc_h, min(jde,jpe),2
1033         if (J .ge. 3 .and. J .le. JDE-2 .and. mod(J,2) .eq. 1) then
1034             pdbdy2dtemp1(i,j,1) = grid%pd(i,j)
1035         endif
1036           END DO
1037         ENDIF
1038
1039
1040       !  There are 2 components to the lateral boundaries.
1041       !  First, there is the starting
1042       !  point of this time period - just the outer few rows and columns.
1043
1044
1045  CALL stuff_bdy_ijk (ubdy3dtemp1, grid%u_bxs, grid%u_bxe, &
1046                                   grid%u_bys, grid%u_bye, &
1047                                   'N', spec_bdy_width  , &
1048                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1049                                   ims , ime , jms , jme , kms , kme , &
1050                                   ips , ipe , jps , jpe , kps , kpe+1 )
1051
1052  CALL stuff_bdy_ijk (vbdy3dtemp1, grid%v_bxs, grid%v_bxe, &
1053                                   grid%v_bys, grid%v_bye, &
1054                                   'N', spec_bdy_width  , &
1055                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1056                                   ims , ime , jms , jme , kms , kme , &
1057                                   ips , ipe , jps , jpe , kps , kpe+1 )
1058
1059  CALL stuff_bdy_ijk (tbdy3dtemp1, grid%t_bxs, grid%t_bxe, &
1060                                   grid%t_bys, grid%t_bye, &
1061                                   'N', spec_bdy_width  , &
1062                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1063                                   ims , ime , jms , jme , kms , kme , &
1064                                   ips , ipe , jps , jpe , kps , kpe+1 )
1065
1066  CALL stuff_bdy_ijk (cwmbdy3dtemp1, grid%cwm_bxs, grid%cwm_bxe, &
1067                                   grid%cwm_bys, grid%cwm_bye, &
1068                                   'N', spec_bdy_width  , &
1069                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1070                                   ims , ime , jms , jme , kms , kme , &
1071                                   ips , ipe , jps , jpe , kps , kpe+1 )
1072
1073  CALL stuff_bdy_ijk (qbdy3dtemp1, grid%q_bxs, grid%q_bxe, &
1074                                   grid%q_bys, grid%q_bye, &
1075                                   'N', spec_bdy_width  , &
1076                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1077                                   ims , ime , jms , jme , kms , kme , &
1078                                   ips , ipe , jps , jpe , kps , kpe+1 )
1079
1080  CALL stuff_bdy_ijk (q2bdy3dtemp1, grid%q2_bxs, grid%q2_bxe, &
1081                                   grid%q2_bys, grid%q2_bye, &
1082                                   'N', spec_bdy_width  , &
1083                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1084                                   ims , ime , jms , jme , kms , kme , &
1085                                   ips , ipe , jps , jpe , kps , kpe+1 )
1086
1087
1088  CALL stuff_bdy_ijk (pdbdy2dtemp1, grid%pd_bxs, grid%pd_bxe, &
1089                                    grid%pd_bys, grid%pd_bye, &
1090                                    'M', spec_bdy_width, &
1091                                    ids , ide+1 , jds , jde+1 , 1 , 1 , &
1092                                    ims , ime , jms , jme , 1 , 1 , &
1093                                    ips , ipe , jps , jpe , 1 , 1 )
1094
1095 !-----------------------------------------------------------------------
1096 !
1097    ELSE IF ( loop .GT. 1 ) THEN
1098 !
1099 !-----------------------------------------------------------------------
1100
1101      call wrf_debug(1,'LOOP>1, so start making non-init boundary conditions')
1102 #if ( HWRF == 1 )
1103
1104      bdytmp_useph: if(grid%use_prep_hybrid) then
1105         call wrf_debug(1,'ALLOCATE PREP_HYBRID BOUNDARY ARRAYS')
1106         !! When running in prep_hybrid mode, we must read in the data here.
1107
1108         ! Allocate boundary arrays:
1109         KB = 2*IDE + JDE - 3
1110         LM = KDE
1111         IM = IDE
1112         JM = JDE
1113         ALLOCATE(TB(KB,LM,2))
1114         ALLOCATE(QB(KB,LM,2))
1115         ALLOCATE(CWMB(KB,LM,2))
1116         ALLOCATE(UB(KB,LM,2))
1117         ALLOCATE(VB(KB,LM,2))
1118         ALLOCATE(PDB(KB,2))
1119         alloc_ph_arrays=.true.
1120
1121         ! Read in the data:
1122         IUNIT_GFS = 900 + LOOP - 1
1123         READ(IUNIT_GFS,iostat=ioerror) PDB,TB,QB,CWMB,UB,VB
1124         if(ioerror/=0) then
1125            write(message,*) 'Unable to read MAKBND output from unit ',IUNIT_GFS
1126            call wrf_error_fatal(message)
1127         endif
1128
1129         ! Now copy the data into the temporary boundary arrays, and
1130         ! switch from IKJ to IJK while we do it.
1131
1132         !!    Southern boundary
1133
1134         IF(JPS.EQ.JDS)THEN
1135            J=1
1136
1137            DO k = kps , MIN(kde,kpe)
1138               N=1
1139               DO i = ips , MIN(ide,ipe)
1140                  tbdy3dtemp2(i,j,k) =   TB(N,k,1)
1141                  qbdy3dtemp2(i,j,k) =   QB(N,k,1)
1142                  cwmbdy3dtemp2(i,j,k) = CWMB(N,k,1)
1143                  q2bdy3dtemp2(i,j,k) =  0.0        !KWON
1144                  write(message,*)'southtend t',k,i,n,tbdy3dtemp2(i,j,k)
1145                  call wrf_debug(10,message)
1146                  write(message,*)'southtend q',k,i,n,qbdy3dtemp2(i,j,k)
1147                  call wrf_debug(10,message)
1148                  if (K .eq. 1 ) then
1149                     write(0,*) 'S boundary values T,Q : ', I,tbdy3dtemp2(i,j,k),  &
1150                          qbdy3dtemp2(i,j,k)
1151                  endif
1152                  N=N+1
1153               END DO
1154            END DO
1155
1156            DO k = kps , MIN(kde,kpe)
1157               N=1
1158               DO i = ips , MIN(ide,ipe)
1159                  ubdy3dtemp2(i,j,k) = UB(N,k,1)
1160                  vbdy3dtemp2(i,j,k) = VB(N,k,1)
1161                  N=N+1
1162               ENDDO
1163            END DO
1164
1165            N=1
1166            DO i = ips , MIN(ide,ipe)
1167               pdbdy2dtemp2(i,j,1) = PDB(N,1)
1168               write(message,*)'southtend p',i,n,pdbdy2dtemp1(i,j,1)
1169               call wrf_debug(10,message)
1170               N=N+1
1171            END DO
1172
1173         ENDIF
1174
1175         !     Northern boundary
1176
1177         IF(JPE.EQ.JDE)THEN
1178
1179            J=MIN(JDE,JPE)
1180            DO k = kps , MIN(kde,kpe)
1181               N=IM+1
1182               DO i = ips , MIN(ide,ipe)
1183                  tbdy3dtemp2(i,j,k) =   TB(N,k,1)
1184                  qbdy3dtemp2(i,j,k) =   QB(N,k,1)
1185                  cwmbdy3dtemp2(i,j,k) = CWMB(N,k,1)
1186                  q2bdy3dtemp2(i,j,k) =  0.0        !KWON
1187                     write(message,*)'northtend t',k,i,n,tbdy3dtemp2(i,j,k)
1188                     call wrf_debug(10,message)
1189                     write(message,*)'northtend q',k,i,n,qbdy3dtemp2(i,j,k)
1190                     call wrf_debug(10,message)
1191                  N=N+1
1192               END DO
1193            END DO
1194
1195            DO k = kps , MIN(kde,kpe)
1196               N=IM
1197               DO i = ips , MIN(ide,ipe)
1198                  ubdy3dtemp2(i,j,k) = UB(N,k,1)
1199                  vbdy3dtemp2(i,j,k) = VB(N,k,1)
1200                  N=N+1
1201               END DO
1202            END DO
1203
1204            N=IM+1
1205            DO i = ips , MIN(ide,ipe)
1206               pdbdy2dtemp2(i,j,1) = PDB(N,1)
1207                  write(message,*)'northtend p',i,n,pdbdy2dtemp1(i,j,1)
1208                  call wrf_debug(10,message)
1209               N=N+1
1210            END DO
1211
1212         ENDIF
1213
1214         !!     Western boundary
1215
1216         IF(IPS.EQ.IDS)THEN
1217            I=1
1218            DO k = kps , MIN(kde,kpe)
1219               N=2*IM+1
1220               inc_h=mod(jps+1,2)
1221               DO j = jps+inc_h, MIN(jde,jpe),2
1222                  if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1223                     tbdy3dtemp2(i,j,k) =   TB(N,k,1)
1224                     qbdy3dtemp2(i,j,k) =   QB(N,k,1)
1225                     cwmbdy3dtemp2(i,j,k) = CWMB(N,k,1)
1226                     q2bdy3dtemp2(i,j,k) =  0.0        !KWON
1227                        write(message,*)'westtend t',k,j,n,tbdy3dtemp2(i,j,k)
1228                        call wrf_debug(10,message)
1229                        write(message,*)'westtend q',k,j,n,qbdy3dtemp2(i,j,k)
1230                        call wrf_debug(10,message)
1231                     N=N+1
1232                  endif
1233               END DO
1234            END DO
1235
1236            DO k = kps , MIN(kde,kpe)
1237               N=2*IM-1
1238               inc_v=mod(jps,2)
1239               DO j = jps+inc_v, MIN(jde,jpe),2
1240                  if (J .ge. 2 .and. J .le. jde-1 .and. mod(J,2) .eq. 0) then
1241                     ubdy3dtemp2(i,j,k) = UB(N,k,1)
1242                     vbdy3dtemp2(i,j,k) = VB(N,k,1)
1243                     N=N+1
1244                  endif
1245               END DO
1246            END DO
1247
1248            N=2*IM+1
1249            inc_h=mod(jps+1,2)
1250            DO j = jps+inc_h, MIN(jde,jpe),2
1251               if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1252                  pdbdy2dtemp2(i,j,1) = PDB(N,1)
1253                     write(message,*)'westtend p',j,n,pdbdy2dtemp1(i,j,1)
1254                     call wrf_debug(10,message)
1255                  N=N+1
1256               endif
1257            END DO
1258
1259         ENDIF
1260
1261         !!     Eastern boundary
1262
1263         IF(IPE.EQ.IDE)THEN
1264
1265            I=MIN(IDE,IPE)
1266
1267            DO k = kps , MIN(kde,kpe)
1268               N=2*IM+(JM/2)
1269               inc_h=mod(jps+1,2)
1270               DO j = jps+inc_h, MIN(jde,jpe),2
1271                  if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1272                     tbdy3dtemp2(i,j,k) =   TB(N,k,1)
1273                     qbdy3dtemp2(i,j,k) =   QB(N,k,1)
1274                     cwmbdy3dtemp2(i,j,k) = CWMB(N,k,1)
1275                     q2bdy3dtemp2(i,j,k) =  0.0        !KWON
1276                        write(message,*)'easttend t',k,j,n,tbdy3dtemp2(i,j,k)
1277                        call wrf_debug(10,message)
1278                        write(message,*)'easttend q',k,j,n,qbdy3dtemp2(i,j,k)
1279                        call wrf_debug(10,message)
1280                     N=N+1
1281                  endif
1282               END DO
1283            END DO
1284
1285            DO k = kps , MIN(kde,kpe)
1286               N=2*IM+(JM/2)-1
1287               inc_v=mod(jps,2)
1288               DO j = jps+inc_v, MIN(jde,jpe),2
1289                  if (J .ge. 2 .and. J .le. jde-1 .and. mod(J,2) .eq. 0) then
1290                     ubdy3dtemp2(i,j,k) = UB(N,k,1)
1291                     vbdy3dtemp2(i,j,k) = VB(N,k,1)
1292                     N=N+1
1293                  endif
1294               END DO
1295            END DO
1296
1297            N=2*IM+(JM/2)
1298            inc_h=mod(jps+1,2)
1299            DO j = jps+inc_h, MIN(jde,jpe),2
1300               if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1301                  pdbdy2dtemp2(i,j,1) = PDB(N,1)
1302                     write(message,*)'easttend p',j,n,pdbdy2dtemp1(i,j,1)
1303                     call wrf_debug(10,message)
1304                  N=N+1
1305               endif
1306            END DO
1307
1308         ENDIF
1309      else
1310 #endif
1311
1312            CALL output_auxinput4 ( id4, grid , config_flags , ierr )
1313
1314 #if ( HWRF == 1 )
1315      endif bdytmp_useph
1316 #endif
1317
1318       write(message,*)' assemble_output loop=',loop,' in IF block'
1319       call wrf_message(message)
1320
1321       !  Open the boundary file.
1322
1323       IF ( loop .eq. 2 ) THEN
1324          CALL construct_filename1( bdyname , 'wrfbdy' , grid%id , 2 )
1325       CALL open_w_dataset ( id, TRIM(bdyname) , grid , config_flags , &
1326                           output_boundary , "DATASET=BOUNDARY", ierr )
1327          IF ( ierr .NE. 0 ) THEN
1328                CALL wrf_error_fatal( 'ideal_hwrf: error opening wrfbdy for writing' )
1329          ENDIF
1330 !         grid%write_metadata = .true.
1331       ELSE
1332 ! what's this do?
1333 !         grid%write_metadata = .true.
1334 !         grid%write_metadata = .false.
1335          CALL domain_clockadvance( grid )
1336       END IF
1337
1338 #if ( HWRF == 1 )
1339      bdytmp_notph: if(.not.grid%use_prep_hybrid) then
1340 #endif
1341 !-----------------------------------------------------------------------
1342 !***  SOUTHERN BOUNDARY
1343 !-----------------------------------------------------------------------
1344 !
1345         IF(JPS==JDS)THEN
1346           J=1
1347           DO k = kps , MIN(kde,kpe)
1348           DO i = ips , MIN(ide,ipe)
1349             ubdy3dtemp2(i,j,k) = grid%u(i,j,k)
1350             vbdy3dtemp2(i,j,k) = grid%v(i,j,k)
1351             tbdy3dtemp2(i,j,k) = grid%t(i,j,k)
1352             qbdy3dtemp2(i,j,k) = grid%q(i,j,k)
1353             cwmbdy3dtemp2(i,j,k) = grid%cwm(i,j,k)
1354             q2bdy3dtemp2(i,j,k) = grid%q2(i,j,k)
1355           END DO
1356           END DO
1357 !
1358           DO i = ips , MIN(ide,ipe)
1359             pdbdy2dtemp2(i,j,1) = grid%pd(i,j)
1360           END DO
1361         ENDIF
1362
1363 !
1364 !-----------------------------------------------------------------------
1365 !***  NORTHERN BOUNDARY
1366 !-----------------------------------------------------------------------
1367 !
1368         IF(JPE==JDE)THEN
1369           J=MIN(JDE,JPE)
1370           DO k = kps , MIN(kde,kpe)
1371           DO i = ips , MIN(ide,ipe)
1372             ubdy3dtemp2(i,j,k) = grid%u(i,j,k)
1373             vbdy3dtemp2(i,j,k) = grid%v(i,j,k)
1374             tbdy3dtemp2(i,j,k) = grid%t(i,j,k)
1375             qbdy3dtemp2(i,j,k) = grid%q(i,j,k)
1376             cwmbdy3dtemp2(i,j,k) = grid%cwm(i,j,k)
1377             q2bdy3dtemp2(i,j,k) = grid%q2(i,j,k)
1378           END DO
1379           END DO
1380
1381           DO i = ips , MIN(ide,ipe)
1382             pdbdy2dtemp2(i,j,1) = grid%pd(i,j)
1383           END DO
1384         ENDIF
1385 !
1386 !-----------------------------------------------------------------------
1387 !***  WESTERN BOUNDARY
1388 !-----------------------------------------------------------------------
1389 !
1390         IF(IPS==IDS)THEN
1391           I=1
1392           DO k = kps , MIN(kde,kpe)
1393           inc_h=mod(jps+1,2)
1394       if(k==1)then
1395         write(message,*)' assemble_ouput loop=',loop,' inc_h=',inc_h,' jps=',jps
1396         call wrf_debug(10,message)
1397       endif
1398           DO j = jps+inc_h, MIN(jde,jpe),2
1399         if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1400             tbdy3dtemp2(i,j,k) = grid%t(i,j,k)
1401       if(k==1)then
1402         write(message,*)' loop=',loop,' i=',i,' j=',j,' tbdy3dtemp1(i,j,k)=',tbdy3dtemp1(i,j,k)
1403         call wrf_debug(10,message)
1404       endif
1405             qbdy3dtemp2(i,j,k) = grid%q(i,j,k)
1406             cwmbdy3dtemp2(i,j,k) = grid%cwm(i,j,k)
1407             q2bdy3dtemp2(i,j,k) = grid%q2(i,j,k)
1408         endif
1409           END DO
1410           END DO
1411 !
1412           DO k = kps , MIN(kde,kpe)
1413           inc_v=mod(jps,2)
1414           DO j = jps+inc_v, MIN(jde,jpe),2
1415         if (J .ge. 2 .and. J .le. jde-1 .and. mod(J,2) .eq. 0) then
1416             ubdy3dtemp2(i,j,k) = grid%u(i,j,k)
1417             vbdy3dtemp2(i,j,k) = grid%v(i,j,k)
1418         endif
1419           END DO
1420           END DO
1421
1422           inc_h=mod(jps+1,2)
1423         DO j = jps+inc_h, MIN(jde,jpe),2
1424         if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1425           pdbdy2dtemp2(i,j,1) = grid%pd(i,j)
1426           write(message,*)' loop=',loop,' i=',i,' j=',j,' pdbdy2dtemp1(i,j)=',pdbdy2dtemp1(i,j,1)
1427           CALL wrf_debug(10,message)
1428         endif
1429           END DO
1430         ENDIF
1431 !
1432 !-----------------------------------------------------------------------
1433 !***  EASTERN BOUNDARY
1434 !-----------------------------------------------------------------------
1435 !
1436         IF(IPE==IDE)THEN
1437           I=MIN(IDE,IPE)
1438
1439           DO k = kps , MIN(kde,kpe)
1440           inc_h=mod(jps+1,2)
1441           DO j = jps+inc_h, MIN(jde,jpe),2
1442         if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1443             tbdy3dtemp2(i,j,k) = grid%t(i,j,k)
1444             qbdy3dtemp2(i,j,k) = grid%q(i,j,k)
1445             cwmbdy3dtemp2(i,j,k) = grid%cwm(i,j,k)
1446             q2bdy3dtemp2(i,j,k) = grid%q2(i,j,k)
1447         endif
1448           END DO
1449           END DO
1450
1451           DO k = kps , MIN(kde,kpe)
1452           inc_v=mod(jps,2)
1453           DO j = jps+inc_v, MIN(jde,jpe),2
1454         if (J .ge. 2 .and. J .le. jde-1 .and. mod(J,2) .eq. 0) then
1455             ubdy3dtemp2(i,j,k) = grid%u(i,j,k)
1456             vbdy3dtemp2(i,j,k) = grid%v(i,j,k)
1457         endif
1458           END DO
1459           END DO
1460
1461           inc_h=mod(jps+1,2)
1462           DO j = jps+inc_h, MIN(jde,jpe),2
1463         if (J .ge. 3 .and. J .le. jde-2 .and. mod(J,2) .eq. 1) then
1464             pdbdy2dtemp2(i,j,1) = grid%pd(i,j)
1465         endif
1466           END DO
1467         ENDIF
1468 #if ( HWRF == 1 )
1469      endif bdytmp_notph
1470 #endif
1471 !-----------------------------------------------------------------------
1472       !  During all of the loops after the first loop,
1473       !  we first compute the boundary
1474       !  tendencies with the current data values
1475       !  (*bdy3dtemp2 arrays) and the previously
1476       !  saved information stored in the *bdy3dtemp1 arrays.
1477
1478
1479       CALL stuff_bdytend_ijk ( ubdy3dtemp2 , ubdy3dtemp1 , REAL(interval_seconds),&
1480                                    grid%u_btxs, grid%u_btxe, &
1481                                    grid%u_btys, grid%u_btye, &
1482                                    'N',  spec_bdy_width      , &
1483                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1484                                    ims , ime , jms , jme , kms , kme , &
1485                                    ips , ipe , jps , jpe , kps , kpe+1 )
1486
1487       CALL stuff_bdytend_ijk ( vbdy3dtemp2 , vbdy3dtemp1 , REAL(interval_seconds),&
1488                                    grid%v_btxs, grid%v_btxe, &
1489                                    grid%v_btys, grid%v_btye, &
1490                                    'N',  spec_bdy_width      , &
1491                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1492                                    ims , ime , jms , jme , kms , kme , &
1493                                    ips , ipe , jps , jpe , kps , kpe+1 )
1494
1495       CALL stuff_bdytend_ijk ( tbdy3dtemp2 , tbdy3dtemp1 , REAL(interval_seconds),&
1496                                    grid%t_btxs, grid%t_btxe, &
1497                                    grid%t_btys, grid%t_btye, &
1498                                    'N',  spec_bdy_width      , &
1499                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1500                                    ims , ime , jms , jme , kms , kme , &
1501                                    ips , ipe , jps , jpe , kps , kpe+1 )
1502
1503       CALL stuff_bdytend_ijk ( cwmbdy3dtemp2 , cwmbdy3dtemp1 , REAL(interval_seconds),&
1504                                    grid%cwm_btxs, grid%cwm_btxe, &
1505                                    grid%cwm_btys, grid%cwm_btye, &
1506                                    'N',  spec_bdy_width      , &
1507                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1508                                    ims , ime , jms , jme , kms , kme , &
1509                                    ips , ipe , jps , jpe , kps , kpe+1 )
1510
1511       CALL stuff_bdytend_ijk ( qbdy3dtemp2 , qbdy3dtemp1 , REAL(interval_seconds),&
1512                                    grid%q_btxs, grid%q_btxe, &
1513                                    grid%q_btys, grid%q_btye, &
1514                                    'N',  spec_bdy_width      , &
1515                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1516                                    ims , ime , jms , jme , kms , kme , &
1517                                    ips , ipe , jps , jpe , kps , kpe+1 )
1518
1519       CALL stuff_bdytend_ijk ( q2bdy3dtemp2 , q2bdy3dtemp1 , REAL(interval_seconds),&
1520                                    grid%q2_btxs, grid%q2_btxe, &
1521                                    grid%q2_btys, grid%q2_btye, &
1522                                    'N',  spec_bdy_width      , &
1523                                    ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1524                                    ims , ime , jms , jme , kms , kme , &
1525                                    ips , ipe , jps , jpe , kps , kpe+1 )
1526
1527       CALL stuff_bdytend_ijk( pdbdy2dtemp2 , pdbdy2dtemp1, REAL(interval_seconds),&
1528                                    grid%pd_btxs, grid%pd_btxe, &
1529                                    grid%pd_btys, grid%pd_btye, &
1530                                    'M',  spec_bdy_width      , &
1531                                    ids , ide+1 , jds , jde+1 , 1 , 1 , &
1532                                    ims , ime   , jms , jme   , 1 , 1 , &
1533                                    ips , ipe   , jps , jpe   , 1 , 1 )
1534
1535
1536
1537       !  Both pieces of the boundary data are now
1538       !  available to be written (initial time and tendency).
1539       !  This looks ugly, these date shifting things.
1540       !  What's it for?  We want the "Times" variable
1541       !  in the lateral BDY file to have the valid times
1542       !  of when the initial fields are written.
1543       !  That's what the loop-2 thingy is for with the start date.
1544       !  We increment the start_date so
1545       !  that the starting time in the attributes is the
1546       !  second time period.  Why you may ask.  I
1547       !  agree, why indeed.
1548
1549       temp24= current_date
1550       temp24b=start_date
1551       start_date = current_date
1552       CALL geth_newdate ( temp19 , temp24b(1:19) , &
1553                          (loop-2) * model_config_rec%interval_seconds )
1554       current_date = temp19 //  '.0000'
1555        CALL domain_clock_set( grid, current_date(1:19) )
1556       write(message,*) 'LBC valid between these times ',current_date, ' ',start_date
1557       CALL wrf_message(message)
1558
1559       CALL output_boundary ( id, grid , config_flags , ierr )
1560       current_date = temp24
1561       start_date = temp24b
1562
1563       !  OK, for all of the loops, we output the initialzation
1564       !  data, which would allow us to
1565       !  start the model at any of the available analysis time periods.
1566
1567 !  WRITE ( loop_char , FMT = '(I4.4)' ) loop
1568 !  CALL open_w_dataset ( id1, 'wrfinput'//loop_char , grid , config_flags , output_input , "DATASET=INPUT", ierr )
1569 !  IF ( ierr .NE. 0 ) THEN
1570 !    CALL wrf_error_fatal( 'ideal_hwrf: error opening wrfinput'//loop_char//' for writing' )
1571 !  ENDIF
1572 !  grid%write_metadata = .true.
1573
1574 !  CALL calc_current_date ( grid%id , 0. )
1575 !  CALL output_input ( id1, grid , config_flags , ierr )
1576 !  CALL close_dataset ( id1 , config_flags , "DATASET=INPUT" )
1577
1578   !  Is this or is this not the last time time?  We can remove some unnecessary
1579   !  stores if it is not.
1580
1581       IF     ( loop .LT. time_loop_max ) THEN
1582
1583          !  We need to save the 3d data to compute a
1584          !  difference during the next loop.  Couple the
1585          !  3d fields with total mu (mub + mu_2) and the
1586          !  stagger-specific map scale factor.
1587          !  We load up the boundary data again for use in the next loop.
1588
1589
1590 !mp     change these limits?????????
1591
1592          DO k = kps , kpe
1593             DO j = jps , jpe
1594                DO i = ips , ipe
1595                   ubdy3dtemp1(i,j,k) = ubdy3dtemp2(i,j,k)
1596                   vbdy3dtemp1(i,j,k) = vbdy3dtemp2(i,j,k)
1597                   tbdy3dtemp1(i,j,k) = tbdy3dtemp2(i,j,k)
1598                   cwmbdy3dtemp1(i,j,k) = cwmbdy3dtemp2(i,j,k)
1599                   qbdy3dtemp1(i,j,k) = qbdy3dtemp2(i,j,k)
1600                   q2bdy3dtemp1(i,j,k) = q2bdy3dtemp2(i,j,k)
1601                END DO
1602             END DO
1603          END DO
1604
1605 !mp     change these limits?????????
1606
1607          DO j = jps , jpe
1608             DO i = ips , ipe
1609                pdbdy2dtemp1(i,j,1) = pdbdy2dtemp2(i,j,1)
1610             END DO
1611          END DO
1612
1613   !  There are 2 components to the lateral boundaries.
1614   !   First, there is the starting
1615   !  point of this time period - just the outer few rows and columns.
1616
1617  CALL stuff_bdy_ijk (ubdy3dtemp1, grid%u_bxs, grid%u_bxe, &
1618                                   grid%u_bys, grid%u_bye, &
1619                                   'N', spec_bdy_width  , &
1620                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1621                                   ims , ime , jms , jme , kms , kme , &
1622                                   ips , ipe , jps , jpe , kps , kpe+1 )
1623
1624  CALL stuff_bdy_ijk (vbdy3dtemp1, grid%v_bxs, grid%v_bxe, &
1625                                   grid%v_bys, grid%v_bye, &
1626                                   'N', spec_bdy_width  , &
1627                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1628                                   ims , ime , jms , jme , kms , kme , &
1629                                   ips , ipe , jps , jpe , kps , kpe+1 )
1630
1631  CALL stuff_bdy_ijk (tbdy3dtemp1, grid%t_bxs, grid%t_bxe, &
1632                                   grid%t_bys, grid%t_bye, &
1633                                   'N', spec_bdy_width  , &
1634                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1635                                   ims , ime , jms , jme , kms , kme , &
1636                                   ips , ipe , jps , jpe , kps , kpe+1 )
1637
1638  CALL stuff_bdy_ijk (cwmbdy3dtemp1, grid%cwm_bxs, grid%cwm_bxe, &
1639                                   grid%cwm_bys, grid%cwm_bye, &
1640                                   'N', spec_bdy_width  , &
1641                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1642                                   ims , ime , jms , jme , kms , kme , &
1643                                   ips , ipe , jps , jpe , kps , kpe+1 )
1644
1645  CALL stuff_bdy_ijk (qbdy3dtemp1, grid%q_bxs, grid%q_bxe, &
1646                                   grid%q_bys, grid%q_bye, &
1647                                   'N', spec_bdy_width  , &
1648                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1649                                   ims , ime , jms , jme , kms , kme , &
1650                                   ips , ipe , jps , jpe , kps , kpe+1 )
1651
1652  CALL stuff_bdy_ijk (q2bdy3dtemp1, grid%q2_bxs, grid%q2_bxe, &
1653                                   grid%q2_bys, grid%q2_bye, &
1654                                   'N', spec_bdy_width  , &
1655                                   ids , ide+1 , jds , jde+1 , kds , kde+1 , &
1656                                   ims , ime , jms , jme , kms , kme , &
1657                                   ips , ipe , jps , jpe , kps , kpe+1 )
1658
1659  CALL stuff_bdy_ijk (pdbdy2dtemp1,grid%pd_bxs, grid%pd_bxe, &
1660                                   grid%pd_bys, grid%pd_bye, &
1661                                   'M', spec_bdy_width  , &
1662                                   ids , ide+1 , jds , jde+1 , 1 , 1 , &
1663                                   ims , ime , jms , jme , 1 , 1 , &
1664                                   ips , ipe , jps , jpe , 1 , 1 )
1665
1666       ELSE IF ( loop .EQ. time_loop_max ) THEN
1667
1668     !  If this is the last time through here, we need to close the files.
1669
1670          CALL close_dataset ( id , config_flags , "DATASET=BOUNDARY" )
1671
1672       END IF
1673
1674    END IF main_loop_test
1675
1676 #if ( HWRF == 1 )
1677   if(alloc_ph_arrays) then
1678      call wrf_debug(1,'DEALLOCATE PREP_HYBRID BOUNARY ARRAYS')
1679      deallocate(TB,QB,CWMB,UB,VB,PDB)
1680   endif
1681 #endif
1682
1683 END SUBROUTINE assemble_output