dyn_em/start_em.F

   1 !-------------------------------------------------------------------
   2
   3    SUBROUTINE start_domain_em ( grid, allowed_to_read &
   4 ! Actual arguments generated from Registry
   5 # include "dummy_new_args.inc"
   6 !
   7 )
   8
   9    USE module_domain, ONLY : domain, wrfu_timeinterval, get_ijk_from_grid, &
  10         domain_setgmtetc
  11    USE module_state_description
  12    USE module_driver_constants
  13    USE module_wrf_error
  14    USE module_model_constants
  15    USE module_bc, ONLY : boundary_condition_check, set_physical_bc2d, set_physical_bc3d, bdyzone
  16    USE module_bc_em, ONLY: lbc_fcx_gcx, set_w_surface
  17    USE module_configure, ONLY : model_to_grid_config_rec, model_config_rec, grid_config_rec_type
  18    USE module_tiles, ONLY : set_tiles
  19 #ifdef DM_PARALLEL
  20    USE module_dm, ONLY : wrf_dm_min_real, wrf_dm_max_real, wrf_dm_maxval, &
  21         ntasks_x, ntasks_y, &
  22         local_communicator_periodic, local_communicator, mytask, ntasks
  23 #else
  24    USE module_dm, ONLY : wrf_dm_min_real, wrf_dm_max_real
  25 #endif
  26    USE module_comm_dm
  27    USE module_llxy, ONLY : proj_cassini
  28    USE module_physics_init
  29    USE module_sf_noahmpdrv, ONLY : groundwater_init
  30    USE module_lightning_driver, ONLY : lightning_init
  31    USE module_fr_fire_driver_wrf, ONLY : fire_driver_em_init
  32    USE module_firebrand_spotting, ONLY : firebrand_spotting_em_init
  33    USE module_stoch, ONLY : setup_rand_perturb, rand_seed, update_stoch, initialize_stoch
  34    USE module_trajectory, ONLY : trajectory_init
  35 #if (WRF_CHEM == 1)
  36    USE module_aerosols_sorgam, ONLY: sum_pm_sorgam
  37    USE module_gocart_aerosols, ONLY: sum_pm_gocart
  38    USE module_mosaic_driver, ONLY: sum_pm_mosaic
  39    USE module_input_tracer, ONLY: initialize_tracer
  40    USE module_aerosols_soa_vbs, only: sum_pm_soa_vbs
  41 #if (WRF_USE_CLM == 1)
  42    USE shr_megan_mod, only : wrf_spc_chk
  43 #endif
  44 #endif
  45    USE module_diag_pld, ONLY : pld
  46    USE module_diag_zld, ONLY : zld
  47    USE module_trad_fields, ONLY : trad_fields
  48
  49 !!debug
  50 !USE module_compute_geop
  51
  52    USE module_model_constants
  53    USE module_avgflx_em, ONLY : zero_avgflx
  54
  55    IMPLICIT NONE
  56    !  Input data.
  57    TYPE (domain)          :: grid
  58
  59    LOGICAL , INTENT(IN)   :: allowed_to_read
  60
  61    !  Definitions of dummy arguments to this routine (generated from Registry).
  62 # include "dummy_new_decl.inc"
  63
  64    !  Structure that contains run-time configuration (namelist) data for domain
  65    TYPE (grid_config_rec_type)              :: config_flags
  66
  67    !  Local data
  68    INTEGER                             ::                       &
  69                                   ids, ide, jds, jde, kds, kde, &
  70                                   ims, ime, jms, jme, kms, kme, &
  71                                   ips, ipe, jps, jpe, kps, kpe, &
  72                                   its, ite, jts, jte, kts, kte, &
  73                                   ij,i,j,k,ii,jj,kk,loop,error,l
  74
  75    INTEGER ::              imsx, imex, jmsx, jmex, kmsx, kmex,    &
  76                            ipsx, ipex, jpsx, jpex, kpsx, kpex,    &
  77                            imsy, imey, jmsy, jmey, kmsy, kmey,    &
  78                            ipsy, ipey, jpsy, jpey, kpsy, kpey
  79
  80    INTEGER     :: i_m
  81
  82    REAL        :: p_top_test, p00, t00, a, tiso, p_surf, pd_surf, temp, tiso_tmp
  83    REAL        :: p_strat, a_strat
  84 #if (WRF_CHEM == 1)
  85    REAL        RGASUNIV ! universal gas constant [ J/mol-K ]
  86    PARAMETER ( RGASUNIV = 8.314510 )
  87       REAL,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: &
  88                       z_at_w,convfac
  89    REAL :: tempfac
  90    INTEGER     :: numgas
  91 #endif
  92
  93    REAL :: qvf1, qvf2, qvf
  94    REAL :: pfu, pfd, phm
  95    REAL :: MPDT
  96    REAL :: spongeweight
  97    LOGICAL :: first_trip_for_this_domain, start_of_simulation, fill_w_flag
  98    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
  99
 100 #if (WRF_CHEM != 1)
 101       REAL,ALLOCATABLE,DIMENSION(:,:,:) :: cldfra_old
 102 #endif
 103
 104    REAL :: lat1 , lat2 , lat3 , lat4
 105    REAL :: lon1 , lon2 , lon3 , lon4
 106    INTEGER :: num_points_lat_lon , iloc , jloc
 107    CHARACTER (LEN=256) :: message, a_message
 108    TYPE(WRFU_TimeInterval) :: stepTime
 109    REAL, DIMENSION(:,:), ALLOCATABLE :: clat_glob
 110    logical :: f_flux  ! flag for computing averaged fluxes in cu_gd
 111
 112    INTEGER :: idex, jdex
 113    INTEGER :: im1,ip1,jm1,jp1
 114    REAL :: hx,hy,pi
 115
 116    REAL :: w_max, w_min
 117    LOGICAL :: w_needs_to_be_set
 118 !
 119 ! Define variables local (topo_wind local vars)
 120    REAL    :: alpha, vfac
 121    INTEGER :: j_save
 122
 123    !  For a global domain
 124
 125    INTEGER :: alloc_status
 126    CHARACTER (LEN=256) :: alloc_err_message
 127
 128 !..Need to fill special height var for setting up initial condition.  G. Thompson
 129    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: z_at_q
 130    REAL, ALLOCATABLE, DIMENSION(:,:,:) :: dz8w   !mchen
 131
 132 !  CCN for MP=18 initializatio
 133    REAL :: ccn_max_val
 134
 135    REAL :: dt_s, dx_km
 136
 137    REAL :: max_mf, max_rot_angle
 138    CALL get_ijk_from_grid ( grid ,                              &
 139                            ids, ide, jds, jde, kds, kde,        &
 140                            ims, ime, jms, jme, kms, kme,        &
 141                            ips, ipe, jps, jpe, kps, kpe,        &
 142                            imsx, imex, jmsx, jmex, kmsx, kmex,  &
 143                            ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
 144                            imsy, imey, jmsy, jmey, kmsy, kmey,  &
 145                            ipsy, ipey, jpsy, jpey, kpsy, kpey )
 146
 147    kts = kps ; kte = kpe     ! note that tile is entire patch
 148    its = ips ; ite = ipe     ! note that tile is entire patch
 149    jts = jps ; jte = jpe    ! note that tile is entire patch
 150 #if (WRF_CHEM != 1)
 151          ALLOCATE(CLDFRA_OLD(IMS:IME,KMS:KME,JMS:JME),STAT=I)  ; CLDFRA_OLD = 0.
 152 #endif
 153    ALLOCATE(dz8w(IMS:IME,KMS:KME,JMS:JME),STAT=I)  ; dz8w = 0.
 154    ALLOCATE(z_at_q(IMS:IME,KMS:KME,JMS:JME),STAT=I)  ; z_at_q = 0.
 155    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 156
 157    IF ( ( MOD (ide-ids,config_flags%parent_grid_ratio) .NE. 0 ) .OR. &
 158         ( MOD (jde-jds,config_flags%parent_grid_ratio) .NE. 0 ) ) THEN
 159       WRITE(message, FMT='(A,I2,":  Both MOD(",I4,"-",I1,",",I2,") and MOD(",I4,"-",I1,",",I2,") must = 0" )') &
 160          "Nested dimensions are illegal for domain ",grid%id,ide,ids,config_flags%parent_grid_ratio,&
 161          jde,jds,config_flags%parent_grid_ratio
 162       CALL wrf_error_fatal ( message )
 163    END IF
 164
 165    IF ( config_flags%polar ) THEN
 166 !write(0,*)"<stdin>",__LINE__,' clat ',ips,ipe,jps,jpe
 167 !do j = jps,jpe
 168 !write(0,*)"<stdin>",__LINE__,' clat ',ids,j,grid%clat(ips,j)
 169 !enddo
 170
 171 #ifdef DM_PARALLEL
 172 ! WARNING:  this might present scaling issues on very large numbers of processors
 173      alloc_err_message = ' '
 174      alloc_err_message(1:12) = 'NO PROBLEMOS'
 175 #if 0
 176      ALLOCATE( clat_glob(ids:ide,jds:jde), STAT=alloc_status, ERRMSG=alloc_err_message )
 177 #else
 178      ALLOCATE( clat_glob(ids:ide,jds:jde), STAT=alloc_status)
 179      alloc_err_message = 'Allocation of space for a global field failed.'
 180 #endif
 181
 182      IF ( alloc_status .NE. 0 ) THEN
 183         CALL wrf_message ( TRIM(alloc_err_message) )
 184         CALL wrf_error_fatal ( 'Error allocating entire domain size of 2d array CLAT for global domain' )
 185      END IF
 186
 187      CALL wrf_patch_to_global_real ( grid%clat, clat_glob, grid%domdesc, 'xy', 'xy', &
 188                                      ids, ide, jds, jde, 1, 1, &
 189                                      ims, ime, jms, jme, 1, 1, &
 190                                      its, ite, jts, jte, 1, 1 )
 191
 192      CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
 193
 194      grid%clat_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
 195
 196      find_j_index_of_fft_filter : DO j = jds , jde-1
 197         IF ( ABS(clat_glob(ids,j)) .LE. config_flags%fft_filter_lat ) THEN
 198            j_save = j
 199            EXIT find_j_index_of_fft_filter
 200         END IF
 201      END DO find_j_index_of_fft_filter
 202
 203      CALL wrf_patch_to_global_real ( grid%msft, clat_glob, grid%domdesc, 'xy', 'xy', &
 204                                      ids, ide, jds, jde, 1, 1, &
 205                                      ims, ime, jms, jme, 1, 1, &
 206                                      its, ite, jts, jte, 1, 1 )
 207
 208      CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
 209
 210      grid%mf_fft = clat_glob(ids,j_save)
 211
 212      grid%mf_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
 213
 214      DEALLOCATE( clat_glob )
 215 #endif
 216    ENDIF
 217
 218 ! here we check to see if the boundary conditions are set properly
 219
 220    CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
 221 ! make sure that topo_wind option has var_sso data available
 222     IF ((config_flags%topo_wind .EQ. 1) .AND. (.NOT. grid%got_var_sso)) THEN
 223       CALL wrf_error_fatal ("topo_wind requires VAR_SSO data")
 224     ENDIF
 225
 226 !kludge - need to stop CG from resetting precip and phys tendencies to zero
 227 !         when we are in here due to a nest being spawned, we want to still
 228 !         recompute the base state, but that is about it
 229    !  This is temporary and will need to be changed when grid%itimestep is removed.
 230
 231    IF ( grid%itimestep .EQ. 0 ) THEN
 232       first_trip_for_this_domain = .TRUE.
 233    ELSE
 234       first_trip_for_this_domain = .FALSE.
 235    END IF
 236
 237    IF ( .not. ( config_flags%restart .or. grid%moved ) ) THEN
 238        grid%itimestep=0
 239    ENDIF
 240
 241    IF ( config_flags%restart .or. grid%moved .or. config_flags%cycling) THEN
 242        first_trip_for_this_domain = .TRUE.
 243    ENDIF
 244
 245    ! Print out the maximum map scale factor on the coarse domain
 246
 247    IF ( ( first_trip_for_this_domain ) .AND. ( grid%id .EQ. 1 ) .AND. &
 248         ( .NOT. config_flags%polar ) ) THEN
 249       max_mf = grid%msft(its,jts)
 250       DO j=jts,MIN(jde-1,jte)
 251          DO i=its,MIN(ide-1,ite)
 252             max_mf = MAX ( max_mf , grid%msft(i,j) )
 253          END DO
 254       END DO
 255 #if ( defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
 256       max_mf = wrf_dm_max_real ( max_mf )
 257 #endif
 258       WRITE ( a_message , FMT='(A,F5.2,A)' ) 'Max map factor in domain 1 = ',max_mf, &
 259                                              '. Scale the dt in the model accordingly.'
 260       CALL wrf_message ( a_message )
 261    END IF
 262
 263 !-----------------------------------------------------------------------
 264 ! Do a gross check on the time step for the most coarse grid (not for
 265 ! adaptive time steps. If the value is significantly larger than reasonable,
 266 ! then stop immediately. Likely the user made a mistake.
 267 ! If the rule of thumb is 6 DX = DT, then anything > 10 is flagged.
 268 ! Only look at domain 1, only do this for the firt time into this routine,
 269 ! only look at real-data cases (where the map projection is not Cartesian),
 270 ! and only look at user-specified time steps (not the adaptive dt option).
 271 !-----------------------------------------------------------------------
 272    IF ( ( grid%id .EQ. 1 ) .AND. &
 273         ( first_trip_for_this_domain ) .AND. &
 274         ( config_flags%map_proj .NE. 0 ) .AND. &
 275         ( .NOT. config_flags%use_adaptive_time_step ) ) THEN
 276       dt_s = REAL(config_flags%time_step) + &
 277              REAL(config_flags%time_step_fract_num) / &
 278              REAL(config_flags%time_step_fract_den)
 279       dx_km = MIN ( config_flags%dx , config_flags%dy ) / 1000.
 280       WRITE (message,*) 'D01: Time step                              = ',dt_s ,' (s)'
 281       CALL wrf_message ( TRIM(message) )
 282       WRITE (message,*) 'D01: Grid Distance                          = ',dx_km ,' (km)'
 283       CALL wrf_message ( TRIM(message) )
 284       WRITE (message,*) 'D01: Grid Distance Ratio dt/dx              = ', dt_s / dx_km, ' (s/km)'
 285       CALL wrf_message ( TRIM(message) )
 286       WRITE (message,*) 'D01: Ratio Including Maximum Map Factor     = ', dt_s / (dx_km / max_mf) , ' (s/km)'
 287       CALL wrf_message ( TRIM(message) )
 288       WRITE (message,*) 'D01: NML defined reasonable_time_step_ratio = ', config_flags%reasonable_time_step_ratio
 289       CALL wrf_message ( TRIM(message) )
 290       IF ( dt_s / dx_km > config_flags%reasonable_time_step_ratio ) THEN
 291          CALL wrf_message ( 'The time step is probably too large for this grid distance, reduce it.' )
 292          WRITE (message, * ) 'If you are sure of your settings, set reasonable_time_step_ratio in namelist.input > ' &
 293                              ,dt_s / (dx_km / max_mf)
 294          CALL wrf_message ( TRIM(message) )
 295          CALL wrf_error_fatal ( '--- ERROR: Time step too large')
 296       END IF
 297    END IF
 298
 299     if(config_flags%cycling) then
 300 ! Clear the buckets for diagnostics at initial time
 301        DO j = jts,min(jte,jde-1)
 302        DO i = its, min(ite,ide-1)
 303             grid%prec_acc_nc(i,j) = 0.
 304             grid%snow_acc_nc(i,j)  = 0.
 305             grid%sfcrunoff  (i,j) = 0.
 306             grid%udrunoff   (i,j) = 0.
 307 ! acsnow, and acrunoff are run-total
 308             grid%acrunoff   (i,j) = 0.
 309             grid%acsnow     (i,j) = 0.
 310        ENDDO
 311        ENDDO
 312     endif
 313
 314 !   --- SETUP AND INITIALIZE STOCHASTIC PERTURBATION SCHEMES ---
 315
 316    IF ( grid%itimestep .EQ. 0 ) THEN
 317       first_trip_for_this_domain = .TRUE.
 318    ELSE
 319       first_trip_for_this_domain = .FALSE.
 320    END IF
 321
 322    IF ( .not. ( config_flags%restart .or. grid%moved .or. config_flags%hrrr_cycling) ) THEN
 323        grid%itimestep=0
 324    ENDIF
 325
 326    IF ( config_flags%restart .or. grid%moved .or. config_flags%hrrr_cycling) THEN
 327        first_trip_for_this_domain = .TRUE.
 328    ENDIF
 329
 330    CALL INITIALIZE_STOCH  (grid, config_flags,                  &
 331                            first_trip_for_this_domain,          &
 332                            ips, ipe, jps, jpe, kps, kpe,        &
 333                            ids, ide, jds, jde, kds, kde,        &
 334                            ims, ime, jms, jme, kms, kme,        &
 335                            its, ite, jts, jte, kts, kte,        &
 336                            imsx, imex, jmsx, jmex, kmsx, kmex,  &
 337                            ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
 338                            imsy, imey, jmsy, jmey, kmsy, kmey,  &
 339                            ipsy, ipey, jpsy, jpey, kpsy, kpey   )
 340
 341 !   --- END SETUP STOCHASTIC PERTURBATION SCHEMES ----------
 342
 343
 344
 345 ! wig: Add a combined exponential+linear weight on the mother boundaries
 346 !      following code changes by Ruby Leung. For the nested grid, there
 347 !      appears to be some problems when a sponge is used. The points where
 348 !      processors meet have problematic values.
 349
 350    CALL lbc_fcx_gcx ( grid%fcx , grid%gcx , grid%spec_bdy_width , &
 351                       grid%spec_zone , grid%relax_zone , grid%dt , config_flags%spec_exp , &
 352                       config_flags%specified , config_flags%nested )
 353
 354    IF ( config_flags%nested ) THEN
 355      grid%dtbc = 0.
 356    ENDIF
 357
 358    IF ( ( grid%id .NE. 1 ) .AND. ( .NOT. config_flags%input_from_file ) ) THEN
 359
 360       !  Every time a domain starts or every time a domain moves, this routine is called.  We want
 361       !  the center (middle) lat/lon of the grid for the metacode.  The lat/lon values are
 362       !  defined at mass points.  Depending on the even/odd points in the SN and WE directions,
 363       !  we end up with the middle point as either 1 point or an average of either 2 or 4 points.
 364       !  Add to this, the need to make sure that we are on the correct patch to retrieve the
 365       !  value of the lat/lon, AND that the lat/lons (for an average) may not all be on the same
 366       !  patch.  Once we find the correct value for lat lon, we need to keep it around on all patches,
 367       !  which is where the wrf_dm_min_real calls come in.
 368       !  If this is the most coarse domain, we do not go in here.  Also, if there is an input file
 369       !  (which has the right values for the middle lat/lon) we do not go in this IF test.
 370
 371       IF      ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
 372          num_points_lat_lon = 1
 373          iloc = ide/2
 374          jloc = jde/2
 375          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 376               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 377             lat1 = grid%xlat (iloc,jloc)
 378             lon1 = grid%xlong(iloc,jloc)
 379          ELSE
 380             lat1 = 99999.
 381             lon1 = 99999.
 382          END IF
 383          lat1 = wrf_dm_min_real ( lat1 )
 384          lon1 = wrf_dm_min_real ( lon1 )
 385          CALL nl_set_cen_lat ( grid%id , lat1 )
 386          CALL nl_set_cen_lon ( grid%id , lon1 )
 387       ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
 388          num_points_lat_lon = 2
 389          iloc = (ide-1)/2
 390          jloc =  jde   /2
 391          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 392               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 393             lat1 = grid%xlat (iloc,jloc)
 394             lon1 = grid%xlong(iloc,jloc)
 395          ELSE
 396             lat1 = 99999.
 397             lon1 = 99999.
 398          END IF
 399          lat1 = wrf_dm_min_real ( lat1 )
 400          lon1 = wrf_dm_min_real ( lon1 )
 401
 402          iloc = (ide+1)/2
 403          jloc =  jde   /2
 404          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 405               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 406             lat2 = grid%xlat (iloc,jloc)
 407             lon2 = grid%xlong(iloc,jloc)
 408          ELSE
 409             lat2 = 99999.
 410             lon2 = 99999.
 411          END IF
 412          lat2 = wrf_dm_min_real ( lat2 )
 413          lon2 = wrf_dm_min_real ( lon2 )
 414
 415          CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
 416          CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
 417       ELSE IF ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
 418          num_points_lat_lon = 2
 419          iloc =  ide   /2
 420          jloc = (jde-1)/2
 421          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 422               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 423             lat1 = grid%xlat (iloc,jloc)
 424             lon1 = grid%xlong(iloc,jloc)
 425          ELSE
 426             lat1 = 99999.
 427             lon1 = 99999.
 428          END IF
 429          lat1 = wrf_dm_min_real ( lat1 )
 430          lon1 = wrf_dm_min_real ( lon1 )
 431
 432          iloc =  ide   /2
 433          jloc = (jde+1)/2
 434          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 435               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 436             lat2 = grid%xlat (iloc,jloc)
 437             lon2 = grid%xlong(iloc,jloc)
 438          ELSE
 439             lat2 = 99999.
 440             lon2 = 99999.
 441          END IF
 442          lat2 = wrf_dm_min_real ( lat2 )
 443          lon2 = wrf_dm_min_real ( lon2 )
 444
 445          CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
 446          CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
 447       ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
 448          num_points_lat_lon = 4
 449          iloc = (ide-1)/2
 450          jloc = (jde-1)/2
 451          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 452               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 453             lat1 = grid%xlat (iloc,jloc)
 454             lon1 = grid%xlong(iloc,jloc)
 455          ELSE
 456             lat1 = 99999.
 457             lon1 = 99999.
 458          END IF
 459          lat1 = wrf_dm_min_real ( lat1 )
 460          lon1 = wrf_dm_min_real ( lon1 )
 461
 462          iloc = (ide+1)/2
 463          jloc = (jde-1)/2
 464          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 465               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 466             lat2 = grid%xlat (iloc,jloc)
 467             lon2 = grid%xlong(iloc,jloc)
 468          ELSE
 469             lat2 = 99999.
 470             lon2 = 99999.
 471          END IF
 472          lat2 = wrf_dm_min_real ( lat2 )
 473          lon2 = wrf_dm_min_real ( lon2 )
 474
 475          iloc = (ide-1)/2
 476          jloc = (jde+1)/2
 477          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 478               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 479             lat3 = grid%xlat (iloc,jloc)
 480             lon3 = grid%xlong(iloc,jloc)
 481          ELSE
 482             lat3 = 99999.
 483             lon3 = 99999.
 484          END IF
 485          lat3 = wrf_dm_min_real ( lat3 )
 486          lon3 = wrf_dm_min_real ( lon3 )
 487
 488          iloc = (ide+1)/2
 489          jloc = (jde+1)/2
 490          IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
 491               ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
 492             lat4 = grid%xlat (iloc,jloc)
 493             lon4 = grid%xlong(iloc,jloc)
 494          ELSE
 495             lat4 = 99999.
 496             lon4 = 99999.
 497          END IF
 498          lat4 = wrf_dm_min_real ( lat4 )
 499          lon4 = wrf_dm_min_real ( lon4 )
 500
 501          CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 + lat3 + lat4 ) * 0.25 )
 502          CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 + lon3 + lon4 ) * 0.25 )
 503       END IF
 504    END IF
 505
 506    !  Make sure that the base-state vales have not changed between what is in the input file
 507    !  and the namelist file.
 508
 509    IF ( .NOT. grid%this_is_an_ideal_run ) THEN
 510       CALL nl_get_p_top_requested  ( 1 , p_top_test )
 511       IF ( grid%p_top .NE. p_top_test ) THEN
 512          CALL wrf_error_fatal ( 'start_em: p_top from the namelist does not match p_top from the input file.' )
 513       END IF
 514    END IF
 515
 516    IF ( config_flags%use_baseparam_fr_nml ) then
 517       CALL nl_get_base_pres        ( 1 , p00        )
 518       CALL nl_get_base_temp        ( 1 , t00        )
 519       CALL nl_get_base_lapse       ( 1 , a          )
 520       CALL nl_get_iso_temp         ( 1 , tiso       )
 521       CALL nl_get_base_lapse_strat ( 1 , a_strat    )
 522       CALL nl_get_base_pres_strat  ( 1 , p_strat    )
 523       IF ( ( t00 .LT. 100. .or. p00 .LT. 10000.) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 524          WRITE(wrf_err_message,*) 'start_em: BAD BASE STATE for T00 or P00 in namelist.input file'
 525          CALL wrf_error_fatal(TRIM(wrf_err_message))
 526       END IF
 527
 528    ELSE
 529    ! get these constants from model data
 530
 531       t00     = grid%t00
 532       p00     = grid%p00
 533       a       = grid%tlp
 534       tiso    = grid%tiso
 535       a_strat = grid%tlp_strat
 536       p_strat = grid%p_strat
 537       IF ( ( t00 .LT. 100. .or. p00 .LT. 10000.) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 538          WRITE(wrf_err_message,*)&
 539          'start_em: did not find base state parameters in wrfinput. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
 540          CALL wrf_error_fatal(TRIM(wrf_err_message))
 541       ENDIF
 542
 543    ENDIF
 544
 545 !  check if tiso in the input file agrees with namelist value
 546
 547    CALL nl_get_iso_temp   ( 1 , tiso_tmp )
 548    IF ( ( tiso_tmp .NE. tiso ) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 549       WRITE(wrf_err_message,*)&
 550       'start_em: namelist iso_temp is not equal to iso_temp in wrfinput. Reset nml value and rerun'
 551       CALL wrf_error_fatal(TRIM(wrf_err_message))
 552    ENDIF
 553
 554    IF ( .NOT. config_flags%restart .AND. &
 555         (( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ))) THEN
 556
 557       IF ( config_flags%map_proj .EQ. 0 ) THEN
 558          CALL wrf_error_fatal ( 'start_domain: Idealized case cannot have a separate nested input file' )
 559       END IF
 560
 561       !  Base state potential temperature and inverse density (alpha = 1/rho) from
 562       !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
 563       !  from equation of state.  The potential temperature is a perturbation from t0.
 564
 565       IF ( grid%c1f(1) .NE. 1. ) THEN
 566          CALL wrf_debug ( 0 , '---- WARNING : Maybe old non-HVC input, setting default 1d array values for TF' )
 567          IF ( grid%hybrid_opt .NE. 0 ) THEN
 568             CALL wrf_error_fatal ( '---- Error : Cannot use old input and try to use hybrid vertical coordinate option' )
 569          END IF
 570          DO k = 1, kte
 571             grid%c1f(k) = 1.
 572             grid%c2f(k) = 0.
 573             grid%c3f(k) = grid%znw(k)
 574             grid%c4f(k) = 0.
 575             grid%c1h(k) = 1.
 576             grid%c2h(k) = 0.
 577             grid%c3h(k) = grid%znu(k)
 578             grid%c4h(k) = 0.
 579          END DO
 580       END IF
 581
 582       !  With newer versions of WRF, the grid%t_2 field might come in as "zero"
 583       !  from an older data set. Yikes, but it can be corrected unless the user
 584       !  has requested moist theta.
 585
 586       IF ( grid%t_2(its,kte-1,jts) .EQ. 0. ) THEN
 587          CALL wrf_debug ( 0 , '---- WARNING : Older v3 input data detected' )
 588          IF ( grid%use_theta_m .NE. 0 ) THEN
 589             CALL wrf_error_fatal ( '---- Error : Cannot use moist theta option with old data' )
 590          END IF
 591          DO j = jts, MIN(jte,jde-1)
 592             DO k = 1, kte-1
 593                DO i = its, MIN(ite,ide-1)
 594                   grid%t_2(i,k,j) = grid%th_phy_m_t0(i,k,j)
 595                END DO
 596             END DO
 597          END DO
 598       END IF
 599
 600       DO j = jts, MIN(jte,jde-1)
 601          DO i = its, MIN(ite,ide-1)
 602             !  Base state pressure is a function of eta level and terrain, only, plus
 603             !  the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
 604             !  temperature, K), A (temperature difference, from 1000 mb to 300 mb, K),
 605             !  tiso (isothermal temperature at tropopause/lower stratosphere),
 606             !  p_strat (pressure at top of isothermal layer), A_strat (lapse rate in
 607             !  stratosphere above isothermal layer)
 608             p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
 609             DO k = 1, kte-1
 610                grid%pb(i,k,j) = grid%c3h(k)*(p_surf - grid%p_top) + grid%c4h(k) + grid%p_top
 611                temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
 612                IF ( grid%pb(i,k,j) .LT. p_strat ) THEN
 613                   temp = tiso + A_strat * LOG ( grid%pb(i,k,j)/p_strat )
 614                ENDIF
 615                grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
 616                grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 617             END DO
 618             !  Base state mu is defined as base state surface pressure minus grid%p_top
 619             grid%MUB(i,j) = p_surf - grid%p_top
 620             !  Integrate base geopotential, starting at terrain elevation.  This assures that
 621             !  the base state is in exact hydrostatic balance with respect to the model equations.
 622             !  This field is on full levels.
 623             grid%phb(i,1,j) = grid%ht(i,j) * g
 624             IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 625                DO kk  = 2,kte
 626                   k = kk - 1
 627                   grid%phb(i,kk,j) = grid%phb(i,kk-1,j) - grid%dnw(kk-1)*(grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))*grid%alb(i,kk-1,j)
 628                END DO
 629             ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 630                DO k  = 2,kte
 631                   pfu = grid%c3f(k  )*grid%MUB(i,j) + grid%c4f(k  ) + grid%p_top
 632                   pfd = grid%c3f(k-1)*grid%MUB(i,j) + grid%c4f(k-1) + grid%p_top
 633                   phm = grid%c3h(k-1)*grid%MUB(i,j) + grid%c4h(k-1) + grid%p_top
 634                   grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
 635                END DO
 636             END IF
 637          END DO
 638       END DO
 639    ENDIF
 640    IF(.not.config_flags%restart)THEN
 641 !  if this is for a nested domain, the defined/interpolated fields are the _2
 642      IF ( first_trip_for_this_domain ) THEN
 643 ! data that is expected to be zero must be explicitly initialized as such
 644 !    grid%h_diabatic = 0.
 645        DO j = jts,min(jte,jde-1)
 646        DO k = kts,kte-1
 647        DO i = its, min(ite,ide-1)
 648            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 649              grid%t_1(i,k,j)=grid%t_2(i,k,j)
 650            ENDIF
 651        ENDDO
 652        ENDDO
 653        ENDDO
 654        DO j = jts,min(jte,jde-1)
 655        DO k = kts,kte
 656        DO i = its, min(ite,ide-1)
 657           grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
 658        ENDDO
 659        ENDDO
 660        ENDDO
 661        DO j = jts,min(jte,jde-1)
 662        DO i = its, min(ite,ide-1)
 663            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 664              grid%MU_1(i,j)=grid%MU_2(i,j)
 665            ENDIF
 666        ENDDO
 667        ENDDO
 668      END IF
 669 !  reconstitute base-state fields
 670     IF(config_flags%max_dom .EQ. 1)THEN
 671 ! with single domain, grid%t_init from wrfinput is OK to use
 672      DO j = jts,min(jte,jde-1)
 673      DO k = kts,kte-1
 674      DO i = its, min(ite,ide-1)
 675        IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 676          grid%pb(i,k,j) = grid%c3h(k  )*grid%MUB(i,j) + grid%c4h(k  ) + grid%p_top
 677          grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 678        ENDIF
 679      ENDDO
 680      ENDDO
 681      ENDDO
 682     ELSE ! there is more than 1 domain
 683      IF ( .NOT. grid%this_is_an_ideal_run ) THEN ! this is a real run
 684         DO j = jts,min(jte,jde-1)
 685         DO k = kts,kte-1
 686         DO i = its, min(ite,ide-1)
 687           IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 688             grid%pb(i,k,j) = grid%c3h(k  )*grid%MUB(i,j) + grid%c4h(k  ) + grid%p_top
 689             temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
 690             IF ( grid%pb(i,k,j) .LT. p_strat ) THEN
 691                temp = tiso + A_strat * LOG ( grid%pb(i,k,j)/p_strat )
 692             ENDIF
 693             grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
 694             grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 695        ENDIF
 696       ENDDO
 697       ENDDO
 698       ENDDO
 699       !Note that rebalance must be set to 1 for vertical nesting
 700       IF ( config_flags%rebalance .EQ. 1 ) THEN
 701         !  Integrate base geopotential, starting at terrain elevation.  This assures that
 702         !  the base state is in exact hydrostatic balance with respect to the model equations.
 703         !  This field is on full levels.
 704        DO j = jts,min(jte,jde-1)
 705        DO i = its, min(ite,ide-1)
 706         grid%phb(i,1,j) = grid%ht(i,j) * g
 707         IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 708            DO kk  = 2,kte
 709               k = kk - 1
 710               grid%phb(i,kk,j) = grid%phb(i,kk-1,j) - grid%dnw(kk-1)*(grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))*grid%alb(i,kk-1,j)
 711            END DO
 712         ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 713            DO k  = 2,kte
 714               pfu = grid%c3f(k  )*grid%MUB(i,j) + grid%c4f(k  ) + grid%p_top
 715               pfd = grid%c3f(k-1)*grid%MUB(i,j) + grid%c4f(k-1) + grid%p_top
 716               phm = grid%c3h(k-1)*grid%MUB(i,j) + grid%c4h(k-1) + grid%p_top
 717               grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
 718             END DO
 719           ENDIF
 720        ENDDO
 721        ENDDO
 722        END IF
 723      ELSE
 724         IF ( config_flags%ideal_init_method .EQ. 1 ) THEN
 725            DO j = jts,min(jte,jde-1)
 726            DO k = kts,kte-1
 727            DO i = its, min(ite,ide-1)
 728              IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 729                grid%pb(i,k,j) = grid%c3h(k  )*grid%MUB(i,j) + grid%c4h(k  ) + grid%p_top
 730                grid%alb(i,k,j) = -grid%rdnw(k)*(grid%phb(i,k+1,j)-grid%phb(i,k,j))/(grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))
 731                grid%t_init(i,k,j) = grid%alb(i,k,j)*(p1000mb/r_d)/((grid%pb(i,k,j)/p1000mb)**cvpm) - t0
 732              ENDIF
 733            ENDDO
 734            ENDDO
 735            ENDDO
 736         ELSE IF ( config_flags%ideal_init_method .EQ. 2 ) THEN
 737            DO j = jts,min(jte,jde-1)
 738            DO k = kts,kte-1
 739            DO i = its, min(ite,ide-1)
 740              IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 741                grid%pb(i,k,j) = grid%c3h(k  )*grid%MUB(i,j) + grid%c4h(k  ) + grid%p_top
 742                grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 743              ENDIF
 744            ENDDO
 745            ENDDO
 746            ENDDO
 747            !  Integrate base geopotential, starting at terrain elevation.  This assures that
 748            !  the base state is in exact hydrostatic balance with respect to the model equations.
 749            !  This field is on full levels.
 750            DO j = jts,min(jte,jde-1)
 751            DO i = its, min(ite,ide-1)
 752             grid%phb(i,1,j) = grid%ht(i,j) * g
 753             DO kk  = 2,kte
 754               k = kk - 1
 755               grid%phb(i,kk,j) = grid%phb(i,kk-1,j) - grid%dnw(kk-1)*(grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))*grid%alb(i,kk-1,j)
 756             END DO
 757            ENDDO
 758            ENDDO
 759         END IF  ! end if initialization for ideal, use method 1 or 2
 760      END IF ! end if this is a real or ideal run
 761     END IF ! end if there is more than 1 domain
 762
 763 !------base state is finished, perturbation values are recalculated below-----
 764   ! For the HRRR application, we want to go through this rebalancing.  They
 765   ! have a clunky way of testig on this, but it is operational, so no need to
 766   ! cause them troubles.
 767   IF ( ( grid%dfi_opt .EQ. DFI_NODFI ) .and. ( config_flags%cycling ) ) THEN
 768        call rebalance_driver_cycl (grid )
 769        DO j = jts,min(jte,jde-1)
 770        DO k = kts,kte
 771        DO i = its, min(ite,ide-1)
 772           grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
 773        ENDDO
 774        ENDDO
 775        ENDDO
 776   !  Everyone else runs thorugh here.  NOTE that we might also want to call the
 777   !  rebalancing ourselves for different set ups.  It does not impact HRRR, so
 778   !  it is an easy decision.
 779   ELSE
 780      ! We request rebalancing for vertical grid nesting, or when the user asks for rebalancing.
 781      ! (Note that rebalance now must be set to 1 for vertical nesting)
 782      ! Rebalance recomputes 1/rho, p, ph_2, ph0, p_hyd
 783      IF ( config_flags%rebalance .EQ. 1 ) THEN
 784        call rebalance_driver_cycl (grid )
 785        DO j = jts,min(jte,jde-1)
 786        DO k = kts,kte
 787        DO i = its, min(ite,ide-1)
 788           grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
 789        ENDDO
 790        ENDDO
 791        ENDDO
 792      END IF
 793      ! Use equations from calc_p_rho_phi to derive p and al from ph: linear in log p
 794     IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 795        DO j=jts,min(jte,jde-1)
 796        DO k=kts,kte-1
 797        DO i=its,min(ite,ide-1)
 798           grid%al(i,k,j)=-1./((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_1(i,j)))*(grid%alb(i,k,j)*(grid%c1h(k)*grid%mu_1(i,j))  &
 799                      +grid%rdnw(k)*(grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)))
 800        ENDDO
 801        ENDDO
 802        ENDDO
 803     ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 804        DO j=jts,min(jte,jde-1)
 805        DO k=kts,kte-1
 806        DO i=its,min(ite,ide-1)
 807           pfu = grid%c3f(k+1)*(grid%MUB(i,j)+grid%MU_1(i,j)) + grid%c4f(k+1) + grid%p_top
 808           pfd = grid%c3f(k  )*(grid%MUB(i,j)+grid%MU_1(i,j)) + grid%c4f(k  ) + grid%p_top
 809           phm = grid%c3h(k  )*(grid%MUB(i,j)+grid%MU_1(i,j)) + grid%c4h(k  ) + grid%p_top
 810           grid%al(i,k,j) = (grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)+grid%phb(i,k+1,j)-grid%phb(i,k,j)) &
 811                             /phm/LOG(pfd/pfu)-grid%alb(i,k,j)
 812        ENDDO
 813        ENDDO
 814        ENDDO
 815     END IF ! which hypsometric option
 816     DO j=jts,min(jte,jde-1)
 817     DO k=kts,kte-1
 818     DO i=its,min(ite,ide-1)
 819 #if ( WRFPLUS == 1 )
 820        IF ( .NOT. config_flags%var4d_run ) THEN
 821           IF ( config_flags%use_theta_m == 0 ) THEN
 822              qvf = 1.+rvovrd*moist(i,k,j,P_QV)
 823           ELSE
 824              qvf = 1
 825           END IF
 826           grid%p(i,k,j)=p1000mb*( (r_d*(t0+grid%t_1(i,k,j))*qvf)/                     &
 827                      (p1000mb*(grid%al(i,k,j)+grid%alb(i,k,j))) )**cpovcv  &
 828                      -grid%pb(i,k,j)
 829        END IF
 830 #else
 831        IF ( config_flags%use_theta_m == 0 ) THEN
 832           qvf = 1.+rvovrd*moist(i,k,j,P_QV)
 833        ELSE
 834           qvf = 1
 835        END IF
 836        grid%p(i,k,j)=p1000mb*( (r_d*(t0+grid%t_1(i,k,j))*qvf)/                     &
 837                   (p1000mb*(grid%al(i,k,j)+grid%alb(i,k,j))) )**cpovcv  &
 838                   -grid%pb(i,k,j)
 839 #endif
 840        grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
 841        grid%alt(i,k,j) = grid%al(i,k,j) + grid%alb(i,k,j)
 842     ENDDO
 843     ENDDO
 844     ENDDO
 845   ENDIF ! Various rebalancing options
 846     IF ( .NOT. grid%this_is_an_ideal_run ) THEN
 847        DO j=jts,min(jte,jde-1)
 848           DO i=its,min(ite,ide-1)
 849              p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
 850              grid%p_hyd_w(i,1,j) = grid%p(i,1,j) + grid%c3f(1)*(p_surf - grid%p_top) + grid%c4f(1) + grid%p_top
 851              DO k=kts+1,kte
 852                 grid%p_hyd_w(i,k,j) = ( 2.*(grid%p(i,k-1,j)+grid%pb(i,k-1,j)) - grid%p_hyd_w(i,k-1,j) )
 853              ENDDO
 854           ENDDO
 855        ENDDO
 856     ELSE
 857        DO j=jts,min(jte,jde-1)
 858           DO i=its,min(ite,ide-1)
 859              p_surf = grid%MUB(i,j)+grid%p_top
 860              grid%p_hyd_w(i,1,j) = grid%p(i,1,j) + grid%c3f(1)*(p_surf - grid%p_top) + grid%c4f(1) + grid%p_top
 861              DO k=kts+1,kte
 862                 grid%p_hyd_w(i,k,j) = ( 2.*(grid%p(i,k-1,j)+grid%pb(i,k-1,j)) - grid%p_hyd_w(i,k-1,j) )
 863              ENDDO
 864           ENDDO
 865        ENDDO
 866     END IF
 867
 868    ENDIF ! first trip for this domain
 869
 870    DO j=jts,min(jte,jde-1)
 871       DO k = kts,kte
 872          DO i = its, min(ite,ide-1)
 873             z_at_q(i,k,j)=(grid%ph_2(i,k,j)+grid%phb(i,k,j))/g
 874          ENDDO
 875       ENDDO
 876    ENDDO
 877
 878    IF ( grid%press_adj .and. ( grid%id .NE. 1 ) .AND. .NOT. ( config_flags%restart ) .AND. &
 879        ( ( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ) ) ) THEN
 880       DO j = jts, MIN(jte,jde-1)
 881          DO i = its, MIN(ite,ide-1)
 882             grid%MU_2(i,j) = grid%MU_2(i,j) + grid%al(i,1,j) / ( grid%alt(i,1,j) * grid%alb(i,1,j) ) * &
 883                                     g * ( grid%ht(i,j) - grid%ht_fine(i,j) )
 884          END DO
 885        END DO
 886        DO j = jts,min(jte,jde-1)
 887        DO i = its, min(ite,ide-1)
 888           grid%MU_1(i,j)=grid%MU_2(i,j)
 889        ENDDO
 890        ENDDO
 891
 892    END IF
 893
 894 ! set GMT outside of phy_init because phy_init may not be called on this
 895 ! process if, for example, it is a moving nest and if this part of the domain is not
 896 ! being initialized (not the leading edge).
 897    CALL domain_setgmtetc( grid, start_of_simulation )
 898    IF ( first_trip_for_this_domain ) THEN
 899
 900    CALL wrf_debug ( 100 , 'start_domain_em: Before call to phy_init' )
 901
 902 ! namelist MPDT does not exist yet, so set it here
 903 ! MPDT is the call frequency for microphysics in minutes (0 means every step)
 904    MPDT = 0.
 905
 906    IF(config_flags%cycling) THEN
 907        start_of_simulation = .true.
 908 !  print *,'cycling, start_of_simulation -->',config_flags%cycling, start_of_simulation
 909        grid%xtime=0.
 910 !   print *,'xtime=',grid%xtime
 911    ENDIF
 912 !-----------------------------------------------------------------------------
 913 ! Adaptive time step: Added by T. Hutchinson, WSI  11/6/07
 914 !
 915 !
 916
 917    IF ( ( grid%use_adaptive_time_step ) .AND. &
 918         ( ( grid%dfi_opt .EQ. DFI_NODFI ) .OR. ( grid%dfi_stage .EQ. DFI_FST ) ) ) THEN
 919
 920 !  Calculate any variables that were not set
 921 !BPR BEGIN
 922 !     This subroutine is called more than once at the first time step for a
 923 !     given domain.  The following if statement is to prevent the code in
 924 !     the if statement from executing more than once per domain at the
 925 !     beginning of the model run (since last_step_update=-1 the first time
 926 !     this is reached and should be =0 after this).
 927 !     Without this if statement, when this code was executed for the second
 928 !     time it can result in grid%dt being set incorrectly.
 929 !     -This is because grid%dt will be set equal to grid%starting_time_step
 930 !      which ignores a possible denominator in the starting time step.
 931 !     -The first time this code is reached is also does this, but then the
 932 !      call to adapt_timestep correct this
 933 !     -Subsequent times this code is reached adapt_timestep will not correct
 934 !      this because it will recognize that it has already been executed for
 935 !      this timestep and exit out before doing the calculation.
 936
 937       if (grid%last_step_updated .NE. grid%itimestep) then
 938
 939       if (grid%starting_time_step == -1) then
 940          grid%starting_time_step = NINT(4 * MIN(grid%dx,grid%dy) / 1000)
 941       endif
 942
 943       grid%time_step = grid%starting_time_step
 944       config_flags%time_step = grid%starting_time_step
 945       model_config_rec%time_step = grid%starting_time_step
 946
 947       if (grid%max_time_step == -1) then
 948          grid%max_time_step = NINT(8 * MIN(grid%dx,grid%dy) / 1000)
 949       endif
 950
 951       if (grid%min_time_step == -1) then
 952          grid%min_time_step = NINT(3 * MIN(grid%dx,grid%dy) / 1000)
 953       endif
 954
 955 !     Set a starting timestep.
 956
 957       grid%dt = grid%starting_time_step
 958
 959 !     Initialize max cfl values
 960
 961       grid%last_max_vert_cfl = 0
 962       grid%last_max_horiz_cfl = 0
 963
 964 !     Check to assure that time_step_sound is to be dynamically set.
 965
 966       CALL nl_set_time_step_sound ( 1 , 0 )
 967       grid%time_step_sound = 0
 968
 969       grid%max_msftx=MAXVAL(grid%msftx)
 970       grid%max_msfty=MAXVAL(grid%msfty)
 971 #ifdef DM_PARALLEL
 972       CALL wrf_dm_maxval(grid%max_msftx, idex, jdex)
 973       CALL wrf_dm_maxval(grid%max_msfty, idex, jdex)
 974 #endif
 975       end if
 976 !BPR END
 977
 978 !     This first call just initializes variables.
 979 !     If a restart, get initialized variables from restart file
 980
 981       IF ( .NOT. ( config_flags%restart ) ) then
 982          CALL adapt_timestep(grid, config_flags)
 983       END IF
 984
 985
 986    END IF  ! adaptive computations
 987
 988 ! End of adaptive time step modifications
 989 !-----------------------------------------------------------------------------
 990
 991
 992    CALL set_tiles ( grid , grid%imask_nostag, ims, ime, jms, jme, ips, ipe, jps, jpe )
 993 !
 994 ! Phy init can do reads and broadcasts when initializing physics -- landuse for example. However, if
 995 ! we're running on a reduced mesh (that is, some tasks don't have any work) we have to at least let them
 996 ! pass through this code so the broadcasts don't hang on the other, active tasks.  Set the number of
 997 ! tiles to a minimum of 1 and assume that the backwards patch ranges (ips=0, ipe=-1) will prevent
 998 ! anything else from happening on the blank tasks.  JM 20080605
 999 !
1000    if ( allowed_to_read ) grid%num_tiles = max(1,grid%num_tiles)
1001 !
1002 ! Phy_init is not necessarily thread-safe; do not multi-thread this loop.
1003 ! The tiling is to handle the fact that we may be masking off part of the computation.
1004 !
1005
1006 #ifdef DM_PARALLEL
1007 if(config_flags%sf_surface_physics.eq.NOAHMPSCHEME.and.config_flags%opt_run.eq.5)then
1008 # include "HALO_EM_HYDRO_NOAHMP_INIT.inc"
1009 endif
1010 #endif
1011
1012     IF ( ( grid%dfi_opt .EQ. DFI_NODFI ) .or. &
1013           ( ( grid%dfi_stage .NE. DFI_BCK ) .and. &
1014           ( grid%dfi_stage .NE. DFI_STARTBCK ) ) ) THEN
1015
1016    DO ij = 1, grid%num_tiles
1017
1018      CALL phy_init (  grid%id , config_flags, grid%DT, grid%RESTART, grid%znw, grid%znu,   &
1019                       grid%p_top, grid%tsk, grid%RADT,grid%BLDT,grid%CUDT, MPDT, &
1020                       grid%rucuten, grid%rvcuten, grid%rthcuten,             &
1021                       grid%rqvcuten, grid%rqrcuten, grid%rqccuten,           &
1022                       grid%rqscuten, grid%rqicuten,                          &
1023                       grid%rushten, grid%rvshten, grid%rthshten,             &
1024                       grid%rqvshten, grid%rqrshten, grid%rqcshten,           &
1025                       grid%rqsshten, grid%rqishten, grid%rqgshten,           &
1026                       grid%rublten,grid%rvblten,grid%rthblten,               &
1027                       grid%rqvblten,grid%rqcblten,grid%rqiblten,             &
1028                       grid%rthraten,grid%rthratenlw,grid%rthratensw,         &
1029                       grid%this_is_an_ideal_run,                             &
1030                       !BSINGH - For WRFCuP scheme(11/12/2013)
1031                       grid%cupflag,grid%cldfra_cup,grid%cldfratend_cup,      & !wig, 18-Sep-2006
1032                       grid%shall,                                            & !wig, 18-Sep-2006
1033                       grid%tcloud_cup,                                       & !rce, 18-apr-2012
1034                       !BSINGH - ENDS
1035                       grid%stepbl,grid%stepra,grid%stepcu,                   &
1036                       grid%w0avg, grid%rainnc, grid%rainc, grid%raincv, grid%rainncv,  &
1037                       grid%snownc, grid%snowncv, grid%graupelnc, grid%graupelncv,  &
1038                       z_at_q, grid%alt, grid%qnwfa2d, grid%qnbca2d,     & ! G. Thompson
1039                       scalar(ims,kms,jms,1), num_scalar,                & ! G. Thompson
1040                       grid%re_cloud, grid%re_ice, grid%re_snow,         & ! G. Thompson
1041                       grid%has_reqc, grid%has_reqi, grid%has_reqs,      & ! G. Thompson
1042                       grid%phb,grid%ph_2,grid%p,grid%pb,                & ! for ntu3m
1043                       moist(:,:,:,P_QV),config_flags%CCNTY,             & ! for ntu3m
1044                       scalar(:,:,:,P_QDCN),scalar(:,:,:,P_QTCN),        & ! for ntu3m
1045                       scalar(:,:,:,P_QCCN),scalar(:,:,:,P_QRCN),        & ! for ntu3m
1046                       scalar(:,:,:,P_QNIN),                             & ! for ntu3m
1047                       grid%re_cloud_gsfc, grid%re_ice_gsfc,             & ! Goddard
1048                       grid%re_snow_gsfc, grid%re_graupel_gsfc,          & ! Goddard
1049                       grid%re_hail_gsfc, grid%re_rain_gsfc,             & ! Goddard
1050                       grid%nca,grid%swrad_scat,                         &
1051                       grid%cldefi,grid%lowlyr,                          &
1052                       grid%mass_flux,                                   &
1053                       grid%rthften, grid%rqvften,                       &
1054 #if (WRF_USE_CLM == 1) && ( WRF_CHEM == 1)
1055                       model_config_rec%megan_specifier, model_config_rec%megan_factors_file,    &
1056                       model_config_rec%megan_mapped_emisfctrs,                      &
1057 #endif
1058                       grid%cldfra,                                      &
1059 #if (WRF_CHEM == 1)
1060                       grid%cldfra_old,                                  &
1061 #else
1062                       cldfra_old,                                       &
1063 #endif
1064                       grid%glw,grid%gsw,grid%emiss,grid%embck,            &
1065                       grid%lu_index,                                      &
1066                       grid%landuse_ISICE, grid%landuse_LUCATS,            &
1067                       grid%landuse_LUSEAS, grid%landuse_ISN,              &
1068                       grid%lu_state,                                      &
1069                       grid%xlat,grid%xlong,grid%xlong_u,grid%xlat_v,grid%albedo,grid%albbck,grid%GMT,grid%JULYR,grid%JULDAY,     &
1070                       grid%levsiz, num_ozmixm, num_aerosolc, grid%paerlev,  &
1071                       grid%alevsiz, grid%no_src_types,                      &
1072                       grid%tmn,grid%xland,grid%znt,grid%z0,grid%ust,grid%mol,grid%pblh,grid%tke_pbl,    &
1073                       grid%exch_h,grid%thc,grid%snowc,grid%mavail,grid%hfx,grid%qfx,grid%rainbl, &
1074                       grid%tslb,grid%zs,grid%dzs,config_flags%num_soil_layers,grid%warm_rain,  &
1075                       grid%adv_moist_cond, grid%is_CAMMGMP_used,                               &
1076                       grid%apr_gr,grid%apr_w,grid%apr_mc,grid%apr_st,grid%apr_as,      &
1077                       grid%apr_capma,grid%apr_capme,grid%apr_capmi,          &
1078                       grid%xice,grid%xicem,grid%vegfra,grid%snow,grid%canwat,grid%smstav,         &
1079                       grid%smstot, grid%sfcrunoff,grid%udrunoff,grid%grdflx,grid%acsnow,      &
1080                       grid%acsnom,grid%ivgtyp,grid%isltyp, grid%sfcevp,grid%smois,     &
1081                       grid%sh2o, grid%snowh, grid%smfr3d,                    &
1082                       grid%snoalb,                 &
1083                       grid%DX,grid%DY,grid%dx2d,grid%area2d,              &
1084                       grid%f_ice_phy,grid%f_rain_phy,grid%f_rimef_phy, &
1085                       grid%mp_restart_state,grid%tbpvs_state,grid%tbpvs0_state,&
1086                       allowed_to_read, grid%moved, start_of_simulation,               &
1087                       grid%LAGDAY,                            &
1088                       grid%U10,grid%V10,grid%U10E,grid%V10E,  &
1089                       ids, ide, jds, jde, kds, kde,           &
1090                       ims, ime, jms, jme, kms, kme,           &
1091                       grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, &
1092                       config_flags%num_urban_ndm,                           & !multi-layer urban
1093                       config_flags%urban_map_zrd,                           & !multi-layer urban
1094                       config_flags%urban_map_zwd,                           & !multi-layer urban
1095                       config_flags%urban_map_gd,                            & !multi-layer urban
1096                       config_flags%urban_map_zd,                            & !multi-layer urban
1097                       config_flags%urban_map_zdf,                           & !multi-layer urban
1098                       config_flags%urban_map_bd,                            & !multi-layer urban
1099                       config_flags%urban_map_wd,                            & !multi-layer urban
1100                       config_flags%urban_map_gbd,                           & !multi-layer urban
1101                       config_flags%urban_map_fbd,                           & !multi-layer urban
1102                       config_flags%urban_map_zgrd,                          & !multi-layer urban
1103                       config_flags%num_urban_hi,                            & !multi-layer urban
1104                       grid%raincv_a,grid%raincv_b,                                    &
1105                       grid%gd_cloud, grid%gd_cloud2,                                  & ! Optional
1106                       grid%gd_cloud_a, grid%gd_cloud2_a,                              & ! Optional
1107                       grid%QC_CU, grid%QI_CU,                                         & ! Optional
1108                       ozmixm,grid%pin,                             &     ! Optional
1109                       grid%aerodm,grid%pina,                       &     ! Optional
1110                       grid%m_ps_1,grid%m_ps_2,grid%m_hybi,aerosolc_1,aerosolc_2,&  ! Optional
1111                       grid%rundgdten,grid%rvndgdten,grid%rthndgdten,         &     ! Optional
1112                       grid%rphndgdten,grid%rqvndgdten,grid%rmundgdten,       &     ! Optional
1113                       grid%SDA_HFX, grid%SDA_QFX, grid%QNORM, grid%HFX_BOTH,grid%QFX_BOTH,      &   ! fasdas
1114                       grid%HFX_FDDA,                                                            &   ! fasdas
1115                       grid%FGDT,grid%stepfg,                        &     ! Optional
1116                       grid%cugd_tten,grid%cugd_ttens,grid%cugd_qvten,        &   ! Optional
1117                       grid%cugd_qvtens,grid%cugd_qcten,                 &   ! Optional
1118                       grid%ISNOWXY, grid%ZSNSOXY, grid%TSNOXY,                                  &   ! Optional Noah-MP
1119                       grid%SNICEXY, grid%SNLIQXY, grid%TVXY, grid%TGXY, grid%CANICEXY,          &   ! Optional Noah-MP
1120                       grid%CANLIQXY, grid%EAHXY, grid%TAHXY, grid%CMXY,                         &   ! Optional Noah-MP
1121                       grid%CHXY, grid%FWETXY, grid%SNEQVOXY, grid%ALBOLDXY, grid%QSNOWXY,       &   ! Optional Noah-MP
1122                       grid%QRAINXY,                                                             &   ! Optional Noah-MP
1123                       grid%WSLAKEXY, grid%ZWTXY, grid%WAXY, grid%WTXY, grid%LFMASSXY,           &   ! Optional Noah-MP
1124                       grid%RTMASSXY,                                                            &   ! Optional Noah-MP
1125                       grid%STMASSXY, grid%WOODXY, grid%STBLCPXY, grid%FASTCPXY,                 &   ! Optional Noah-MP
1126                       grid%GRAINXY, grid%GDDXY,                                                 &   ! Optional Noah-MP
1127                       grid%CROPTYPE, grid%CROPCAT,                                              &   ! Optional Noah-MP
1128                       config_flags%opt_crop,                                                    &
1129                       config_flags%opt_irr, config_flags%opt_irrm,                              &   ! Optional Noah-MP
1130                       grid%IRNUMSI, grid%IRNUMMI, grid%IRNUMFI, grid%IRWATSI, grid%IRWATMI,     &   ! Optional Noah-MP
1131                       grid%IRWATFI, grid%IRELOSS, grid%IRSIVOL, grid%IRMIVOL, grid%IRFIVOL,     &   ! Optional Noah-MP
1132                       grid%IRRSPLH, grid%qtdrain,                                               &   ! Optional Noah-MP
1133                       grid%gecros_state,                                                        &   ! Optional gecros crop
1134                       grid%XSAIXY,grid%LAI,                                                     &   ! Optional Noah-MP
1135                       grid%T2MVXY, grid%T2MBXY, grid%CHSTARXY,                                  &   ! Optional Noah-MP
1136                       grid%SMOISEQ  ,grid%SMCWTDXY ,grid%RECHXY, grid%DEEPRECHXY, grid%AREAXY,  & ! Optional Noah-MP
1137                       config_flags%wtddt ,grid%stepwtd ,grid%QRFSXY ,grid%QSPRINGSXY ,grid%QSLATXY, & ! Optional Noah-MP
1138                       grid%FDEPTHXY, grid%RIVERBEDXY, grid%EQZWT, grid%RIVERCONDXY, grid%PEXPXY, & ! Optional Noah-MP
1139                       grid%rechclim  ,                                                           & ! Optional Noah-MP
1140                       grid%msftx, grid%msfty,                              &
1141                       grid%DZR, grid%DZB, grid%DZG,                          & !Optional urban
1142                       grid%TR_URB2D,grid%TB_URB2D,grid%TG_URB2D,grid%TC_URB2D,    & !Optional urban
1143                       grid%QC_URB2D, grid%XXXR_URB2D,grid%XXXB_URB2D,        & !Optional urban
1144                       grid%XXXG_URB2D, grid%XXXC_URB2D,                 & !Optional urban
1145                       grid%TRL_URB3D, grid%TBL_URB3D, grid%TGL_URB3D,        & !Optional urban
1146                       grid%SH_URB2D, grid%LH_URB2D, grid%G_URB2D, grid%RN_URB2D,  & !Optional urban
1147                       grid%TS_URB2D, grid%FRC_URB2D, grid%UTYPE_URB2D,      & !Optional urban
1148                       grid%CMCR_URB2D,grid%TGR_URB2D,grid%TGRL_URB3D,grid%SMR_URB3D, & !Optional urban
1149                       grid%DRELR_URB2D,grid%DRELB_URB2D,grid%DRELG_URB2D,    & !Optional urban
1150                       grid%FLXHUMR_URB2D,grid%FLXHUMB_URB2D,grid%FLXHUMG_URB2D,   & !Optional urban
1151                       grid%TRB_URB4D,grid%TW1_URB4D,grid%TW2_URB4D,grid%TGB_URB4D,grid%TLEV_URB3D,  & !multi-layer urban
1152                       grid%QLEV_URB3D,grid%TW1LEV_URB3D,grid%TW2LEV_URB3D,        & !multi-layer urban
1153                       grid%TGLEV_URB3D,grid%TFLEV_URB3D,grid%SF_AC_URB3D,         & !multi-layer urban
1154                       grid%LF_AC_URB3D,grid%CM_AC_URB3D,grid%SFVENT_URB3D,grid%LFVENT_URB3D,   & !multi-layer urban
1155                       grid%SFWIN1_URB3D,grid%SFWIN2_URB3D, & !multi-layer urban
1156                       grid%SFW1_URB3D,grid%SFW2_URB3D,grid%SFR_URB3D,grid%SFG_URB3D,     & !multi-layer urban
1157                       grid%EP_PV_URB3D,grid%T_PV_URB3D,                                  & !GRZ
1158                       grid%TRV_URB4D,grid%QR_URB4D,grid%QGR_URB3D,grid%TGR_URB3D,        & !GRZ
1159                       grid%DRAIN_URB4D,grid%DRAINGR_URB3D,grid%SFRV_URB3D,               & !GRZ
1160                       grid%LFRV_URB3D,grid%DGR_URB3D,grid%DG_URB3D,grid%LFR_URB3D,     & !GRZ
1161                       grid%LFG_URB3D,                                                  & !GRZ
1162                       grid%LP_URB2D,grid%HI_URB2D,grid%LB_URB2D,grid%HGT_URB2D,        & !multi-layer urban
1163                       grid%MH_URB2D,grid%STDH_URB2D,grid%LF_URB2D,                     & !SLUCM
1164                       grid%A_U_BEP,grid%A_V_BEP,grid%A_T_BEP,grid%A_Q_BEP,             & !multi-layer urban
1165                       grid%A_E_BEP,grid%B_U_BEP,grid%B_V_BEP,grid%B_T_BEP,             & !multi-layer urban
1166                       grid%B_Q_BEP,grid%B_E_BEP,grid%DLG_BEP,                          & !multi-layer urban
1167                       grid%DL_U_BEP,grid%SF_BEP,grid%VL_BEP,                           & !multi-layer urban
1168                       grid%TML,grid%T0ML,grid%HML,grid%H0ML,grid%HUML,grid%HVML,grid%TMOML,     & !Optional oml
1169                       grid%lakedepth2d,  grid%savedtke12d,  grid%snowdp2d,   grid%h2osno2d,       & !lake
1170                       grid%snl2d,        grid%t_grnd2d,     grid%t_lake3d,   grid%lake_icefrac3d, & !lake
1171                       grid%z_lake3d,     grid%dz_lake3d,    grid%t_soisno3d, grid%h2osoi_ice3d,   & !lake
1172                       grid%h2osoi_liq3d, grid%h2osoi_vol3d, grid%z3d,        grid%dz3d,           & !lake
1173                       grid%zi3d,         grid%watsat3d,     grid%csol3d,     grid%tkmg3d,         & !lake
1174                       grid%tkdry3d,      grid%tksatu3d,     grid%lake2d,                            & !lake
1175                       config_flags%lakedepth_default,        config_flags%lake_min_elev, grid%lake_depth,     & !lake
1176                       grid%lakemask,        grid%lakeflag,  grid%LAKE_DEPTH_FLAG, grid%use_lakedepth,     & !lake
1177                       grid%water_depth,  grid%BATHYMETRY_FLAG, grid%shalwater_z0,grid%shalwater_depth, & ! bathymetry
1178                       config_flags%sf_surface_mosaic, config_flags%mosaic_cat, config_flags%num_land_cat, & ! Noah tiling
1179                       config_flags%maxpatch,           &   ! start of CLM variables
1180               grid%numc,grid%nump,grid%snl,grid%snowdp,&   !
1181               grid%wtc,grid%wtp,&
1182               grid%h2osno,grid%t_grnd,grid%t_veg,grid%h2ocan, &
1183               grid%h2ocan_col,grid%t2m_max,grid%t2m_min,&
1184               grid%t_ref2m,&
1185               grid%h2osoi_liq_s1,&
1186               grid%h2osoi_liq_s2,grid%h2osoi_liq_s3,&
1187               grid%h2osoi_liq_s4, &
1188               grid%h2osoi_liq_s5, &
1189               grid%h2osoi_liq1,grid%h2osoi_liq2,&
1190               grid%h2osoi_liq3,grid%h2osoi_liq4,&
1191               grid%h2osoi_liq5,grid%h2osoi_liq6, &
1192               grid%h2osoi_liq7,grid%h2osoi_liq8,&
1193               grid%h2osoi_liq9,   &
1194               grid%h2osoi_liq10, &
1195               grid%h2osoi_ice_s1,grid%h2osoi_ice_s2,&
1196               grid%h2osoi_ice_s3, &
1197               grid%h2osoi_ice_s4, &
1198               grid%h2osoi_ice_s5, &
1199               grid%h2osoi_ice1,&
1200               grid%h2osoi_ice2,&
1201               grid%h2osoi_ice3,grid%h2osoi_ice4,&
1202               grid%h2osoi_ice5, &
1203               grid%h2osoi_ice6, &
1204               grid%h2osoi_ice7,grid%h2osoi_ice8,&
1205               grid%h2osoi_ice9,grid%h2osoi_ice10,&
1206               grid%t_soisno_s1,grid%t_soisno_s2, &
1207               grid%t_soisno_s3,grid%t_soisno_s4, &
1208               grid%t_soisno_s5,grid%t_soisno1,&
1209               grid%t_soisno2,grid%t_soisno3,&
1210               grid%t_soisno4,grid%t_soisno5, &
1211               grid%t_soisno6,grid%t_soisno7,&
1212               grid%t_soisno8,grid%t_soisno9,&
1213               grid%t_soisno10,grid%dzsnow1,grid%dzsnow2,grid%dzsnow3,grid%dzsnow4,&
1214               grid%dzsnow5,grid%snowrds1,grid%snowrds2,grid%snowrds3,grid%snowrds4,grid%snowrds5, &
1215               grid%t_lake1,grid%t_lake2,&
1216               grid%t_lake3,grid%t_lake4, &
1217               grid%t_lake5,grid%t_lake6, &
1218               grid%t_lake7,grid%t_lake8, &
1219               grid%t_lake9,grid%t_lake10,&
1220               grid%h2osoi_vol1,grid%h2osoi_vol2, &
1221               grid%h2osoi_vol3,grid%h2osoi_vol4, &
1222               grid%h2osoi_vol5,&
1223               grid%h2osoi_vol6,&
1224               grid%h2osoi_vol7,grid%h2osoi_vol8,&
1225               grid%h2osoi_vol9,grid%h2osoi_vol10,&
1226               grid%ht,                          &
1227               grid%ALBEDOsubgrid,grid%LHsubgrid,&
1228               grid%HFXsubgrid,grid%LWUPsubgrid,&
1229               grid%Q2subgrid,grid%SABVsubgrid,  &
1230               grid%SABGsubgrid,grid%NRAsubgrid,&
1231               grid%SWUPsubgrid,grid%lhsoi, &
1232               grid%lhveg, grid%lhtran, & !end of CLM variables
1233                       grid%TSK_SAVE,                        & !Optional fractional seaice
1234                       grid%itimestep, grid%fdob,            &
1235                       t00, p00, a,                      & ! for obs_nudge base state
1236                       grid%TYR, grid%TYRA, grid%TDLY, grid%TLAG, grid%NYEAR, grid%NDAY,grid%tmn_update, &
1237                       grid%achfx, grid%aclhf, grid%acgrdflx,                  &
1238                       config_flags%nssl_cccn,                                 &
1239                       config_flags%nssl_alphah, config_flags%nssl_alphahl,    &
1240                       config_flags%nssl_cnoh, config_flags%nssl_cnohl,        &
1241                       config_flags%nssl_cnor, config_flags%nssl_cnos,         &
1242                       config_flags%nssl_rho_qh, config_flags%nssl_rho_qhl,    &
1243                       config_flags%nssl_rho_qs,                               &
1244                       config_flags%nssl_ipelec,                               &
1245                       config_flags%nssl_isaund                               &
1246                       ,grid%RQCNCUTEN, grid%RQINCUTEN,grid%rliq             &      !mchen add for cammpmg
1247                       ,grid%cldfra_dp,grid%cldfra_sh                        & ! ckay for subgrid cloud
1248                       ,grid%te_temf,grid%cf3d_temf,grid%wm_temf        & ! WA
1249                       ,grid%massflux_EDKF, grid%entr_EDKF, grid%detr_EDKF                &
1250                       ,grid%thl_up,grid%thv_up,grid%rt_up                                &
1251                       ,grid%rv_up,grid%rc_up,grid%u_up,grid%v_up,grid%frac_up            &
1252                       ,grid%RDCASHTEN, grid%RQCDCSHTEN                                   &
1253                       ,grid%cldareaa, grid%cldareab, grid%cldliqa, grid%cldliqb          &
1254                       ,grid%ca_rad, grid%cw_rad                                          &
1255                       ,grid%pblmax, grid%wub, grid%ltopb, grid%clddpthb, grid%cldtopb    &
1256                       ,grid%capesave, grid%ainckfsa, grid%radsave                        &
1257                       ,grid%rainsh, grid%rainshvb, grid%kdcldtop, grid%kdcldbas          &
1258                       ,grid%xtime1, grid%PBLHAVG, grid%TKEAVG                            &
1259                       ,grid%ccn_conc                                          & ! RAS
1260                       ,grid%QKE                                               &!JOE-for mynn
1261                       ,grid%pek_pbl,grid%pep_pbl                              & ! EEPS
1262                       ,grid%ghi_accum                                         & ! Solar diagnostics
1263                       ,grid%landusef,grid%landusef2,grid%mosaic_cat_index                            & ! danli mosaic
1264                       ,grid%TSK_mosaic,grid%TSLB_mosaic,grid%SMOIS_mosaic,grid%SH2O_mosaic           & ! danli mosaic
1265                       ,grid%CANWAT_mosaic,grid%SNOW_mosaic,grid%SNOWH_mosaic,grid%SNOWC_mosaic       & ! danli mosaic
1266                       ,grid%ALBEDO_mosaic,grid%ALBBCK_mosaic, grid%EMISS_mosaic                      & ! danli mosaic
1267                       ,grid%EMBCK_mosaic, grid%ZNT_mosaic, grid%Z0_mosaic                            & ! danli mosaic
1268                       ,grid%TR_URB2D_mosaic,grid%TB_URB2D_mosaic                                     & ! danli mosaic
1269                       ,grid%TG_URB2D_mosaic,grid%TC_URB2D_mosaic                                     & ! danli mosaic
1270                       ,grid%QC_URB2D_mosaic                                                          & ! danli mosaic
1271                       ,grid%TRL_URB3D_mosaic,grid%TBL_URB3D_mosaic                                   & ! danli mosaic
1272                       ,grid%TGL_URB3D_mosaic                                                         & ! danli mosaic
1273                       ,grid%SH_URB2D_mosaic,grid%LH_URB2D_mosaic                                     & ! danli mosaic
1274                       ,grid%G_URB2D_mosaic,grid%RN_URB2D_mosaic                                      & ! danli mosaic
1275                       ,grid%TS_URB2D_mosaic                                                          & ! danli mosaic
1276                       ,grid%TS_RUL2D_mosaic                                                          & ! danli mosaic
1277                       ,grid%irr_rand_field,config_flags%irr_ph,config_flags%irr_freq                 &
1278                       ,grid%QR_CU, grid%QS_CU                                                        &
1279                       ,grid%NC_CU, grid%NI_CU, grid%NR_CU, grid%NS_CU,grid%CCN_CU                    &
1280                       ,grid%alevsiz_cu,grid%num_months,grid%no_src_types_cu,grid%aeromcu             &
1281                       ,grid%aeropcu,grid%EFCG,grid%EFCS,grid%EFIG,grid%EFIS,grid%EFSG,grid%EFSS      &
1282                       )
1283        ENDDO ! loop of tiles for phy_init
1284
1285        if(.not.grid%restart) then
1286
1287          if(config_flags%sf_surface_physics .eq. NOAHMPSCHEME .and. config_flags%opt_run .eq. 5) then
1288               call groundwater_init ( grid ,                                                            &
1289                  config_flags%num_soil_layers, grid%dzs, grid%isltyp, grid%ivgtyp, config_flags%wtddt , &
1290                  grid%fdepthxy   , grid%ht     , grid%riverbedxy, grid%eqzwt     ,                      &
1291                  grid%rivercondxy, grid%pexpxy , grid%areaxy    , grid%zwtxy     ,                      &
1292                  grid%smois      , grid%sh2o   , grid%smoiseq   , grid%smcwtdxy  ,                      &
1293                  grid%qlatxy     , grid%qslatxy, grid%qrfxy     , grid%qrfsxy    ,                      &
1294                  grid%deeprechxy , grid%rechxy , grid%qspringxy , grid%qspringsxy,                      &
1295                  grid%rechclim   ,                                                                      &
1296                  ids,ide, jds,jde, kds,kde,                                                             &
1297                  ims,ime, jms,jme, kms,kme,                                                             &
1298                  ips,ipe, jps,jpe, kps,kpe,                                                             &
1299                  its,ite, jts,jte, kts,kte )
1300          endif
1301
1302        endif
1303
1304     ENDIF   ! no phy_init for the backwards part of the DFI
1305
1306    CALL wrf_debug ( 100 , 'start_domain_em: After call to phy_init' )
1307
1308    IF (config_flags%do_avgflx_em .EQ. 1) THEN
1309       WRITE ( message , FMT = '("start_em: initializing avgflx on domain ",I3)' ) &
1310            & grid%id
1311       CALL wrf_message(trim(message))
1312       grid%avgflx_count = 0
1313       f_flux = config_flags%do_avgflx_cugd .EQ. 1  ! WA 9/24/10
1314       DO ij = 1, grid%num_tiles
1315          call wrf_debug(200,'In start_em, before zero_avgflx call')
1316          if (.not. grid%restart) call zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
1317               &   ids, ide, jds, jde, kds, kde,           &
1318               &   ims, ime, jms, jme, kms, kme,           &
1319               &   grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, f_flux, &
1320               &   grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
1321               &   grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
1322          call wrf_debug(200,'In start_em, after zero_avgflx call')
1323       ENDDO
1324    ENDIF
1325
1326    call wrf_debug(100,'start_em: calling lightning_init')
1327    CALL lightning_init ( id=grid%id, itimestep=grid%itimestep, restart=grid%restart, dt=grid%dt, dx=grid%dx   &
1328                     ! Namelist control options
1329                         ,cu_physics=config_flags%cu_physics,mp_physics=config_flags%mp_physics    &
1330                         ,do_radar_ref=config_flags%do_radar_ref                                   &
1331                         ,lightning_option=config_flags%lightning_option                           &
1332                         ,lightning_dt=config_flags%lightning_dt                                   &
1333                         ,lightning_start_seconds=config_flags%lightning_start_seconds             &
1334                         ,ltngacttime=grid%ltngacttime                                             &
1335                         ,iccg_prescribed_num=config_flags%iccg_prescribed_num                     &
1336                         ,iccg_prescribed_den=config_flags%iccg_prescribed_den                     &
1337                         ,cellcount_method=config_flags%cellcount_method                           &
1338                     ! Order dependent args for domain, mem, and tile dims
1339                         ,ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde             &
1340                         ,ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme             &
1341                         ,its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte             &
1342                     ! IC and CG flash rates and accumulated flash count
1343                         ,ic_flashcount=grid%ic_flashcount, ic_flashrate=grid%ic_flashrate        &
1344                         ,cg_flashcount=grid%cg_flashcount, cg_flashrate=grid%cg_flashrate        &
1345 #if (WRF_CHEM == 1)
1346                         ,lnox_opt=config_flags%lnox_opt           &
1347                         ,lnox_passive=config_flags%lnox_passive   &
1348                         ,lnox_total=tracer(:,:,:,p_lnox_total)    &
1349                         ,lnox_ic=tracer(:,:,:,p_lnox_ic)          &
1350                         ,lnox_cg=tracer(:,:,:,p_lnox_cg)          &
1351 #endif
1352                     )
1353
1354     call wrf_debug(100,'start_em: after calling lightning_init')
1355
1356    END IF ! restart
1357
1358 #if 0
1359 #include "CYCLE_TEST.inc"
1360 #endif
1361
1362 !
1363 !
1364
1365  !  set physical boundary conditions for all initialized variables
1366
1367 !-----------------------------------------------------------------------
1368 !  Stencils for patch communications  (WCS, 29 June 2001)
1369 !  Note:  the size of this halo exchange reflects the
1370 !         fact that we are carrying the uncoupled variables
1371 !         as state variables in the mass coordinate model, as
1372 !         opposed to the coupled variables as in the height
1373 !         coordinate model.
1374 !
1375 !                           * * * * *
1376 !         *        * * *    * * * * *
1377 !       * + *      * + *    * * + * *
1378 !         *        * * *    * * * * *
1379 !                           * * * * *
1380 !
1381 !j  grid%u_1                          x
1382 !j  grid%u_2                          x
1383 !j  grid%v_1                          x
1384 !j  grid%v_2                          x
1385 !j  grid%w_1                          x
1386 !j  grid%w_2                          x
1387 !j  grid%t_1                          x
1388 !j  grid%t_2                          x
1389 !j grid%ph_1                          x
1390 !j grid%ph_2                          x
1391 !
1392 !j  grid%t_init                       x
1393 !
1394 !j  grid%phb   x
1395 !j  grid%ph0   x
1396 !j  grid%php   x
1397 !j   grid%pb   x
1398 !j   grid%al   x
1399 !j  grid%alt   x
1400 !j  grid%alb   x
1401 !
1402 !  the following are 2D (xy) variables
1403 !
1404 !j grid%mu_1                          x
1405 !j grid%mu_2                          x
1406 !j grid%mub   x
1407 !j grid%mu0   x
1408 !j grid%ht    x
1409 !j grid%msftx x
1410 !j grid%msfty x
1411 !j grid%msfux x
1412 !j grid%msfuy x
1413 !j grid%msfvx x
1414 !j grid%msfvy x
1415 !j grid%sina  x
1416 !j grid%cosa  x
1417 !j grid%e     x
1418 !j grid%f     x
1419 !
1420 !  4D variables
1421 !
1422 ! moist                        x
1423 !  chem                        x
1424 !scalar                        x
1425
1426 !--------------------------------------------------------------
1427
1428 #ifdef DM_PARALLEL
1429 # include "HALO_EM_INIT_1.inc"
1430 # include "HALO_EM_INIT_2.inc"
1431 # include "HALO_EM_INIT_3.inc"
1432 # include "HALO_EM_INIT_4.inc"
1433 # include "HALO_EM_INIT_5.inc"
1434       IF ( config_flags%sf_ocean_physics .EQ. PWP3DSCHEME ) THEN
1435 # include "HALO_EM_INIT_6.inc"
1436       END IF
1437 # include "PERIOD_BDY_EM_INIT.inc"
1438 # include "PERIOD_BDY_EM_MOIST.inc"
1439 # include "PERIOD_BDY_EM_TKE.inc"
1440 # include "PERIOD_BDY_EM_SCALAR.inc"
1441 # include "PERIOD_BDY_EM_CHEM.inc"
1442 if(config_flags%sf_surface_physics.eq.NOAHMPSCHEME.and.config_flags%opt_run.eq.5)then
1443 # include "HALO_EM_HYDRO_NOAHMP_INIT.inc"
1444 endif
1445 #endif
1446
1447
1448    CALL set_physical_bc3d( grid%u_1 , 'U' , config_flags ,                  &
1449                          ids , ide , jds , jde , kds , kde ,        &
1450                          ims , ime , jms , jme , kms , kme ,        &
1451                          its , ite , jts , jte , kts , kte ,        &
1452                          its , ite , jts , jte , kts , kte )
1453    CALL set_physical_bc3d( grid%u_2 , 'U' , config_flags ,                  &
1454                          ids , ide , jds , jde , kds , kde ,        &
1455                          ims , ime , jms , jme , kms , kme ,        &
1456                          its , ite , jts , jte , kts , kte ,        &
1457                          its , ite , jts , jte , kts , kte )
1458
1459    CALL set_physical_bc3d( grid%v_1 , 'V' , config_flags ,                  &
1460                          ids , ide , jds , jde , kds , kde ,        &
1461                          ims , ime , jms , jme , kms , kme ,        &
1462                          its , ite , jts , jte , kts , kte ,        &
1463                          its , ite , jts , jte , kts , kte )
1464    CALL set_physical_bc3d( grid%v_2 , 'V' , config_flags ,                  &
1465                          ids , ide , jds , jde , kds , kde ,        &
1466                          ims , ime , jms , jme , kms , kme ,        &
1467                          its , ite , jts , jte , kts , kte ,        &
1468                          its , ite , jts , jte , kts , kte )
1469
1470 ! set kinematic condition for w
1471
1472    CALL set_physical_bc2d( grid%ht , 'r' , config_flags ,                &
1473                            ids , ide , jds , jde , &
1474                            ims , ime , jms , jme , &
1475                            its , ite , jts , jte , &
1476                            its , ite , jts , jte   )
1477
1478 ! Set the w at the surface.  If this is the start of a forecast, or if this is a cycled run
1479 ! and the start of that forecast, we define the w field.  However, if this a restart, then
1480 ! we leave w alone.  Initial value is V dot grad(topo) at surface, then smoothly decaying
1481 ! above that.
1482
1483    IF ( ( start_of_simulation .OR. config_flags%cycling ) .AND. ( .NOT. config_flags%restart ) ) THEN
1484
1485 ! If W already exists (not zero), then we leave it alone.  How to do this?  We find the
1486 ! max/min on this node at the surface.  If parallel, we collect the max/min from all procs.
1487 ! If the max/min throughout the entire domain at the surface is identically 0, then we say
1488 ! that the W field is NOT initialized, and we run the set_w_surface routines for the
1489 ! two time levels of W.  If the field is already initialized, we do NOT run those two
1490 ! routines.
1491
1492       w_max = grid%w_2(its,1,jts)
1493       w_min = grid%w_2(its,1,jts)
1494       DO j = jts, MIN(jte,jde-1)
1495          DO i = its, MIN(ite,ide-1)
1496          w_max = MAX ( w_max , grid%w_2(i,1,j) )
1497          w_min = MIN ( w_min , grid%w_2(i,1,j) )
1498          END DO
1499       END DO
1500 #ifdef DM_PARALLEL
1501       w_max = wrf_dm_max_real ( w_max )
1502       w_min = wrf_dm_min_real ( w_min )
1503 #endif
1504
1505       IF ( ( ABS(w_max) .LT. 1.E-6 ) .AND. &
1506            ( ABS(w_min) .LT. 1.E-6 ) ) THEN
1507          w_needs_to_be_set = .TRUE.
1508          ! Force use_input_w to false, if there is no input w
1509 #if ( WRFPLUS == 1 )
1510          IF ( config_flags%use_input_w ) THEN
1511             config_flags%use_input_w = .FALSE.
1512             CALL nl_set_use_input_w (grid%id, .FALSE.)
1513          ENDIF
1514 #endif
1515       ELSE
1516          IF ( config_flags%use_input_w ) THEN
1517             w_needs_to_be_set = .FALSE.
1518          ELSE
1519             w_needs_to_be_set = .TRUE.
1520          END IF
1521       END IF
1522
1523       IF ( w_needs_to_be_set ) THEN
1524
1525          ! setting kinematic condition for w at the surface
1526
1527          fill_w_flag = .true.
1528          CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
1529                               grid%w_1, grid%ht, grid%u_1, grid%v_1, grid%cf1, &
1530                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
1531                               ids, ide, jds, jde, kds, kde,                    &
1532                               ims, ime, jms, jme, kms, kme,                    &
1533                               its, ite, jts, jte, kts, kte                     )
1534          CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
1535                               grid%w_2, grid%ht, grid%u_2, grid%v_2, grid%cf1, &
1536                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
1537                               ids, ide, jds, jde, kds, kde,                    &
1538                               ims, ime, jms, jme, kms, kme,                    &
1539                               its, ite, jts, jte, kts, kte                     )
1540       END IF
1541    END IF
1542
1543    IF ( .NOT. config_flags%restart ) THEN
1544 ! set up slope-radiation constant arrays based on topography
1545 ! paj: also computes laplacian for topo_wind
1546 !
1547        DO j = jts,min(jte,jde-1)
1548        DO i = its, min(ite,ide-1)
1549           if (config_flags%periodic_x ) then
1550             im1 = i-1
1551             ip1 = i+1
1552           else
1553             im1 = max(i-1,ids)
1554             ip1 = min(i+1,ide-1)
1555           endif
1556
1557           if (config_flags%periodic_y ) then
1558             jm1 = j-1
1559             jp1 = j+1
1560           else
1561             jm1 = max(j-1,jds)
1562             jp1 = min(j+1,jde-1)
1563           endif
1564           grid%toposlpx(i,j)=(grid%ht(ip1,j)-grid%ht(im1,j))*grid%msftx(i,j)*grid%rdx/(ip1-im1)
1565           grid%toposlpy(i,j)=(grid%ht(i,jp1)-grid%ht(i,jm1))*grid%msfty(i,j)*grid%rdy/(jp1-jm1)
1566           grid%lap_hgt(i,j)=(grid%ht(ip1,j)+grid%ht(im1,j)+grid%ht(i,jp1)+grid%ht(i,jm1)-grid%ht(i,j)*4.)/4.
1567              hx = grid%toposlpx(i,j)
1568              hy = grid%toposlpy(i,j)
1569              pi = 4.*atan(1.)
1570              grid%slope(i,j) = atan((hx**2+hy**2)**.5)
1571              if (grid%slope(i,j).lt.1.e-4) then
1572                grid%slope(i,j) = 0.
1573                grid%slp_azi(i,j) = 0.
1574              else
1575                grid%slp_azi(i,j) = atan2(hx,hy)+pi
1576
1577 ! Rotate slope azimuth to lat-lon grid
1578                if (grid%cosa(i,j).ge.0) then
1579                  grid%slp_azi(i,j) = grid%slp_azi(i,j) - asin(grid%sina(i,j))
1580                else
1581                  grid%slp_azi(i,j) = grid%slp_azi(i,j) - (pi - asin(grid%sina(i,j)))
1582                endif
1583              endif
1584        ENDDO
1585        ENDDO
1586 !
1587 ! paj: computes ct and ct2 for topo_wind
1588 !
1589        grid%ctopo=1.
1590        grid%ctopo2=1.
1591 !
1592        if (config_flags%topo_wind.eq.1) then
1593          DO j = jts,min(jte,jde-1)
1594          DO i = its, min(ite,ide-1)
1595             if(grid%xland(i,j).lt.1.5)grid%ctopo(i,j)=sqrt(grid%var_sso(i,j))
1596             if (grid%ctopo(i,j).le.2.718) then
1597             grid%ctopo(i,j)=1.
1598             else
1599             grid%ctopo(i,j)=alog(grid%ctopo(i,j))
1600             endif
1601 !
1602             if (grid%lap_hgt(i,j).gt.-10.) then
1603               grid%ctopo(i,j)=grid%ctopo(i,j)
1604             else
1605                if (grid%lap_hgt(i,j).ge.-20) then
1606                  alpha=(grid%lap_hgt(i,j)+20.)/10.
1607                  grid%ctopo(i,j)=alpha*grid%ctopo(i,j)+(1-alpha)
1608                else
1609                   if (grid%lap_hgt(i,j).ge.-30.) then
1610                   grid%ctopo(i,j)=(grid%lap_hgt(i,j)+30.)/10.
1611                   grid%ctopo2(i,j)=(grid%lap_hgt(i,j)+30.)/10.
1612                   else
1613                   grid%ctopo(i,j)=0.
1614                   grid%ctopo2(i,j)=0.
1615                   endif
1616                endif
1617             endif
1618          ENDDO
1619          ENDDO
1620 !       if topo_wind==2 (D.Ovens, C.Mass)
1621        else if (config_flags%topo_wind.eq.2) then
1622          DO j = jts,min(jte,jde-1)
1623          DO i = its, min(ite,ide-1)
1624            if (grid%xland(i,j).lt.1.5) then
1625               vfac = amin1(1.575,(grid%var2d(i,j)*0.4/200.+1.175))
1626               vfac = vfac * vfac
1627            else                      ! over water, leave it alone
1628               vfac = 1.
1629            endif
1630 !          print *, j,i,grid%ctopo(i,j),vfac
1631            grid%ctopo(i,j)=grid%ctopo(i,j)*vfac
1632          ENDDO
1633          ENDDO
1634        endif
1635
1636    END IF
1637
1638    CALL set_physical_bc3d( grid%w_1 , 'W' , config_flags ,                 &
1639                          ids , ide , jds , jde , kds , kde ,        &
1640                          ims , ime , jms , jme , kms , kme ,        &
1641                          its , ite , jts , jte , kts , kte ,        &
1642                          its , ite , jts , jte , kts , kte )
1643    CALL set_physical_bc3d( grid%w_2 , 'W' , config_flags ,                 &
1644                          ids , ide , jds , jde , kds , kde ,        &
1645                          ims , ime , jms , jme , kms , kme ,        &
1646                          its , ite , jts , jte , kts , kte ,        &
1647                          its , ite , jts , jte , kts , kte )
1648
1649    CALL set_physical_bc3d( grid%ph_1 , 'W' , config_flags ,                 &
1650                          ids , ide , jds , jde , kds , kde ,        &
1651                          ims , ime , jms , jme , kms , kme ,        &
1652                          its , ite , jts , jte , kts , kte ,        &
1653                          its , ite , jts , jte , kts , kte )
1654
1655    CALL set_physical_bc3d( grid%ph_2 , 'W' , config_flags ,                 &
1656                          ids , ide , jds , jde , kds , kde ,        &
1657                          ims , ime , jms , jme , kms , kme ,        &
1658                          its , ite , jts , jte , kts , kte ,        &
1659                          its , ite , jts , jte , kts , kte )
1660
1661    CALL set_physical_bc3d( grid%t_1 , 't' , config_flags ,                 &
1662                          ids , ide , jds , jde , kds , kde ,        &
1663                          ims , ime , jms , jme , kms , kme ,        &
1664                          its , ite , jts , jte , kts , kte ,        &
1665                          its , ite , jts , jte , kts , kte )
1666
1667    CALL set_physical_bc3d( grid%t_2 , 't' , config_flags ,                 &
1668                          ids , ide , jds , jde , kds , kde ,        &
1669                          ims , ime , jms , jme , kms , kme ,        &
1670                          its , ite , jts , jte , kts , kte ,        &
1671                          its , ite , jts , jte , kts , kte )
1672
1673    CALL set_physical_bc2d( grid%mu_1, 't' , config_flags ,   &
1674                            ids , ide , jds , jde ,  &
1675                            ims , ime , jms , jme ,  &
1676                            its , ite , jts , jte ,  &
1677                            its , ite , jts , jte   )
1678    CALL set_physical_bc2d( grid%mu_2, 't' , config_flags ,   &
1679                            ids , ide , jds , jde ,  &
1680                            ims , ime , jms , jme ,  &
1681                            its , ite , jts , jte ,  &
1682                            its , ite , jts , jte   )
1683    CALL set_physical_bc2d( grid%mub , 't' , config_flags ,   &
1684                            ids , ide , jds , jde ,  &
1685                            ims , ime , jms , jme ,  &
1686                            its , ite , jts , jte ,  &
1687                            its , ite , jts , jte   )
1688 !  CALL set_physical_bc2d( grid%mu0 , 't' , config_flags ,   &
1689 !                          ids , ide , jds , jde ,  &
1690 !                          ims , ime , jms , jme ,  &
1691 !                          its , ite , jts , jte ,  &
1692 !                          its , ite , jts , jte   )
1693
1694
1695    CALL set_physical_bc3d( grid%phb , 'W' , config_flags ,                 &
1696                          ids , ide , jds , jde , kds , kde ,        &
1697                          ims , ime , jms , jme , kms , kme ,        &
1698                          its , ite , jts , jte , kts , kte ,        &
1699                          its , ite , jts , jte , kts , kte )
1700    CALL set_physical_bc3d( grid%ph0 , 'W' , config_flags ,                 &
1701                          ids , ide , jds , jde , kds , kde ,        &
1702                          ims , ime , jms , jme , kms , kme ,        &
1703                          its , ite , jts , jte , kts , kte ,        &
1704                          its , ite , jts , jte , kts , kte )
1705    CALL set_physical_bc3d( grid%php , 'W' , config_flags ,                 &
1706                          ids , ide , jds , jde , kds , kde ,        &
1707                          ims , ime , jms , jme , kms , kme ,        &
1708                          its , ite , jts , jte , kts , kte ,        &
1709                          its , ite , jts , jte , kts , kte )
1710
1711    CALL set_physical_bc3d( grid%pb , 't' , config_flags ,                 &
1712                          ids , ide , jds , jde , kds , kde ,        &
1713                          ims , ime , jms , jme , kms , kme ,        &
1714                          its , ite , jts , jte , kts , kte ,        &
1715                          its , ite , jts , jte , kts , kte )
1716    CALL set_physical_bc3d( grid%al , 't' , config_flags ,                 &
1717                          ids , ide , jds , jde , kds , kde ,        &
1718                          ims , ime , jms , jme , kms , kme ,        &
1719                          its , ite , jts , jte , kts , kte ,        &
1720                          its , ite , jts , jte , kts , kte )
1721    CALL set_physical_bc3d( grid%alt , 't' , config_flags ,                 &
1722                          ids , ide , jds , jde , kds , kde ,        &
1723                          ims , ime , jms , jme , kms , kme ,        &
1724                          its , ite , jts , jte , kts , kte ,        &
1725                          its , ite , jts , jte , kts , kte )
1726    CALL set_physical_bc3d( grid%alb , 't' , config_flags ,                 &
1727                          ids , ide , jds , jde , kds , kde ,        &
1728                          ims , ime , jms , jme , kms , kme ,        &
1729                          its , ite , jts , jte , kts , kte ,        &
1730                          its , ite , jts , jte , kts , kte )
1731    CALL set_physical_bc3d(grid%t_init, 't' , config_flags ,                 &
1732                          ids , ide , jds , jde , kds , kde ,        &
1733                          ims , ime , jms , jme , kms , kme ,        &
1734                          its , ite , jts , jte , kts , kte ,        &
1735                          its , ite , jts , jte , kts , kte )
1736    CALL set_physical_bc3d(grid%tke_2, 't' , config_flags ,                 &
1737                          ids , ide , jds , jde , kds , kde ,        &
1738                          ims , ime , jms , jme , kms , kme ,        &
1739                          its , ite , jts , jte , kts , kte ,        &
1740                          its , ite , jts , jte , kts , kte )
1741
1742    IF (num_moist > 0) THEN
1743
1744 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1745
1746       loop_3d_m   : DO loop = 1 , num_moist
1747          CALL set_physical_bc3d( moist(:,:,:,loop) , 'r' , config_flags ,                 &
1748                                  ids , ide , jds , jde , kds , kde ,        &
1749                                  ims , ime , jms , jme , kms , kme ,        &
1750                                  its , ite , jts , jte , kts , kte ,        &
1751                                  its , ite , jts , jte , kts , kte )
1752       END DO loop_3d_m
1753
1754    ENDIF
1755
1756
1757 ! Some initializations for the simple ccn field.
1758
1759    IF ( f_qnn ) THEN
1760       IF ( config_flags%mp_physics == wdm5scheme .or. config_flags%mp_physics == wdm6scheme ) THEN
1761          ! NO OP
1762       ELSE IF ( config_flags%mp_physics == nssl_2momccn ) THEN
1763          grid%ccn_conc =  config_flags%nssl_cccn/1.225
1764       ELSE
1765          ! NO OP
1766       END IF
1767       ccn_max_val = MAXVAL(scalar(its:MIN(ite,ide-1),kts:kte-1,jts:MIN(jte,jde-1),p_qnn))
1768 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1769       ccn_max_val = wrf_dm_max_real ( ccn_max_val )
1770 #endif
1771       IF ( ccn_max_val < 1.0 ) THEN ! initialization of ccn not already done
1772          DO j=jts,MIN(jte,jde-1)
1773             DO k=kts,kte
1774                DO i=its,MIN(ite,ide-1)
1775                   scalar(i,k,j,p_qnn) = grid%ccn_conc
1776                END DO
1777             END DO
1778          END DO
1779       END IF
1780    END IF
1781
1782    IF (num_scalar > 0) THEN
1783
1784 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1785
1786       loop_3d_s   : DO loop = 1 , num_scalar
1787          CALL set_physical_bc3d( scalar(:,:,:,loop) , 'r' , config_flags ,                 &
1788                                  ids , ide , jds , jde , kds , kde ,        &
1789                                  ims , ime , jms , jme , kms , kme ,        &
1790                                  its , ite , jts , jte , kts , kte ,        &
1791                                  its , ite , jts , jte , kts , kte )
1792       END DO loop_3d_s
1793
1794    ENDIF
1795
1796
1797 #if (WRF_CHEM == 1)
1798         if(config_flags%tracer_opt > 0 )then
1799            call initialize_tracer (tracer,config_flags%chem_in_opt, &
1800                                config_flags%tracer_opt,num_tracer,  &
1801                                ids,ide, jds,jde, kds,kde,      & ! domain dims
1802                                ims,ime, jms,jme, kms,kme,      & ! memory dims
1803                                ips,ipe, jps,jpe, kps,kpe,      & ! patch  dims
1804                                its,ite, jts,jte, kts,kte )
1805
1806         endif
1807 !
1808 !   we do this here, so we only have one chem_init routine for either core....
1809 !
1810          do j=jts,min(jte,jde-1)
1811             do i=its,min(ite,ide-1)
1812                do k=kts,kte
1813                  z_at_w(i,k,j)=(grid%ph_2(i,k,j)+grid%phb(i,k,j))/g
1814                enddo
1815                do k=kts,min(kte,kde-1)
1816                  tempfac=(grid%t_1(i,k,j) + t0)*((grid%p(i,k,j) + grid%pb(i,k,j))/p1000mb)**rcp
1817                  convfac(i,k,j) = (grid%p(i,k,j)+grid%pb(i,k,j))/rgasuniv/tempfac
1818                enddo
1819             enddo
1820          enddo
1821
1822          CALL chem_init (grid%id,chem,emis_ant,scalar,grid%dt,grid%bioemdt,grid%photdt,     &
1823                          grid%chemdt,                                       &
1824                          grid%stepbioe,grid%stepphot,grid%stepchem,grid%stepfirepl,      &
1825                          grid%plumerisefire_frq,z_at_w,grid%xlat,grid%xlong,g,           &
1826                          grid%aerwrf,config_flags,grid,            &
1827                          grid%alt,grid%t_1,grid%p,convfac,grid%ttday,grid%tcosz,         &
1828                          grid%julday,grid%gmt,&
1829                          grid%tauaer1,grid%tauaer2,grid%tauaer3,grid%tauaer4,                   &
1830                          grid%gaer1,grid%gaer2,grid%gaer3,grid%gaer4,                           &
1831                          grid%waer1,grid%waer2,grid%waer3,grid%waer4,                           &
1832                          grid%l2aer,grid%l3aer,grid%l4aer,grid%l5aer,grid%l6aer,grid%l7aer,     &
1833                          grid%extaerlw1,grid%extaerlw2,grid%extaerlw3,grid%extaerlw4,           &
1834                          grid%extaerlw5,grid%extaerlw6,grid%extaerlw7,grid%extaerlw8,           &
1835                          grid%extaerlw9,grid%extaerlw10,grid%extaerlw11,grid%extaerlw12,        &
1836                          grid%extaerlw13,grid%extaerlw14,grid%extaerlw15,grid%extaerlw16,       &
1837                          grid%tauaerlw1,grid%tauaerlw2,grid%tauaerlw3,grid%tauaerlw4,           &
1838                          grid%tauaerlw5,grid%tauaerlw6,grid%tauaerlw7,grid%tauaerlw8,           &
1839                          grid%tauaerlw9,grid%tauaerlw10,grid%tauaerlw11,grid%tauaerlw12,        &
1840                          grid%tauaerlw13,grid%tauaerlw14,grid%tauaerlw15,grid%tauaerlw16,       &
1841                          grid%dgnum4d, grid%dgnumwet4d, grid%dgnum_a1, grid%dgnum_a2,           &
1842                          grid%dgnum_a3, grid%dgnumwet_a1, grid%dgnumwet_a2, grid%dgnumwet_a3,   &
1843                          grid%pm2_5_dry,grid%pm2_5_water,grid%pm2_5_dry_ec,                &
1844                          grid%tsoa,grid%asoa,grid%bsoa,                                    &
1845                          grid%last_chem_time_year,grid%last_chem_time_month,               &
1846                          grid%last_chem_time_day,grid%last_chem_time_hour,                 &
1847                          grid%last_chem_time_minute,grid%last_chem_time_second,            &
1848                          grid%chem_in_opt,grid%kemit,grid%num_vert_mix,numgas,             &
1849                          ids , ide , jds , jde , kds , kde ,                &
1850                          ims , ime , jms , jme , kms , kme ,                &
1851                          its , ite , jts , jte , kts , kte                  )
1852 #if( WRF_USE_CLM == 1)
1853 !-------------------------------------------------------------------
1854 !   check that the megan emission species are in chem option
1855 !-------------------------------------------------------------------
1856          if (config_flags%chem_opt > 0) then
1857            write(message,'('' start_domain_em: numgas = '',i4)') numgas
1858            CALL wrf_message( trim(message) )
1859            CALL wrf_spc_chk( numgas )
1860          endif
1861 #endif
1862
1863 !
1864 ! calculate initial pm
1865 !
1866 !       print *,'calculating initial pm'
1867         select case (config_flags%chem_opt)
1868         case (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2)
1869            call sum_pm_gocart (                                             &
1870                 grid%alt, chem, grid%pm2_5_dry, grid%pm2_5_dry_ec,grid%pm10,&
1871                 ids,ide, jds,jde, kds,kde,                                  &
1872                 ims,ime, jms,jme, kms,kme,                                  &
1873                 its,ite, jts,jte, kts,kte-1                                 )
1874         case (RADM2SORG,RADM2SORG_AQ,RADM2SORG_AQCHEM,RADM2SORG_KPP,RACMSORG_AQ,RACMSORG_AQCHEM_KPP, &
1875               RACM_ESRLSORG_AQCHEM_KPP,RACMSORG_KPP)
1876            call sum_pm_sorgam (                                             &
1877                 grid%alt, chem, grid%h2oaj, grid%h2oai,                                    &
1878                 grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1879                 grid%dust_opt,ids,ide, jds,jde, kds,kde,                    &
1880                 ims,ime, jms,jme, kms,kme,                                  &
1881                 its,ite, jts,jte, kts,kte-1                                 )
1882          case (RACM_SOA_VBS_KPP,RACM_SOA_VBS_AQCHEM_KPP)
1883              call sum_pm_soa_vbs (                                           &
1884                  grid%alt, chem, grid%h2oaj, grid%h2oai,                                   &
1885                  grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                &
1886                  config_flags%dust_opt,ids,ide, jds,jde, kds,kde,           &
1887                  ims,ime, jms,jme, kms,kme,                                 &
1888                  its,ite, jts,jte, kts,kte-1                                )
1889
1890         case (CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ, &
1891               CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, CBMZ_MOSAIC_DMS_4BIN_AQ,       &
1892               CBMZ_MOSAIC_DMS_8BIN_AQ, CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP)
1893            call sum_pm_mosaic (                                             &
1894                 grid%alt, chem,                                                  &
1895                 grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1896                 ids,ide, jds,jde, kds,kde,                                  &
1897                 ims,ime, jms,jme, kms,kme,                                  &
1898                 its,ite, jts,jte, kts,kte-1                                 )
1899
1900         case default
1901            do j=jts,min(jte,jde-1)
1902               do k=kts,min(kte,kde-1)
1903                  do i=its,min(ite,ide-1)
1904                     grid%pm2_5_dry(i,k,j)    = 0.
1905                     grid%pm2_5_water(i,k,j)  = 0.
1906                     grid%pm2_5_dry_ec(i,k,j) = 0.
1907                     grid%pm10(i,k,j)         = 0.
1908                     grid%tsoa(i,k,j)         = 0.
1909                  enddo
1910               enddo
1911            enddo
1912         end select
1913
1914         ! initialize advective tendency diagnostics for chem
1915         if ( grid%itimestep .eq. 0 .and. config_flags%chemdiag .eq. USECHEMDIAG ) then
1916           grid%advh_ct(:,:,:,:) = 0.
1917           grid%advz_ct(:,:,:,:) = 0.
1918         endif
1919 #endif
1920
1921         ! initialize advective tendency diagnostics for non-chem
1922         if ( grid%itimestep .eq. 0 .and. config_flags%tenddiag .eq. USETENDDIAG ) then
1923           advh_t(:,:,:,:) = 0.
1924           advz_t(:,:,:,:) = 0.
1925         endif
1926
1927    IF (num_chem >= PARAM_FIRST_SCALAR ) THEN
1928 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1929
1930       loop_3d_c   : DO loop = PARAM_FIRST_SCALAR , num_chem
1931          CALL set_physical_bc3d( chem(:,:,:,loop) , 'r' , config_flags ,                 &
1932                                  ids , ide , jds , jde , kds , kde ,        &
1933                                  ims , ime , jms , jme , kms , kme ,        &
1934                                  its , ite , jts , jte , kts , kte ,        &
1935                                  its , ite , jts , jte , kts , kte )
1936       END DO loop_3d_c
1937
1938    ENDIF
1939
1940    CALL set_physical_bc2d( grid%msftx , 'r' , config_flags ,              &
1941                          ids , ide , jds , jde , &
1942                          ims , ime , jms , jme , &
1943                          its , ite , jts , jte , &
1944                          its , ite , jts , jte   )
1945    CALL set_physical_bc2d( grid%msfty , 'r' , config_flags ,              &
1946                          ids , ide , jds , jde , &
1947                          ims , ime , jms , jme , &
1948                          its , ite , jts , jte , &
1949                          its , ite , jts , jte   )
1950    CALL set_physical_bc2d( grid%msfux , 'x' , config_flags ,              &
1951                          ids , ide , jds , jde , &
1952                          ims , ime , jms , jme , &
1953                          its , ite , jts , jte , &
1954                          its , ite , jts , jte   )
1955    CALL set_physical_bc2d( grid%msfuy , 'x' , config_flags ,              &
1956                          ids , ide , jds , jde , &
1957                          ims , ime , jms , jme , &
1958                          its , ite , jts , jte , &
1959                          its , ite , jts , jte   )
1960    CALL set_physical_bc2d( grid%msfvx , 'y' , config_flags ,              &
1961                          ids , ide , jds , jde , &
1962                          ims , ime , jms , jme , &
1963                          its , ite , jts , jte , &
1964                          its , ite , jts , jte   )
1965    CALL set_physical_bc2d( grid%msfvy , 'y' , config_flags ,              &
1966                          ids , ide , jds , jde , &
1967                          ims , ime , jms , jme , &
1968                          its , ite , jts , jte , &
1969                          its , ite , jts , jte   )
1970    CALL set_physical_bc2d( grid%sina , 'r' , config_flags ,              &
1971                          ids , ide , jds , jde , &
1972                          ims , ime , jms , jme , &
1973                          its , ite , jts , jte , &
1974                          its , ite , jts , jte   )
1975    CALL set_physical_bc2d( grid%cosa , 'r' , config_flags ,              &
1976                          ids , ide , jds , jde , &
1977                          ims , ime , jms , jme , &
1978                          its , ite , jts , jte , &
1979                          its , ite , jts , jte   )
1980    CALL set_physical_bc2d( grid%e , 'r' , config_flags ,                 &
1981                          ids , ide , jds , jde , &
1982                          ims , ime , jms , jme , &
1983                          its , ite , jts , jte , &
1984                          its , ite , jts , jte   )
1985    CALL set_physical_bc2d( grid%f , 'r' , config_flags ,                 &
1986                          ids , ide , jds , jde , &
1987                          ims , ime , jms , jme , &
1988                          its , ite , jts , jte , &
1989                          its , ite , jts , jte   )
1990
1991 #if (WRF_CHEM != 1)
1992       DEALLOCATE(CLDFRA_OLD)
1993 #endif
1994 #ifdef DM_PARALLEL
1995 # include "HALO_EM_INIT_1.inc"
1996 # include "HALO_EM_INIT_2.inc"
1997 # include "HALO_EM_INIT_3.inc"
1998 # include "HALO_EM_INIT_4.inc"
1999 # include "HALO_EM_INIT_5.inc"
2000 # include "PERIOD_BDY_EM_INIT.inc"
2001 # include "PERIOD_BDY_EM_MOIST.inc"
2002 # include "PERIOD_BDY_EM_TKE.inc"
2003 # include "PERIOD_BDY_EM_SCALAR.inc"
2004 # include "PERIOD_BDY_EM_CHEM.inc"
2005 #endif
2006
2007 DEALLOCATE(z_at_q)
2008
2009    IF (config_flags%p_lev_diags == PRESS_DIAGS ) THEN
2010     CALL wrf_debug ( 200 , ' PLD: pressure level diags' )
2011     CALL pld (                                  &
2012                  !  Input data for computing
2013                   U=grid%u_2                    &
2014                  ,V=grid%v_2                    &
2015                  ,W=grid%w_2                    &
2016                  ,t=grid%t_2                    &
2017                  ,qv=moist(:,:,:,P_QV)          &
2018                  ,zp=grid%ph_2                  &
2019                  ,zb=grid%phb                   &
2020                  ,pp=grid%p                     &
2021                  ,pb=grid%pb                    &
2022                  ,p=grid%p_hyd                  &
2023                  ,pw=grid%p_hyd_w               &
2024                  !  Map factors, coriolis for diags
2025                  ,msfux=grid%msfux              &
2026                  ,msfuy=grid%msfuy              &
2027                  ,msfvx=grid%msfvx              &
2028                  ,msfvy=grid%msfvy              &
2029                  ,msftx=grid%msftx              &
2030                  ,msfty=grid%msfty              &
2031                  ,f=grid%f                      &
2032                  ,e=grid%e                      &
2033                  !  Namelist info
2034                  ,use_tot_or_hyd_p=config_flags%use_tot_or_hyd_p &
2035                  ,extrap_below_grnd=config_flags%extrap_below_grnd &
2036                  ,missing=config_flags%p_lev_missing &
2037                  !  The diagnostics, mostly output variables
2038                  ,num_press_levels=config_flags%num_press_levels &
2039                  ,max_press_levels=max_plevs    &
2040                  ,press_levels=model_config_rec%press_levels &
2041                  ,p_pl  = grid%p_pl             &
2042                  ,u_pl  = grid%u_pl             &
2043                  ,v_pl  = grid%v_pl             &
2044                  ,t_pl  = grid%t_pl             &
2045                  ,rh_pl = grid%rh_pl            &
2046                  ,ght_pl= grid%ght_pl           &
2047                  ,s_pl  = grid%s_pl             &
2048                  ,td_pl = grid%td_pl            &
2049                  ,q_pl = grid%q_pl              &
2050                  !  Dimension arguments
2051                  ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde    &
2052                  ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme    &
2053                  ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte    )
2054    ENDIF
2055
2056    IF (config_flags%z_lev_diags == Z_DIAGS ) THEN
2057      CALL wrf_debug ( 200 , ' ZLD: height level and AGL diags' )
2058      CALL zld (                                  &
2059                  !  Input data for computing
2060                   U=grid%u_2                    &
2061                  ,V=grid%v_2                    &
2062                  ,W=grid%w_2                    &
2063                  ,t=grid%t_2                    &
2064                  ,qv=moist(:,:,:,P_QV)          &
2065                  ,zp=grid%ph_2                  &
2066                  ,zb=grid%phb                   &
2067                  ,pp=grid%p                     &
2068                  ,pb=grid%pb                    &
2069                  ,p=grid%p_hyd                  &
2070                  ,pw=grid%p_hyd_w               &
2071                  !  Map factors, coriolis for diags
2072                  ,msfux=grid%msfux              &
2073                  ,msfuy=grid%msfuy              &
2074                  ,msfvx=grid%msfvx              &
2075                  ,msfvy=grid%msfvy              &
2076                  ,msftx=grid%msftx              &
2077                  ,msfty=grid%msfty              &
2078                  ,f=grid%f                      &
2079                  ,e=grid%e                      &
2080                  ,ht=grid%ht                    &
2081                  !  Namelist info
2082                  ,use_tot_or_hyd_p=config_flags%use_tot_or_hyd_p &
2083                  ,extrap_below_grnd=config_flags%extrap_below_grnd &
2084                  ,missing=config_flags%z_lev_missing &
2085                  !  The diagnostics, mostly output variables
2086                  ,num_z_levels=config_flags%num_z_levels &
2087                  ,max_z_levels=max_zlevs    &
2088                  ,z_levels=model_config_rec%z_levels &
2089                  ,z_zl  = grid%z_zl             &
2090                  ,u_zl  = grid%u_zl             &
2091                  ,v_zl  = grid%v_zl             &
2092                  ,t_zl  = grid%t_zl             &
2093                  ,rh_zl = grid%rh_zl            &
2094                  ,ght_zl= grid%ght_zl           &
2095                  ,s_zl  = grid%s_zl             &
2096                  ,td_zl = grid%td_zl            &
2097                  ,q_zl = grid%q_zl              &
2098                  ,p_zl = grid%p_zl              &
2099                  !  Dimension arguments
2100                  ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde    &
2101                  ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme    &
2102                  ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte    )
2103    ENDIF
2104
2105    IF ( config_flags%diag_nwp2 == DO_TRAD_FIELDS ) THEN
2106       CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: TRAD_FIELDS' )
2107       DO j = jts, MIN(jte,jde-1)
2108          DO k = 1, kte-1
2109             DO i = its, MIN(ite,ide-1)
2110                dz8w(i,k,j) = (grid%phb(i,k+1,j)+grid%ph_2(i,k+1,j))/g - (grid%phb(i,k,j)+grid%ph_2(i,k,j))/g
2111                grid%rho(i,k,j) = 1./( grid%al(i,k,j)+grid%alb(i,k,j) )
2112             END DO
2113          END DO
2114       END DO
2115       CALL trad_fields ( &
2116                !  Input data for computing
2117                U=grid%u_2                                           &
2118                ,V=grid%v_2                                          &
2119                ,W=grid%w_2                                          &
2120                ,t=grid%th_phy_m_t0                                  &
2121                ,qv=moist(:,:,:,P_QV)                                &
2122                ,zp=grid%ph_2                                        &
2123                ,zb=grid%phb                                         &
2124                ,pp=grid%p                                           &
2125                ,pb=grid%pb                                          &
2126                ,p=grid%p_hyd                                        &
2127                ,pw=grid%p_hyd_w                                     &
2128                !  Map factors, coriolis for diags
2129                ,msfux=grid%msfux                                    &
2130                ,msfuy=grid%msfuy                                    &
2131                ,msfvx=grid%msfvx                                    &
2132                ,msfvy=grid%msfvy                                    &
2133                ,msftx=grid%msftx                                    &
2134                ,msfty=grid%msfty                                    &
2135                ,f=grid%f                                            &
2136                ,e=grid%e                                            &
2137                ,sina=grid%sina                                      &
2138                ,cosa=grid%cosa                                      &
2139                !  Terrestrial data
2140                ,rho=grid%rho                                        &
2141                ,dz8w=dz8w                                           &
2142                ,ht=grid%ht                                          &
2143                !  Namelist info
2144                ,use_theta_m=config_flags%use_theta_m                &
2145                ,qc=moist(:,:,:,P_QC)                                &
2146                ,rainc=grid%rainc                                    &
2147                ,rainnc=grid%rainnc                                  &
2148                ,snownc=grid%snownc                                  &
2149                ,graupelnc=grid%graupelnc                            &
2150                ,hailnc=grid%hailnc                                  &
2151                !  The diagnostic output variables
2152                ,sealevelp=grid%sealevelp                            &
2153                ,temperature=grid%temperature                        &
2154                ,pressure=grid%pressure                              &
2155                ,geoheight=grid%geoheight                            &
2156                ,umet=grid%umet                                      &
2157                ,vmet=grid%vmet                                      &
2158                ,speed=grid%speed                                    &
2159                ,dir=grid%dir                                        &
2160                ,psfc=grid%psfc                                      &
2161                ,rain=grid%rain                                      &
2162                ,liqrain=grid%liqrain                                &
2163                ,tpw=grid%tpw                                        &
2164                ,potential_t=grid%potential_t                        &
2165                ,rh=grid%rh                                          &
2166                !  Various indexes for declarations, loop bounds
2167                ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde   &
2168                ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme   &
2169                ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe   &
2170                ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte   )
2171    ENDIF
2172
2173
2174 !  FIRE
2175 if(config_flags%ifire.eq.2)then
2176
2177    call fire_driver_em_init ( grid , config_flags    &
2178             ,ids,ide, kds,kde, jds,jde                &
2179             ,ims,ime, kms,kme, jms,jme                &
2180             ,ips,ipe, kps,kpe, jps,jpe )
2181
2182    CALL wrf_debug ( 100 , 'start_domain_em: After call to fire_driver_em_init' )
2183
2184 !-----------------------------------------------------------------------
2185 ! fire spotting (passive Lagrangian particle transport, tracks
2186 ! firebrand physics properties)
2187 !-----------------------------------------------------------------------
2188
2189    IF (config_flags%fs_firebrand_gen_lim > 0) THEN
2190
2191        IF (config_flags%max_dom == grid%id) THEN
2192
2193            WRITE (message,*) 'SPFire_init: In inner most grid_id: ', grid%id, ' of ', config_flags%max_dom
2194            CALL wrf_debug ( 100 ,  message)
2195            CALL wrf_debug ( 100 , 'start_em: calling firebrand_spotting_em_init ...' )
2196
2197            CALL firebrand_spotting_em_init ( &
2198                grid          = grid,                &
2199                cf            = config_flags,        &
2200                fs_p_id       = grid%fs_p_id,        &
2201                fs_p_src      = grid%fs_p_src,       &
2202                fs_p_dt       = grid%fs_p_dt,        &
2203                fs_p_x        = grid%fs_p_x,         &
2204                fs_p_y        = grid%fs_p_y,         &
2205                fs_p_z        = grid%fs_p_z,         &
2206                fs_p_mass     = grid%fs_p_mass,      &
2207                fs_p_diam     = grid%fs_p_diam,      &
2208                fs_p_effd     = grid%fs_p_effd,      &
2209                fs_p_temp     = grid%fs_p_temp,      &
2210                fs_p_tvel     = grid%fs_p_tvel,      &
2211                fs_last_gen_dt= grid%fs_last_gen_dt, &
2212                fs_gen_idmax  = grid%fs_gen_idmax,   &
2213                fs_count_landed_all   = grid%fs_count_landed_all, &
2214                fs_count_landed_hist  = grid%fs_count_landed_hist,&
2215                fs_landing_mask       = grid%fs_landing_mask,     &
2216                fs_spotting_lkhd      = grid%fs_spotting_lkhd,    &
2217                fs_frac_landed        = grid%fs_frac_landed,      &
2218                fs_count_reset        = grid%fs_count_reset)
2219
2220            CALL wrf_debug ( 100 , 'start_em: firebrand_spotting_em_init ok ' )
2221        ENDIF
2222
2223    ENDIF
2224 !-----------------------------------------------------------------------
2225
2226 endif
2227
2228 #if ( WRFPLUS != 1 )
2229 if( grid%traj_opt /= no_trajectory ) then
2230      call trajectory_init( grid, config_flags, &
2231                            ims,ime, jms,jme, kms,kme )
2232      CALL wrf_debug ( 100 , 'start_domain_em: After call to trajectory_init' )
2233 endif
2234 #endif
2235
2236
2237      CALL wrf_debug ( 100 , 'start_domain_em: Returning' )
2238
2239      RETURN
2240
2241    END SUBROUTINE start_domain_em
2242
2243 !-------------------------------------------------------------------
2244
2245    SUBROUTINE rebalance_driver_cycl ( grid )
2246
2247       USE module_domain, ONLY : domain
2248       IMPLICIT NONE
2249
2250       TYPE (domain)          :: grid
2251
2252       CALL rebalance_cycl( grid &
2253 !
2254 #include "actual_new_args.inc"
2255 !
2256       )
2257
2258    END SUBROUTINE rebalance_driver_cycl
2259
2260 !---------------------------------------------------------------------
2261
2262    SUBROUTINE rebalance_cycl ( grid  &
2263 !
2264 #include "dummy_new_args.inc"
2265 !
2266                         )
2267       USE module_domain, ONLY : domain
2268       USE module_configure, ONLY : grid_config_rec_type, model_config_rec
2269       USE module_model_constants
2270       USE module_state_description
2271       USE module_driver_constants, ONLY: DATA_ORDER_XYZ, DATA_ORDER_YXZ, DATA_ORDER_ZXY, &
2272                                          DATA_ORDER_ZYX, DATA_ORDER_XZY, DATA_ORDER_YZX, &
2273                                          DATA_ORDER_XY, DATA_ORDER_YX, model_data_order
2274
2275 #ifdef DM_PARALLEL
2276       USE module_comm_dm, ONLY : &
2277                            HALO_EM_INIT_1_sub   &
2278                           ,HALO_EM_INIT_2_sub   &
2279                           ,HALO_EM_INIT_3_sub   &
2280                           ,HALO_EM_INIT_4_sub   &
2281                           ,HALO_EM_INIT_5_sub   &
2282                           ,HALO_EM_VINTERP_UV_1_sub
2283 #endif
2284 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
2285    USE module_dm, ONLY : ntasks_x, ntasks_y, ntasks, mytask, local_communicator
2286 #endif
2287
2288       IMPLICIT NONE
2289
2290       TYPE (domain)          :: grid
2291
2292 #include "dummy_new_decl.inc"
2293
2294       TYPE (grid_config_rec_type)              :: config_flags
2295
2296       REAL :: p_surf ,  pd_surf, p_surf_int , pb_int , ht_hold
2297       REAL :: qvf , qvf1 , qvf2, qtot
2298       REAL :: pfu, pfd, phm
2299       REAL :: z0, z1, z2, w1, w2
2300
2301       !  Local domain indices and counters.
2302
2303       INTEGER :: n_moist
2304
2305       INTEGER                             ::                       &
2306                                      ids, ide, jds, jde, kds, kde, &
2307                                      ims, ime, jms, jme, kms, kme, &
2308                                      its, ite, jts, jte, kts, kte, &
2309                                      ips, ipe, jps, jpe, kps, kpe, &
2310                                      i, j, k, kk, ispe, ktf
2311
2312       SELECT CASE ( model_data_order )
2313          CASE ( DATA_ORDER_ZXY )
2314             kds = grid%sd31 ; kde = grid%ed31 ;
2315             ids = grid%sd32 ; ide = grid%ed32 ;
2316             jds = grid%sd33 ; jde = grid%ed33 ;
2317
2318             kms = grid%sm31 ; kme = grid%em31 ;
2319             ims = grid%sm32 ; ime = grid%em32 ;
2320             jms = grid%sm33 ; jme = grid%em33 ;
2321
2322             kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
2323             its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
2324             jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
2325
2326          CASE ( DATA_ORDER_XYZ )
2327             ids = grid%sd31 ; ide = grid%ed31 ;
2328             jds = grid%sd32 ; jde = grid%ed32 ;
2329             kds = grid%sd33 ; kde = grid%ed33 ;
2330
2331             ims = grid%sm31 ; ime = grid%em31 ;
2332             jms = grid%sm32 ; jme = grid%em32 ;
2333             kms = grid%sm33 ; kme = grid%em33 ;
2334
2335             its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
2336             jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
2337             kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
2338
2339          CASE ( DATA_ORDER_XZY )
2340             ids = grid%sd31 ; ide = grid%ed31 ;
2341             kds = grid%sd32 ; kde = grid%ed32 ;
2342             jds = grid%sd33 ; jde = grid%ed33 ;
2343
2344             ims = grid%sm31 ; ime = grid%em31 ;
2345             kms = grid%sm32 ; kme = grid%em32 ;
2346             jms = grid%sm33 ; jme = grid%em33 ;
2347
2348             its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
2349             kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
2350             jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
2351
2352       END SELECT
2353
2354             ktf=MIN(kte,kde-1)
2355
2356       DO j=jts,jte
2357          DO i=its,ite
2358             grid%ph_2(i,1,j) = 0.
2359          END DO
2360       END DO
2361
2362       !  Base state potential temperature and inverse density (alpha = 1/rho) from
2363       !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
2364       !  from equation of state.  The potential temperature is a perturbation from t0.
2365
2366          n_moist = num_moist
2367
2368       print *,'n_moist,PARAM_FIRST_SCALAR',n_moist,PARAM_FIRST_SCALAR
2369
2370        DO j = jts, MIN(jte,jde-1)
2371          DO i = its, MIN(ite,ide-1)
2372
2373         IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
2374
2375                !  Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2376                !  equation) down from the top to get the pressure perturbation.  First get the pressure
2377                !  perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2378
2379                kk = kte - 1
2380                k=kk+1
2381
2382                qtot = 0.
2383                DO ispe=PARAM_FIRST_SCALAR,n_moist
2384                  qtot = qtot + 0.5*(moist(i,kk,j,ispe)+moist(i,kk,j,ispe))
2385                ENDDO
2386                qvf2 = 1./(1.+qtot)
2387                qvf1 = qtot*qvf2
2388
2389                grid%p(i,kk,j) = - 0.5*((grid%c1f(k)*grid%Mu_2(i,j))+qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/grid%rdnw(kk)/qvf2
2390                IF ( grid%use_theta_m == 0) THEN
2391                   qvf = 1.+rvovrd*moist(i,kk,j,P_QV)
2392                ELSE
2393                   qvf = 1.
2394                ENDIF
2395                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf*         &
2396                       (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2397                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2398                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2399
2400
2401                !  Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2402                !  inverse density fields (total and perturbation).
2403
2404         DO kk=kte-2,kts,-1
2405              k = kk + 1
2406
2407              qtot = 0.
2408              DO ispe=PARAM_FIRST_SCALAR,n_moist
2409                qtot = qtot + 0.5*(  moist(i,kk  ,j,ispe) + moist(i,kk+1,j,ispe) )
2410              ENDDO
2411                qvf2 = 1./(1.+qtot)
2412                qvf1 = qtot*qvf2
2413                grid%p(i,kk,j) = grid%p(i,kk+1,j) - ((grid%c1f(k)*grid%Mu_2(i,j)) +       &
2414                                qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/qvf2/grid%rdn(kk+1)
2415                IF ( grid%use_theta_m == 0) THEN
2416                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2417                ELSE
2418                   qvf = 1.
2419                ENDIF
2420                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
2421                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2422                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2423                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2424         ENDDO
2425
2426                !  This is the hydrostatic equation used in the model after the
2427                !  small timesteps.  In
2428                !  the model, grid%al (inverse density) is computed from the
2429                !  geopotential.
2430
2431             IF (grid%hypsometric_opt == 1) THEN
2432                   DO kk  = 2,kte
2433                      k = kk - 1
2434                      grid%ph_2(i,kk,j) = grid%ph_2(i,kk-1,j) - &
2435                                    grid%dnw(kk-1) * ( ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_2(i,j)))*grid%al(i,kk-1,j) &
2436                                  + (grid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,kk-1,j) )
2437                      grid%ph0(i,kk,j) = grid%ph_2(i,kk,j) + grid%phb(i,kk,j)
2438                   END DO
2439             ELSE IF (grid%hypsometric_opt == 2) THEN
2440                 ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is
2441                 ! dry pressure.
2442                 ! Note that al*p approximates Rd*T and dLOG(p) does z.
2443                 ! Here T varies mostly linear with z, the first-order
2444                 ! integration produces better result.
2445
2446                   grid%ph_2(i,1,j) = grid%phb(i,1,j)
2447                   DO k = 2,kte
2448                      pfu = grid%c3f(k  )*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k  ) + grid%p_top
2449                      pfd = grid%c3f(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k-1) + grid%p_top
2450                      phm = grid%c3h(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4h(k-1) + grid%p_top
2451                      grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
2452                   END DO
2453                   DO k = 1,kte
2454                      grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
2455                   END DO
2456                DO k = 1,kte
2457                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2458                END DO
2459             END IF ! hypsometric option
2460        ELSE  ! n_moist
2461          kk = kte - 1
2462                k = kk + 1
2463                qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk,j,P_QV))
2464                qvf2 = 1./(1.+qvf1)
2465                qvf1 = qvf1*qvf2
2466                grid%p(i,kk,j) = - 0.5*((grid%c1f(k)*grid%Mu_2(i,j))+qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/grid%rdnw(kk)/qvf2
2467                IF ( grid%use_theta_m == 0) THEN
2468                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2469                ELSE
2470                   qvf = 1.
2471                ENDIF
2472                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf    &
2473                            *(((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2474                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2475                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2476             !  Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2477             !  inverse density fields (total and perturbation).
2478            DO kk=kte-2,kts,-1
2479                k = kk + 1
2480                qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk+1,j,P_QV))
2481                qvf2 = 1./(1.+qvf1)
2482                qvf1 = qvf1*qvf2
2483                grid%p(i,kk,j) = grid%p(i,kk+1,j) - ((grid%c1f(k)*grid%Mu_2(i,j)) + qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/qvf2/grid%rdn(kk+1)
2484                IF ( grid%use_theta_m == 0) THEN
2485                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2486                ELSE
2487                   qvf = 1.
2488                ENDIF
2489                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
2490                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2491                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2492                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2493            ENDDO
2494                !  This is the hydrostatic equation used in the model after the small timesteps.  In
2495                !  the model, grid%al (inverse density) is computed from the geopotential.
2496             IF (grid%hypsometric_opt == 1) THEN
2497                DO k  = 2,kte
2498                   grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
2499                                 grid%dnw(k-1) * ( ((grid%c1h(k)*grid%mub(i,j)+grid%c2h(k))+(grid%c1h(k)*grid%mu_2(i,j)))*grid%al(i,k-1,j) &
2500                               + (grid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,k-1,j) )
2501                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2502                END DO
2503             ELSE IF (grid%hypsometric_opt == 2) THEN
2504              ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry
2505              ! pressure.
2506              ! Note that al*p approximates Rd*T and dLOG(p) does z.
2507              ! Here T varies mostly linear with z, the first-order integration
2508              ! produces better result.
2509                grid%ph_2(i,1,j) = grid%phb(i,1,j)
2510                DO k = 2,kte
2511                   pfu = grid%c3f(k  )*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k  ) + grid%p_top
2512                   pfd = grid%c3f(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k-1) + grid%p_top
2513                   phm = grid%c3h(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4h(k-1) + grid%p_top
2514                   grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
2515                END DO
2516
2517                DO k = 1,kte
2518                   grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
2519                END DO
2520
2521                DO k = 1,kte
2522                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2523                END DO
2524
2525             END IF ! hypsometric
2526        ENDIF ! nmoist
2527
2528 ! update surface pressure PSFC (needed for post-processing):
2529                z0 = grid%ph0(i,1,j)/g
2530                z1 = 0.5*(grid%ph0(i,1,j)+grid%ph0(i,2,j))/g
2531                z2 = 0.5*(grid%ph0(i,2,j)+grid%ph0(i,3,j))/g
2532                w1 = (z0 - z2)/(z1 - z2)
2533                w2 = 1. - w1
2534                grid%psfc(i,j) = w1*(grid%p(i,1,j)+grid%pb(i,1,j))+w2*(grid%p(i,2,j)+grid%pb(i,2,j))
2535
2536          END DO !i
2537         ENDDO !j
2538
2539       ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2540 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
2541 # include "HALO_EM_INIT_1.inc"
2542 # include "HALO_EM_INIT_2.inc"
2543 # include "HALO_EM_INIT_3.inc"
2544 # include "HALO_EM_INIT_4.inc"
2545 # include "HALO_EM_INIT_5.inc"
2546 #endif
2547    END SUBROUTINE rebalance_cycl
2548
2549 !---------------------------------------------------------------------