main/real_nmm.F

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