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                       ,grid%RQCNCUTEN, grid%RQINCUTEN,grid%rliq             &      !mchen add for cammpmg
1239                       ,grid%cldfra_dp,grid%cldfra_sh                        & ! ckay for subgrid cloud
1240                       ,grid%te_temf,grid%cf3d_temf,grid%wm_temf        & ! WA
1241                       ,grid%massflux_EDKF, grid%entr_EDKF, grid%detr_EDKF                &
1242                       ,grid%thl_up,grid%thv_up,grid%rt_up                                &
1243                       ,grid%rv_up,grid%rc_up,grid%u_up,grid%v_up,grid%frac_up            &
1244                       ,grid%RDCASHTEN, grid%RQCDCSHTEN                                   &
1245                       ,grid%cldareaa, grid%cldareab, grid%cldliqa, grid%cldliqb          &
1246                       ,grid%ca_rad, grid%cw_rad                                          &
1247                       ,grid%pblmax, grid%wub, grid%ltopb, grid%clddpthb, grid%cldtopb    &
1248                       ,grid%capesave, grid%ainckfsa, grid%radsave                        &
1249                       ,grid%rainsh, grid%rainshvb, grid%kdcldtop, grid%kdcldbas          &
1250                       ,grid%xtime1, grid%PBLHAVG, grid%TKEAVG                            &
1251                       ,grid%ccn_conc                                          & ! RAS
1252                       ,grid%QKE                                               &!JOE-for mynn
1253                       ,grid%pek_pbl,grid%pep_pbl                              & ! EEPS
1254                       ,grid%ghi_accum                                         & ! Solar diagnostics
1255                       ,grid%landusef,grid%landusef2,grid%mosaic_cat_index                            & ! danli mosaic
1256                       ,grid%TSK_mosaic,grid%TSLB_mosaic,grid%SMOIS_mosaic,grid%SH2O_mosaic           & ! danli mosaic
1257                       ,grid%CANWAT_mosaic,grid%SNOW_mosaic,grid%SNOWH_mosaic,grid%SNOWC_mosaic       & ! danli mosaic
1258                       ,grid%ALBEDO_mosaic,grid%ALBBCK_mosaic, grid%EMISS_mosaic                      & ! danli mosaic
1259                       ,grid%EMBCK_mosaic, grid%ZNT_mosaic, grid%Z0_mosaic                            & ! danli mosaic
1260                       ,grid%TR_URB2D_mosaic,grid%TB_URB2D_mosaic                                     & ! danli mosaic
1261                       ,grid%TG_URB2D_mosaic,grid%TC_URB2D_mosaic                                     & ! danli mosaic
1262                       ,grid%QC_URB2D_mosaic                                                          & ! danli mosaic
1263                       ,grid%TRL_URB3D_mosaic,grid%TBL_URB3D_mosaic                                   & ! danli mosaic
1264                       ,grid%TGL_URB3D_mosaic                                                         & ! danli mosaic
1265                       ,grid%SH_URB2D_mosaic,grid%LH_URB2D_mosaic                                     & ! danli mosaic
1266                       ,grid%G_URB2D_mosaic,grid%RN_URB2D_mosaic                                      & ! danli mosaic
1267                       ,grid%TS_URB2D_mosaic                                                          & ! danli mosaic
1268                       ,grid%TS_RUL2D_mosaic                                                          & ! danli mosaic
1269                       ,grid%irr_rand_field,config_flags%irr_ph,config_flags%irr_freq                 &
1270                       ,grid%QR_CU, grid%QS_CU                                                        &
1271                       ,grid%NC_CU, grid%NI_CU, grid%NR_CU, grid%NS_CU,grid%CCN_CU                    &
1272                       ,grid%alevsiz_cu,grid%num_months,grid%no_src_types_cu,grid%aeromcu             &
1273                       ,grid%aeropcu,grid%EFCG,grid%EFCS,grid%EFIG,grid%EFIS,grid%EFSG,grid%EFSS      &
1274                       )
1275        ENDDO ! loop of tiles for phy_init
1276
1277        if(.not.grid%restart) then
1278
1279          if(config_flags%sf_surface_physics .eq. NOAHMPSCHEME .and. config_flags%opt_run .eq. 5) then
1280               call groundwater_init ( grid ,                                                            &
1281                  config_flags%num_soil_layers, grid%dzs, grid%isltyp, grid%ivgtyp, config_flags%wtddt , &
1282                  grid%fdepthxy   , grid%ht     , grid%riverbedxy, grid%eqzwt     ,                      &
1283                  grid%rivercondxy, grid%pexpxy , grid%areaxy    , grid%zwtxy     ,                      &
1284                  grid%smois      , grid%sh2o   , grid%smoiseq   , grid%smcwtdxy  ,                      &
1285                  grid%qlatxy     , grid%qslatxy, grid%qrfxy     , grid%qrfsxy    ,                      &
1286                  grid%deeprechxy , grid%rechxy , grid%qspringxy , grid%qspringsxy,                      &
1287                  grid%rechclim   ,                                                                      &
1288                  ids,ide, jds,jde, kds,kde,                                                             &
1289                  ims,ime, jms,jme, kms,kme,                                                             &
1290                  ips,ipe, jps,jpe, kps,kpe,                                                             &
1291                  its,ite, jts,jte, kts,kte )
1292          endif
1293
1294        endif
1295
1296     ENDIF   ! no phy_init for the backwards part of the DFI
1297
1298    CALL wrf_debug ( 100 , 'start_domain_em: After call to phy_init' )
1299
1300    IF (config_flags%do_avgflx_em .EQ. 1) THEN
1301       WRITE ( message , FMT = '("start_em: initializing avgflx on domain ",I3)' ) &
1302            & grid%id
1303       CALL wrf_message(trim(message))
1304       grid%avgflx_count = 0
1305       f_flux = config_flags%do_avgflx_cugd .EQ. 1  ! WA 9/24/10
1306       DO ij = 1, grid%num_tiles
1307          call wrf_debug(200,'In start_em, before zero_avgflx call')
1308          if (.not. grid%restart) call zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
1309               &   ids, ide, jds, jde, kds, kde,           &
1310               &   ims, ime, jms, jme, kms, kme,           &
1311               &   grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, f_flux, &
1312               &   grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
1313               &   grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
1314          call wrf_debug(200,'In start_em, after zero_avgflx call')
1315       ENDDO
1316    ENDIF
1317
1318    call wrf_debug(100,'start_em: calling lightning_init')
1319    CALL lightning_init ( id=grid%id, itimestep=grid%itimestep, restart=grid%restart, dt=grid%dt, dx=grid%dx   &
1320                     ! Namelist control options
1321                         ,cu_physics=config_flags%cu_physics,mp_physics=config_flags%mp_physics    &
1322                         ,do_radar_ref=config_flags%do_radar_ref                                   &
1323                         ,lightning_option=config_flags%lightning_option                           &
1324                         ,lightning_dt=config_flags%lightning_dt                                   &
1325                         ,lightning_start_seconds=config_flags%lightning_start_seconds             &
1326                         ,ltngacttime=grid%ltngacttime                                             &
1327                         ,iccg_prescribed_num=config_flags%iccg_prescribed_num                     &
1328                         ,iccg_prescribed_den=config_flags%iccg_prescribed_den                     &
1329                         ,cellcount_method=config_flags%cellcount_method                           &
1330                     ! Order dependent args for domain, mem, and tile dims
1331                         ,ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde             &
1332                         ,ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme             &
1333                         ,its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte             &
1334                     ! IC and CG flash rates and accumulated flash count
1335                         ,ic_flashcount=grid%ic_flashcount, ic_flashrate=grid%ic_flashrate        &
1336                         ,cg_flashcount=grid%cg_flashcount, cg_flashrate=grid%cg_flashrate        &
1337 #if (WRF_CHEM == 1)
1338                         ,lnox_opt=config_flags%lnox_opt           &
1339                         ,lnox_passive=config_flags%lnox_passive   &
1340                         ,lnox_total=tracer(:,:,:,p_lnox_total)    &
1341                         ,lnox_ic=tracer(:,:,:,p_lnox_ic)          &
1342                         ,lnox_cg=tracer(:,:,:,p_lnox_cg)          &
1343 #endif
1344                     )
1345
1346     call wrf_debug(100,'start_em: after calling lightning_init')
1347
1348    END IF ! restart
1349
1350 #if 0
1351 #include "CYCLE_TEST.inc"
1352 #endif
1353
1354 !
1355 !
1356
1357  !  set physical boundary conditions for all initialized variables
1358
1359 !-----------------------------------------------------------------------
1360 !  Stencils for patch communications  (WCS, 29 June 2001)
1361 !  Note:  the size of this halo exchange reflects the
1362 !         fact that we are carrying the uncoupled variables
1363 !         as state variables in the mass coordinate model, as
1364 !         opposed to the coupled variables as in the height
1365 !         coordinate model.
1366 !
1367 !                           * * * * *
1368 !         *        * * *    * * * * *
1369 !       * + *      * + *    * * + * *
1370 !         *        * * *    * * * * *
1371 !                           * * * * *
1372 !
1373 !j  grid%u_1                          x
1374 !j  grid%u_2                          x
1375 !j  grid%v_1                          x
1376 !j  grid%v_2                          x
1377 !j  grid%w_1                          x
1378 !j  grid%w_2                          x
1379 !j  grid%t_1                          x
1380 !j  grid%t_2                          x
1381 !j grid%ph_1                          x
1382 !j grid%ph_2                          x
1383 !
1384 !j  grid%t_init                       x
1385 !
1386 !j  grid%phb   x
1387 !j  grid%ph0   x
1388 !j  grid%php   x
1389 !j   grid%pb   x
1390 !j   grid%al   x
1391 !j  grid%alt   x
1392 !j  grid%alb   x
1393 !
1394 !  the following are 2D (xy) variables
1395 !
1396 !j grid%mu_1                          x
1397 !j grid%mu_2                          x
1398 !j grid%mub   x
1399 !j grid%mu0   x
1400 !j grid%ht    x
1401 !j grid%msftx x
1402 !j grid%msfty x
1403 !j grid%msfux x
1404 !j grid%msfuy x
1405 !j grid%msfvx x
1406 !j grid%msfvy x
1407 !j grid%sina  x
1408 !j grid%cosa  x
1409 !j grid%e     x
1410 !j grid%f     x
1411 !
1412 !  4D variables
1413 !
1414 ! moist                        x
1415 !  chem                        x
1416 !scalar                        x
1417
1418 !--------------------------------------------------------------
1419
1420 #ifdef DM_PARALLEL
1421 # include "HALO_EM_INIT_1.inc"
1422 # include "HALO_EM_INIT_2.inc"
1423 # include "HALO_EM_INIT_3.inc"
1424 # include "HALO_EM_INIT_4.inc"
1425 # include "HALO_EM_INIT_5.inc"
1426       IF ( config_flags%sf_ocean_physics .EQ. PWP3DSCHEME ) THEN
1427 # include "HALO_EM_INIT_6.inc"
1428       END IF
1429 # include "PERIOD_BDY_EM_INIT.inc"
1430 # include "PERIOD_BDY_EM_MOIST.inc"
1431 # include "PERIOD_BDY_EM_TKE.inc"
1432 # include "PERIOD_BDY_EM_SCALAR.inc"
1433 # include "PERIOD_BDY_EM_CHEM.inc"
1434 if(config_flags%sf_surface_physics.eq.NOAHMPSCHEME.and.config_flags%opt_run.eq.5)then
1435 # include "HALO_EM_HYDRO_NOAHMP_INIT.inc"
1436 endif
1437 #endif
1438
1439
1440    CALL set_physical_bc3d( grid%u_1 , 'U' , config_flags ,                  &
1441                          ids , ide , jds , jde , kds , kde ,        &
1442                          ims , ime , jms , jme , kms , kme ,        &
1443                          its , ite , jts , jte , kts , kte ,        &
1444                          its , ite , jts , jte , kts , kte )
1445    CALL set_physical_bc3d( grid%u_2 , 'U' , config_flags ,                  &
1446                          ids , ide , jds , jde , kds , kde ,        &
1447                          ims , ime , jms , jme , kms , kme ,        &
1448                          its , ite , jts , jte , kts , kte ,        &
1449                          its , ite , jts , jte , kts , kte )
1450
1451    CALL set_physical_bc3d( grid%v_1 , 'V' , config_flags ,                  &
1452                          ids , ide , jds , jde , kds , kde ,        &
1453                          ims , ime , jms , jme , kms , kme ,        &
1454                          its , ite , jts , jte , kts , kte ,        &
1455                          its , ite , jts , jte , kts , kte )
1456    CALL set_physical_bc3d( grid%v_2 , 'V' , config_flags ,                  &
1457                          ids , ide , jds , jde , kds , kde ,        &
1458                          ims , ime , jms , jme , kms , kme ,        &
1459                          its , ite , jts , jte , kts , kte ,        &
1460                          its , ite , jts , jte , kts , kte )
1461
1462 ! set kinematic condition for w
1463
1464    CALL set_physical_bc2d( grid%ht , 'r' , config_flags ,                &
1465                            ids , ide , jds , jde , &
1466                            ims , ime , jms , jme , &
1467                            its , ite , jts , jte , &
1468                            its , ite , jts , jte   )
1469
1470 ! Set the w at the surface.  If this is the start of a forecast, or if this is a cycled run
1471 ! and the start of that forecast, we define the w field.  However, if this a restart, then
1472 ! we leave w alone.  Initial value is V dot grad(topo) at surface, then smoothly decaying
1473 ! above that.
1474
1475    IF ( ( start_of_simulation .OR. config_flags%cycling ) .AND. ( .NOT. config_flags%restart ) ) THEN
1476
1477 ! If W already exists (not zero), then we leave it alone.  How to do this?  We find the
1478 ! max/min on this node at the surface.  If parallel, we collect the max/min from all procs.
1479 ! If the max/min throughout the entire domain at the surface is identically 0, then we say
1480 ! that the W field is NOT initialized, and we run the set_w_surface routines for the
1481 ! two time levels of W.  If the field is already initialized, we do NOT run those two
1482 ! routines.
1483
1484       w_max = grid%w_2(its,1,jts)
1485       w_min = grid%w_2(its,1,jts)
1486       DO j = jts, MIN(jte,jde-1)
1487          DO i = its, MIN(ite,ide-1)
1488          w_max = MAX ( w_max , grid%w_2(i,1,j) )
1489          w_min = MIN ( w_min , grid%w_2(i,1,j) )
1490          END DO
1491       END DO
1492 #ifdef DM_PARALLEL
1493       w_max = wrf_dm_max_real ( w_max )
1494       w_min = wrf_dm_min_real ( w_min )
1495 #endif
1496
1497       IF ( ( ABS(w_max) .LT. 1.E-6 ) .AND. &
1498            ( ABS(w_min) .LT. 1.E-6 ) ) THEN
1499          w_needs_to_be_set = .TRUE.
1500          ! Force use_input_w to false, if there is no input w
1501 #if ( WRFPLUS == 1 )
1502          IF ( config_flags%use_input_w ) THEN
1503             config_flags%use_input_w = .FALSE.
1504             CALL nl_set_use_input_w (grid%id, .FALSE.)
1505          ENDIF
1506 #endif
1507       ELSE
1508          IF ( config_flags%use_input_w ) THEN
1509             w_needs_to_be_set = .FALSE.
1510          ELSE
1511             w_needs_to_be_set = .TRUE.
1512          END IF
1513       END IF
1514
1515       IF ( w_needs_to_be_set ) THEN
1516
1517          ! setting kinematic condition for w at the surface
1518
1519          fill_w_flag = .true.
1520          CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
1521                               grid%w_1, grid%ht, grid%u_1, grid%v_1, grid%cf1, &
1522                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
1523                               ids, ide, jds, jde, kds, kde,                    &
1524                               ims, ime, jms, jme, kms, kme,                    &
1525                               its, ite, jts, jte, kts, kte                     )
1526          CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
1527                               grid%w_2, grid%ht, grid%u_2, grid%v_2, grid%cf1, &
1528                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
1529                               ids, ide, jds, jde, kds, kde,                    &
1530                               ims, ime, jms, jme, kms, kme,                    &
1531                               its, ite, jts, jte, kts, kte                     )
1532       END IF
1533    END IF
1534
1535    IF ( .NOT. config_flags%restart ) THEN
1536 ! set up slope-radiation constant arrays based on topography
1537 ! paj: also computes laplacian for topo_wind
1538 !
1539        DO j = jts,min(jte,jde-1)
1540        DO i = its, min(ite,ide-1)
1541           if (config_flags%periodic_x ) then
1542             im1 = i-1
1543             ip1 = i+1
1544           else
1545             im1 = max(i-1,ids)
1546             ip1 = min(i+1,ide-1)
1547           endif
1548
1549           if (config_flags%periodic_y ) then
1550             jm1 = j-1
1551             jp1 = j+1
1552           else
1553             jm1 = max(j-1,jds)
1554             jp1 = min(j+1,jde-1)
1555           endif
1556           grid%toposlpx(i,j)=(grid%ht(ip1,j)-grid%ht(im1,j))*grid%msftx(i,j)*grid%rdx/(ip1-im1)
1557           grid%toposlpy(i,j)=(grid%ht(i,jp1)-grid%ht(i,jm1))*grid%msfty(i,j)*grid%rdy/(jp1-jm1)
1558           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.
1559              hx = grid%toposlpx(i,j)
1560              hy = grid%toposlpy(i,j)
1561              pi = 4.*atan(1.)
1562              grid%slope(i,j) = atan((hx**2+hy**2)**.5)
1563              if (grid%slope(i,j).lt.1.e-4) then
1564                grid%slope(i,j) = 0.
1565                grid%slp_azi(i,j) = 0.
1566              else
1567                grid%slp_azi(i,j) = atan2(hx,hy)+pi
1568
1569 ! Rotate slope azimuth to lat-lon grid
1570                if (grid%cosa(i,j).ge.0) then
1571                  grid%slp_azi(i,j) = grid%slp_azi(i,j) - asin(grid%sina(i,j))
1572                else
1573                  grid%slp_azi(i,j) = grid%slp_azi(i,j) - (pi - asin(grid%sina(i,j)))
1574                endif
1575              endif
1576        ENDDO
1577        ENDDO
1578 !
1579 ! paj: computes ct and ct2 for topo_wind
1580 !
1581        grid%ctopo=1.
1582        grid%ctopo2=1.
1583 !
1584        if (config_flags%topo_wind.eq.1) then
1585          DO j = jts,min(jte,jde-1)
1586          DO i = its, min(ite,ide-1)
1587             if(grid%xland(i,j).lt.1.5)grid%ctopo(i,j)=sqrt(grid%var_sso(i,j))
1588             if (grid%ctopo(i,j).le.2.718) then
1589             grid%ctopo(i,j)=1.
1590             else
1591             grid%ctopo(i,j)=alog(grid%ctopo(i,j))
1592             endif
1593 !
1594             if (grid%lap_hgt(i,j).gt.-10.) then
1595               grid%ctopo(i,j)=grid%ctopo(i,j)
1596             else
1597                if (grid%lap_hgt(i,j).ge.-20) then
1598                  alpha=(grid%lap_hgt(i,j)+20.)/10.
1599                  grid%ctopo(i,j)=alpha*grid%ctopo(i,j)+(1-alpha)
1600                else
1601                   if (grid%lap_hgt(i,j).ge.-30.) then
1602                   grid%ctopo(i,j)=(grid%lap_hgt(i,j)+30.)/10.
1603                   grid%ctopo2(i,j)=(grid%lap_hgt(i,j)+30.)/10.
1604                   else
1605                   grid%ctopo(i,j)=0.
1606                   grid%ctopo2(i,j)=0.
1607                   endif
1608                endif
1609             endif
1610          ENDDO
1611          ENDDO
1612 !       if topo_wind==2 (D.Ovens, C.Mass)
1613        else if (config_flags%topo_wind.eq.2) then
1614          DO j = jts,min(jte,jde-1)
1615          DO i = its, min(ite,ide-1)
1616            if (grid%xland(i,j).lt.1.5) then
1617               vfac = amin1(1.575,(grid%var2d(i,j)*0.4/200.+1.175))
1618               vfac = vfac * vfac
1619            else                      ! over water, leave it alone
1620               vfac = 1.
1621            endif
1622 !          print *, j,i,grid%ctopo(i,j),vfac
1623            grid%ctopo(i,j)=grid%ctopo(i,j)*vfac
1624          ENDDO
1625          ENDDO
1626        endif
1627
1628    END IF
1629
1630    CALL set_physical_bc3d( grid%w_1 , 'W' , config_flags ,                 &
1631                          ids , ide , jds , jde , kds , kde ,        &
1632                          ims , ime , jms , jme , kms , kme ,        &
1633                          its , ite , jts , jte , kts , kte ,        &
1634                          its , ite , jts , jte , kts , kte )
1635    CALL set_physical_bc3d( grid%w_2 , 'W' , config_flags ,                 &
1636                          ids , ide , jds , jde , kds , kde ,        &
1637                          ims , ime , jms , jme , kms , kme ,        &
1638                          its , ite , jts , jte , kts , kte ,        &
1639                          its , ite , jts , jte , kts , kte )
1640
1641    CALL set_physical_bc3d( grid%ph_1 , 'W' , config_flags ,                 &
1642                          ids , ide , jds , jde , kds , kde ,        &
1643                          ims , ime , jms , jme , kms , kme ,        &
1644                          its , ite , jts , jte , kts , kte ,        &
1645                          its , ite , jts , jte , kts , kte )
1646
1647    CALL set_physical_bc3d( grid%ph_2 , 'W' , config_flags ,                 &
1648                          ids , ide , jds , jde , kds , kde ,        &
1649                          ims , ime , jms , jme , kms , kme ,        &
1650                          its , ite , jts , jte , kts , kte ,        &
1651                          its , ite , jts , jte , kts , kte )
1652
1653    CALL set_physical_bc3d( grid%t_1 , 't' , config_flags ,                 &
1654                          ids , ide , jds , jde , kds , kde ,        &
1655                          ims , ime , jms , jme , kms , kme ,        &
1656                          its , ite , jts , jte , kts , kte ,        &
1657                          its , ite , jts , jte , kts , kte )
1658
1659    CALL set_physical_bc3d( grid%t_2 , 't' , config_flags ,                 &
1660                          ids , ide , jds , jde , kds , kde ,        &
1661                          ims , ime , jms , jme , kms , kme ,        &
1662                          its , ite , jts , jte , kts , kte ,        &
1663                          its , ite , jts , jte , kts , kte )
1664
1665    CALL set_physical_bc2d( grid%mu_1, 't' , config_flags ,   &
1666                            ids , ide , jds , jde ,  &
1667                            ims , ime , jms , jme ,  &
1668                            its , ite , jts , jte ,  &
1669                            its , ite , jts , jte   )
1670    CALL set_physical_bc2d( grid%mu_2, 't' , config_flags ,   &
1671                            ids , ide , jds , jde ,  &
1672                            ims , ime , jms , jme ,  &
1673                            its , ite , jts , jte ,  &
1674                            its , ite , jts , jte   )
1675    CALL set_physical_bc2d( grid%mub , 't' , config_flags ,   &
1676                            ids , ide , jds , jde ,  &
1677                            ims , ime , jms , jme ,  &
1678                            its , ite , jts , jte ,  &
1679                            its , ite , jts , jte   )
1680 !  CALL set_physical_bc2d( grid%mu0 , 't' , config_flags ,   &
1681 !                          ids , ide , jds , jde ,  &
1682 !                          ims , ime , jms , jme ,  &
1683 !                          its , ite , jts , jte ,  &
1684 !                          its , ite , jts , jte   )
1685
1686
1687    CALL set_physical_bc3d( grid%phb , 'W' , config_flags ,                 &
1688                          ids , ide , jds , jde , kds , kde ,        &
1689                          ims , ime , jms , jme , kms , kme ,        &
1690                          its , ite , jts , jte , kts , kte ,        &
1691                          its , ite , jts , jte , kts , kte )
1692    CALL set_physical_bc3d( grid%ph0 , 'W' , config_flags ,                 &
1693                          ids , ide , jds , jde , kds , kde ,        &
1694                          ims , ime , jms , jme , kms , kme ,        &
1695                          its , ite , jts , jte , kts , kte ,        &
1696                          its , ite , jts , jte , kts , kte )
1697    CALL set_physical_bc3d( grid%php , 'W' , config_flags ,                 &
1698                          ids , ide , jds , jde , kds , kde ,        &
1699                          ims , ime , jms , jme , kms , kme ,        &
1700                          its , ite , jts , jte , kts , kte ,        &
1701                          its , ite , jts , jte , kts , kte )
1702
1703    CALL set_physical_bc3d( grid%pb , 't' , config_flags ,                 &
1704                          ids , ide , jds , jde , kds , kde ,        &
1705                          ims , ime , jms , jme , kms , kme ,        &
1706                          its , ite , jts , jte , kts , kte ,        &
1707                          its , ite , jts , jte , kts , kte )
1708    CALL set_physical_bc3d( grid%al , 't' , config_flags ,                 &
1709                          ids , ide , jds , jde , kds , kde ,        &
1710                          ims , ime , jms , jme , kms , kme ,        &
1711                          its , ite , jts , jte , kts , kte ,        &
1712                          its , ite , jts , jte , kts , kte )
1713    CALL set_physical_bc3d( grid%alt , 't' , config_flags ,                 &
1714                          ids , ide , jds , jde , kds , kde ,        &
1715                          ims , ime , jms , jme , kms , kme ,        &
1716                          its , ite , jts , jte , kts , kte ,        &
1717                          its , ite , jts , jte , kts , kte )
1718    CALL set_physical_bc3d( grid%alb , 't' , config_flags ,                 &
1719                          ids , ide , jds , jde , kds , kde ,        &
1720                          ims , ime , jms , jme , kms , kme ,        &
1721                          its , ite , jts , jte , kts , kte ,        &
1722                          its , ite , jts , jte , kts , kte )
1723    CALL set_physical_bc3d(grid%t_init, 't' , config_flags ,                 &
1724                          ids , ide , jds , jde , kds , kde ,        &
1725                          ims , ime , jms , jme , kms , kme ,        &
1726                          its , ite , jts , jte , kts , kte ,        &
1727                          its , ite , jts , jte , kts , kte )
1728    CALL set_physical_bc3d(grid%tke_2, 't' , config_flags ,                 &
1729                          ids , ide , jds , jde , kds , kde ,        &
1730                          ims , ime , jms , jme , kms , kme ,        &
1731                          its , ite , jts , jte , kts , kte ,        &
1732                          its , ite , jts , jte , kts , kte )
1733
1734    IF (num_moist > 0) THEN
1735
1736 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1737
1738       loop_3d_m   : DO loop = 1 , num_moist
1739          CALL set_physical_bc3d( moist(:,:,:,loop) , 'r' , config_flags ,                 &
1740                                  ids , ide , jds , jde , kds , kde ,        &
1741                                  ims , ime , jms , jme , kms , kme ,        &
1742                                  its , ite , jts , jte , kts , kte ,        &
1743                                  its , ite , jts , jte , kts , kte )
1744       END DO loop_3d_m
1745
1746    ENDIF
1747
1748
1749 ! Some initializations for the simple ccn field.
1750
1751    IF ( f_qnn ) THEN
1752       IF ( config_flags%mp_physics == wdm5scheme .or. config_flags%mp_physics == wdm6scheme ) THEN
1753          ! NO OP
1754       ELSE IF ( config_flags%mp_physics == nssl_2mom ) THEN
1755          IF ( config_flags%nssl_ccn_is_ccna == 0 ) THEN
1756            grid%ccn_conc =  config_flags%nssl_cccn/1.225
1757          ELSE
1758            grid%ccn_conc = 0
1759          ENDIF
1760       ELSE
1761          ! NO OP
1762       END IF
1763       ccn_max_val = MAXVAL(scalar(its:MIN(ite,ide-1),kts:kte-1,jts:MIN(jte,jde-1),p_qnn))
1764 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1765       ccn_max_val = wrf_dm_max_real ( ccn_max_val )
1766 #endif
1767       IF ( ccn_max_val < 1.0 ) THEN ! initialization of ccn not already done
1768          DO j=jts,MIN(jte,jde-1)
1769             DO k=kts,kte
1770                DO i=its,MIN(ite,ide-1)
1771                   scalar(i,k,j,p_qnn) = grid%ccn_conc
1772                END DO
1773             END DO
1774          END DO
1775       END IF
1776    END IF
1777
1778    IF (num_scalar > 0) THEN
1779
1780 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1781
1782       loop_3d_s   : DO loop = 1 , num_scalar
1783          CALL set_physical_bc3d( scalar(:,:,:,loop) , 'r' , config_flags ,                 &
1784                                  ids , ide , jds , jde , kds , kde ,        &
1785                                  ims , ime , jms , jme , kms , kme ,        &
1786                                  its , ite , jts , jte , kts , kte ,        &
1787                                  its , ite , jts , jte , kts , kte )
1788       END DO loop_3d_s
1789
1790    ENDIF
1791
1792
1793 #if (WRF_CHEM == 1)
1794         if(config_flags%tracer_opt > 0 )then
1795            call initialize_tracer (tracer,config_flags%chem_in_opt, &
1796                                config_flags%tracer_opt,num_tracer,  &
1797                                ids,ide, jds,jde, kds,kde,      & ! domain dims
1798                                ims,ime, jms,jme, kms,kme,      & ! memory dims
1799                                ips,ipe, jps,jpe, kps,kpe,      & ! patch  dims
1800                                its,ite, jts,jte, kts,kte )
1801
1802         endif
1803 !
1804 !   we do this here, so we only have one chem_init routine for either core....
1805 !
1806          do j=jts,min(jte,jde-1)
1807             do i=its,min(ite,ide-1)
1808                do k=kts,kte
1809                  z_at_w(i,k,j)=(grid%ph_2(i,k,j)+grid%phb(i,k,j))/g
1810                enddo
1811                do k=kts,min(kte,kde-1)
1812                  tempfac=(grid%t_1(i,k,j) + t0)*((grid%p(i,k,j) + grid%pb(i,k,j))/p1000mb)**rcp
1813                  convfac(i,k,j) = (grid%p(i,k,j)+grid%pb(i,k,j))/rgasuniv/tempfac
1814                enddo
1815             enddo
1816          enddo
1817
1818          CALL chem_init (grid%id,chem,emis_ant,scalar,grid%dt,grid%bioemdt,grid%photdt,     &
1819                          grid%chemdt,                                       &
1820                          grid%stepbioe,grid%stepphot,grid%stepchem,grid%stepfirepl,      &
1821                          grid%plumerisefire_frq,z_at_w,grid%xlat,grid%xlong,g,           &
1822                          grid%aerwrf,config_flags,grid,            &
1823                          grid%alt,grid%t_1,grid%p,convfac,grid%ttday,grid%tcosz,         &
1824                          grid%julday,grid%gmt,&
1825                          grid%tauaer1,grid%tauaer2,grid%tauaer3,grid%tauaer4,                   &
1826                          grid%gaer1,grid%gaer2,grid%gaer3,grid%gaer4,                           &
1827                          grid%waer1,grid%waer2,grid%waer3,grid%waer4,                           &
1828                          grid%l2aer,grid%l3aer,grid%l4aer,grid%l5aer,grid%l6aer,grid%l7aer,     &
1829                          grid%extaerlw1,grid%extaerlw2,grid%extaerlw3,grid%extaerlw4,           &
1830                          grid%extaerlw5,grid%extaerlw6,grid%extaerlw7,grid%extaerlw8,           &
1831                          grid%extaerlw9,grid%extaerlw10,grid%extaerlw11,grid%extaerlw12,        &
1832                          grid%extaerlw13,grid%extaerlw14,grid%extaerlw15,grid%extaerlw16,       &
1833                          grid%tauaerlw1,grid%tauaerlw2,grid%tauaerlw3,grid%tauaerlw4,           &
1834                          grid%tauaerlw5,grid%tauaerlw6,grid%tauaerlw7,grid%tauaerlw8,           &
1835                          grid%tauaerlw9,grid%tauaerlw10,grid%tauaerlw11,grid%tauaerlw12,        &
1836                          grid%tauaerlw13,grid%tauaerlw14,grid%tauaerlw15,grid%tauaerlw16,       &
1837                          grid%dgnum4d, grid%dgnumwet4d, grid%dgnum_a1, grid%dgnum_a2,           &
1838                          grid%dgnum_a3, grid%dgnumwet_a1, grid%dgnumwet_a2, grid%dgnumwet_a3,   &
1839                          grid%pm2_5_dry,grid%pm2_5_water,grid%pm2_5_dry_ec,                &
1840                          grid%tsoa,grid%asoa,grid%bsoa,                                    &
1841                          grid%last_chem_time_year,grid%last_chem_time_month,               &
1842                          grid%last_chem_time_day,grid%last_chem_time_hour,                 &
1843                          grid%last_chem_time_minute,grid%last_chem_time_second,            &
1844                          grid%chem_in_opt,grid%kemit,grid%num_vert_mix,numgas,             &
1845                          ids , ide , jds , jde , kds , kde ,                &
1846                          ims , ime , jms , jme , kms , kme ,                &
1847                          its , ite , jts , jte , kts , kte                  )
1848 #if( WRF_USE_CLM == 1)
1849 !-------------------------------------------------------------------
1850 !   check that the megan emission species are in chem option
1851 !-------------------------------------------------------------------
1852          if (config_flags%chem_opt > 0) then
1853            write(message,'('' start_domain_em: numgas = '',i4)') numgas
1854            CALL wrf_message( trim(message) )
1855            CALL wrf_spc_chk( numgas )
1856          endif
1857 #endif
1858
1859 !
1860 ! calculate initial pm
1861 !
1862 !       print *,'calculating initial pm'
1863         select case (config_flags%chem_opt)
1864         case (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2)
1865            call sum_pm_gocart (                                             &
1866                 grid%alt, chem, grid%pm2_5_dry, grid%pm2_5_dry_ec,grid%pm10,&
1867                 ids,ide, jds,jde, kds,kde,                                  &
1868                 ims,ime, jms,jme, kms,kme,                                  &
1869                 its,ite, jts,jte, kts,kte-1                                 )
1870         case (RADM2SORG,RADM2SORG_AQ,RADM2SORG_AQCHEM,RADM2SORG_KPP,RACMSORG_AQ,RACMSORG_AQCHEM_KPP, &
1871               RACM_ESRLSORG_AQCHEM_KPP,RACMSORG_KPP)
1872            call sum_pm_sorgam (                                             &
1873                 grid%alt, chem, grid%h2oaj, grid%h2oai,                                    &
1874                 grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1875                 grid%dust_opt,ids,ide, jds,jde, kds,kde,                    &
1876                 ims,ime, jms,jme, kms,kme,                                  &
1877                 its,ite, jts,jte, kts,kte-1                                 )
1878          case (RACM_SOA_VBS_KPP,RACM_SOA_VBS_AQCHEM_KPP)
1879              call sum_pm_soa_vbs (                                           &
1880                  grid%alt, chem, grid%h2oaj, grid%h2oai,                                   &
1881                  grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                &
1882                  config_flags%dust_opt,ids,ide, jds,jde, kds,kde,           &
1883                  ims,ime, jms,jme, kms,kme,                                 &
1884                  its,ite, jts,jte, kts,kte-1                                )
1885
1886         case (CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ, &
1887               CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, CBMZ_MOSAIC_DMS_4BIN_AQ,       &
1888               CBMZ_MOSAIC_DMS_8BIN_AQ, CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP)
1889            call sum_pm_mosaic (                                             &
1890                 grid%alt, chem,                                                  &
1891                 grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1892                 ids,ide, jds,jde, kds,kde,                                  &
1893                 ims,ime, jms,jme, kms,kme,                                  &
1894                 its,ite, jts,jte, kts,kte-1                                 )
1895
1896         case default
1897            do j=jts,min(jte,jde-1)
1898               do k=kts,min(kte,kde-1)
1899                  do i=its,min(ite,ide-1)
1900                     grid%pm2_5_dry(i,k,j)    = 0.
1901                     grid%pm2_5_water(i,k,j)  = 0.
1902                     grid%pm2_5_dry_ec(i,k,j) = 0.
1903                     grid%pm10(i,k,j)         = 0.
1904                     grid%tsoa(i,k,j)         = 0.
1905                  enddo
1906               enddo
1907            enddo
1908         end select
1909
1910         ! initialize advective tendency diagnostics for chem
1911         if ( grid%itimestep .eq. 0 .and. config_flags%chemdiag .eq. USECHEMDIAG ) then
1912           grid%advh_ct(:,:,:,:) = 0.
1913           grid%advz_ct(:,:,:,:) = 0.
1914         endif
1915 #endif
1916
1917         ! initialize advective tendency diagnostics for non-chem
1918         if ( grid%itimestep .eq. 0 .and. config_flags%tenddiag .eq. USETENDDIAG ) then
1919           advh_t(:,:,:,:) = 0.
1920           advz_t(:,:,:,:) = 0.
1921         endif
1922
1923    IF (num_chem >= PARAM_FIRST_SCALAR ) THEN
1924 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1925
1926       loop_3d_c   : DO loop = PARAM_FIRST_SCALAR , num_chem
1927          CALL set_physical_bc3d( chem(:,:,:,loop) , 'r' , config_flags ,                 &
1928                                  ids , ide , jds , jde , kds , kde ,        &
1929                                  ims , ime , jms , jme , kms , kme ,        &
1930                                  its , ite , jts , jte , kts , kte ,        &
1931                                  its , ite , jts , jte , kts , kte )
1932       END DO loop_3d_c
1933
1934    ENDIF
1935
1936    CALL set_physical_bc2d( grid%msftx , 'r' , config_flags ,              &
1937                          ids , ide , jds , jde , &
1938                          ims , ime , jms , jme , &
1939                          its , ite , jts , jte , &
1940                          its , ite , jts , jte   )
1941    CALL set_physical_bc2d( grid%msfty , 'r' , config_flags ,              &
1942                          ids , ide , jds , jde , &
1943                          ims , ime , jms , jme , &
1944                          its , ite , jts , jte , &
1945                          its , ite , jts , jte   )
1946    CALL set_physical_bc2d( grid%msfux , 'x' , config_flags ,              &
1947                          ids , ide , jds , jde , &
1948                          ims , ime , jms , jme , &
1949                          its , ite , jts , jte , &
1950                          its , ite , jts , jte   )
1951    CALL set_physical_bc2d( grid%msfuy , 'x' , config_flags ,              &
1952                          ids , ide , jds , jde , &
1953                          ims , ime , jms , jme , &
1954                          its , ite , jts , jte , &
1955                          its , ite , jts , jte   )
1956    CALL set_physical_bc2d( grid%msfvx , 'y' , config_flags ,              &
1957                          ids , ide , jds , jde , &
1958                          ims , ime , jms , jme , &
1959                          its , ite , jts , jte , &
1960                          its , ite , jts , jte   )
1961    CALL set_physical_bc2d( grid%msfvy , 'y' , config_flags ,              &
1962                          ids , ide , jds , jde , &
1963                          ims , ime , jms , jme , &
1964                          its , ite , jts , jte , &
1965                          its , ite , jts , jte   )
1966    CALL set_physical_bc2d( grid%sina , 'r' , config_flags ,              &
1967                          ids , ide , jds , jde , &
1968                          ims , ime , jms , jme , &
1969                          its , ite , jts , jte , &
1970                          its , ite , jts , jte   )
1971    CALL set_physical_bc2d( grid%cosa , 'r' , config_flags ,              &
1972                          ids , ide , jds , jde , &
1973                          ims , ime , jms , jme , &
1974                          its , ite , jts , jte , &
1975                          its , ite , jts , jte   )
1976    CALL set_physical_bc2d( grid%e , 'r' , config_flags ,                 &
1977                          ids , ide , jds , jde , &
1978                          ims , ime , jms , jme , &
1979                          its , ite , jts , jte , &
1980                          its , ite , jts , jte   )
1981    CALL set_physical_bc2d( grid%f , 'r' , config_flags ,                 &
1982                          ids , ide , jds , jde , &
1983                          ims , ime , jms , jme , &
1984                          its , ite , jts , jte , &
1985                          its , ite , jts , jte   )
1986
1987 #if (WRF_CHEM != 1)
1988       DEALLOCATE(CLDFRA_OLD)
1989 #endif
1990 #ifdef DM_PARALLEL
1991 # include "HALO_EM_INIT_1.inc"
1992 # include "HALO_EM_INIT_2.inc"
1993 # include "HALO_EM_INIT_3.inc"
1994 # include "HALO_EM_INIT_4.inc"
1995 # include "HALO_EM_INIT_5.inc"
1996 # include "PERIOD_BDY_EM_INIT.inc"
1997 # include "PERIOD_BDY_EM_MOIST.inc"
1998 # include "PERIOD_BDY_EM_TKE.inc"
1999 # include "PERIOD_BDY_EM_SCALAR.inc"
2000 # include "PERIOD_BDY_EM_CHEM.inc"
2001 #endif
2002
2003 DEALLOCATE(z_at_q)
2004
2005    IF (config_flags%p_lev_diags == PRESS_DIAGS ) THEN
2006     CALL wrf_debug ( 200 , ' PLD: pressure level diags' )
2007     CALL pld (                                  &
2008                  !  Input data for computing
2009                   U=grid%u_2                    &
2010                  ,V=grid%v_2                    &
2011                  ,W=grid%w_2                    &
2012                  ,t=grid%t_2                    &
2013                  ,qv=moist(:,:,:,P_QV)          &
2014                  ,zp=grid%ph_2                  &
2015                  ,zb=grid%phb                   &
2016                  ,pp=grid%p                     &
2017                  ,pb=grid%pb                    &
2018                  ,p=grid%p_hyd                  &
2019                  ,pw=grid%p_hyd_w               &
2020                  !  Map factors, coriolis for diags
2021                  ,msfux=grid%msfux              &
2022                  ,msfuy=grid%msfuy              &
2023                  ,msfvx=grid%msfvx              &
2024                  ,msfvy=grid%msfvy              &
2025                  ,msftx=grid%msftx              &
2026                  ,msfty=grid%msfty              &
2027                  ,f=grid%f                      &
2028                  ,e=grid%e                      &
2029                  !  Namelist info
2030                  ,use_tot_or_hyd_p=config_flags%use_tot_or_hyd_p &
2031                  ,extrap_below_grnd=config_flags%extrap_below_grnd &
2032                  ,missing=config_flags%p_lev_missing &
2033                  !  The diagnostics, mostly output variables
2034                  ,num_press_levels=config_flags%num_press_levels &
2035                  ,max_press_levels=max_plevs    &
2036                  ,press_levels=model_config_rec%press_levels &
2037                  ,p_pl  = grid%p_pl             &
2038                  ,u_pl  = grid%u_pl             &
2039                  ,v_pl  = grid%v_pl             &
2040                  ,t_pl  = grid%t_pl             &
2041                  ,rh_pl = grid%rh_pl            &
2042                  ,ght_pl= grid%ght_pl           &
2043                  ,s_pl  = grid%s_pl             &
2044                  ,td_pl = grid%td_pl            &
2045                  ,q_pl = grid%q_pl              &
2046                  !  Dimension arguments
2047                  ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde    &
2048                  ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme    &
2049                  ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte    )
2050    ENDIF
2051
2052    IF (config_flags%z_lev_diags == Z_DIAGS ) THEN
2053      CALL wrf_debug ( 200 , ' ZLD: height level and AGL diags' )
2054      CALL zld (                                  &
2055                  !  Input data for computing
2056                   U=grid%u_2                    &
2057                  ,V=grid%v_2                    &
2058                  ,W=grid%w_2                    &
2059                  ,t=grid%t_2                    &
2060                  ,qv=moist(:,:,:,P_QV)          &
2061                  ,zp=grid%ph_2                  &
2062                  ,zb=grid%phb                   &
2063                  ,pp=grid%p                     &
2064                  ,pb=grid%pb                    &
2065                  ,p=grid%p_hyd                  &
2066                  ,pw=grid%p_hyd_w               &
2067                  !  Map factors, coriolis for diags
2068                  ,msfux=grid%msfux              &
2069                  ,msfuy=grid%msfuy              &
2070                  ,msfvx=grid%msfvx              &
2071                  ,msfvy=grid%msfvy              &
2072                  ,msftx=grid%msftx              &
2073                  ,msfty=grid%msfty              &
2074                  ,f=grid%f                      &
2075                  ,e=grid%e                      &
2076                  ,ht=grid%ht                    &
2077                  !  Namelist info
2078                  ,use_tot_or_hyd_p=config_flags%use_tot_or_hyd_p &
2079                  ,extrap_below_grnd=config_flags%extrap_below_grnd &
2080                  ,missing=config_flags%z_lev_missing &
2081                  !  The diagnostics, mostly output variables
2082                  ,num_z_levels=config_flags%num_z_levels &
2083                  ,max_z_levels=max_zlevs    &
2084                  ,z_levels=model_config_rec%z_levels &
2085                  ,z_zl  = grid%z_zl             &
2086                  ,u_zl  = grid%u_zl             &
2087                  ,v_zl  = grid%v_zl             &
2088                  ,t_zl  = grid%t_zl             &
2089                  ,rh_zl = grid%rh_zl            &
2090                  ,ght_zl= grid%ght_zl           &
2091                  ,s_zl  = grid%s_zl             &
2092                  ,td_zl = grid%td_zl            &
2093                  ,q_zl = grid%q_zl              &
2094                  ,p_zl = grid%p_zl              &
2095                  !  Dimension arguments
2096                  ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde    &
2097                  ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme    &
2098                  ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte    )
2099    ENDIF
2100
2101    IF ( config_flags%diag_nwp2 == DO_TRAD_FIELDS ) THEN
2102       CALL wrf_debug ( 100 , '--> CALL DIAGNOSTICS PACKAGE: TRAD_FIELDS' )
2103       DO j = jts, MIN(jte,jde-1)
2104          DO k = 1, kte-1
2105             DO i = its, MIN(ite,ide-1)
2106                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
2107                grid%rho(i,k,j) = 1./( grid%al(i,k,j)+grid%alb(i,k,j) )
2108             END DO
2109          END DO
2110       END DO
2111       CALL trad_fields ( &
2112                !  Input data for computing
2113                U=grid%u_2                                           &
2114                ,V=grid%v_2                                          &
2115                ,W=grid%w_2                                          &
2116                ,t=grid%th_phy_m_t0                                  &
2117                ,qv=moist(:,:,:,P_QV)                                &
2118                ,zp=grid%ph_2                                        &
2119                ,zb=grid%phb                                         &
2120                ,pp=grid%p                                           &
2121                ,pb=grid%pb                                          &
2122                ,p=grid%p_hyd                                        &
2123                ,pw=grid%p_hyd_w                                     &
2124                !  Map factors, coriolis for diags
2125                ,msfux=grid%msfux                                    &
2126                ,msfuy=grid%msfuy                                    &
2127                ,msfvx=grid%msfvx                                    &
2128                ,msfvy=grid%msfvy                                    &
2129                ,msftx=grid%msftx                                    &
2130                ,msfty=grid%msfty                                    &
2131                ,f=grid%f                                            &
2132                ,e=grid%e                                            &
2133                ,sina=grid%sina                                      &
2134                ,cosa=grid%cosa                                      &
2135                !  Terrestrial data
2136                ,rho=grid%rho                                        &
2137                ,dz8w=dz8w                                           &
2138                ,ht=grid%ht                                          &
2139                !  Namelist info
2140                ,use_theta_m=config_flags%use_theta_m                &
2141                ,qc=moist(:,:,:,P_QC)                                &
2142                ,rainc=grid%rainc                                    &
2143                ,rainnc=grid%rainnc                                  &
2144                ,snownc=grid%snownc                                  &
2145                ,graupelnc=grid%graupelnc                            &
2146                ,hailnc=grid%hailnc                                  &
2147                !  The diagnostic output variables
2148                ,sealevelp=grid%sealevelp                            &
2149                ,temperature=grid%temperature                        &
2150                ,pressure=grid%pressure                              &
2151                ,geoheight=grid%geoheight                            &
2152                ,umet=grid%umet                                      &
2153                ,vmet=grid%vmet                                      &
2154                ,speed=grid%speed                                    &
2155                ,dir=grid%dir                                        &
2156                ,psfc=grid%psfc                                      &
2157                ,rain=grid%rain                                      &
2158                ,liqrain=grid%liqrain                                &
2159                ,tpw=grid%tpw                                        &
2160                ,potential_t=grid%potential_t                        &
2161                ,rh=grid%rh                                          &
2162                !  Various indexes for declarations, loop bounds
2163                ,IDS=ids,IDE=ide, JDS=jds,JDE=jde, KDS=kds,KDE=kde   &
2164                ,IMS=ims,IME=ime, JMS=jms,JME=jme, KMS=kms,KME=kme   &
2165                ,IPS=ips,IPE=ipe, JPS=jps,JPE=jpe, KPS=kps,KPE=kpe   &
2166                ,ITS=its,ITE=ite, JTS=jts,JTE=jte, KTS=kts,KTE=kte   )
2167    ENDIF
2168
2169
2170 !  FIRE
2171 if(config_flags%ifire.eq.2)then
2172
2173    call fire_driver_em_init ( grid , config_flags    &
2174             ,ids,ide, kds,kde, jds,jde                &
2175             ,ims,ime, kms,kme, jms,jme                &
2176             ,ips,ipe, kps,kpe, jps,jpe )
2177
2178    CALL wrf_debug ( 100 , 'start_domain_em: After call to fire_driver_em_init' )
2179
2180 !-----------------------------------------------------------------------
2181 ! fire spotting (passive Lagrangian particle transport, tracks
2182 ! firebrand physics properties)
2183 !-----------------------------------------------------------------------
2184
2185    IF (config_flags%fs_firebrand_gen_lim > 0) THEN
2186
2187        IF (config_flags%max_dom == grid%id) THEN
2188
2189            WRITE (message,*) 'SPFire_init: In inner most grid_id: ', grid%id, ' of ', config_flags%max_dom
2190            CALL wrf_debug ( 100 ,  message)
2191            CALL wrf_debug ( 100 , 'start_em: calling firebrand_spotting_em_init ...' )
2192
2193            CALL firebrand_spotting_em_init ( &
2194                grid          = grid,                &
2195                cf            = config_flags,        &
2196                fs_p_id       = grid%fs_p_id,        &
2197                fs_p_src      = grid%fs_p_src,       &
2198                fs_p_dt       = grid%fs_p_dt,        &
2199                fs_p_x        = grid%fs_p_x,         &
2200                fs_p_y        = grid%fs_p_y,         &
2201                fs_p_z        = grid%fs_p_z,         &
2202                fs_p_mass     = grid%fs_p_mass,      &
2203                fs_p_diam     = grid%fs_p_diam,      &
2204                fs_p_effd     = grid%fs_p_effd,      &
2205                fs_p_temp     = grid%fs_p_temp,      &
2206                fs_p_tvel     = grid%fs_p_tvel,      &
2207                fs_last_gen_dt= grid%fs_last_gen_dt, &
2208                fs_gen_idmax  = grid%fs_gen_idmax,   &
2209                fs_count_landed_all   = grid%fs_count_landed_all, &
2210                fs_count_landed_hist  = grid%fs_count_landed_hist,&
2211                fs_landing_mask       = grid%fs_landing_mask,     &
2212                fs_spotting_lkhd      = grid%fs_spotting_lkhd,    &
2213                fs_frac_landed        = grid%fs_frac_landed,      &
2214                fs_count_reset        = grid%fs_count_reset)
2215
2216            CALL wrf_debug ( 100 , 'start_em: firebrand_spotting_em_init ok ' )
2217        ENDIF
2218
2219    ENDIF
2220 !-----------------------------------------------------------------------
2221
2222 endif
2223
2224 #if ( WRFPLUS != 1 )
2225 if( grid%traj_opt /= no_trajectory ) then
2226      call trajectory_init( grid, config_flags, &
2227                            ims,ime, jms,jme, kms,kme )
2228      CALL wrf_debug ( 100 , 'start_domain_em: After call to trajectory_init' )
2229 endif
2230 #endif
2231
2232
2233      CALL wrf_debug ( 100 , 'start_domain_em: Returning' )
2234
2235      RETURN
2236
2237    END SUBROUTINE start_domain_em
2238
2239 !-------------------------------------------------------------------
2240
2241    SUBROUTINE rebalance_driver_cycl ( grid )
2242
2243       USE module_domain, ONLY : domain
2244       IMPLICIT NONE
2245
2246       TYPE (domain)          :: grid
2247
2248       CALL rebalance_cycl( grid &
2249 !
2250 #include "actual_new_args.inc"
2251 !
2252       )
2253
2254    END SUBROUTINE rebalance_driver_cycl
2255
2256 !---------------------------------------------------------------------
2257
2258    SUBROUTINE rebalance_cycl ( grid  &
2259 !
2260 #include "dummy_new_args.inc"
2261 !
2262                         )
2263       USE module_domain, ONLY : domain
2264       USE module_configure, ONLY : grid_config_rec_type, model_config_rec
2265       USE module_model_constants
2266       USE module_state_description
2267       USE module_driver_constants, ONLY: DATA_ORDER_XYZ, DATA_ORDER_YXZ, DATA_ORDER_ZXY, &
2268                                          DATA_ORDER_ZYX, DATA_ORDER_XZY, DATA_ORDER_YZX, &
2269                                          DATA_ORDER_XY, DATA_ORDER_YX, model_data_order
2270
2271 #ifdef DM_PARALLEL
2272       USE module_comm_dm, ONLY : &
2273                            HALO_EM_INIT_1_sub   &
2274                           ,HALO_EM_INIT_2_sub   &
2275                           ,HALO_EM_INIT_3_sub   &
2276                           ,HALO_EM_INIT_4_sub   &
2277                           ,HALO_EM_INIT_5_sub   &
2278                           ,HALO_EM_VINTERP_UV_1_sub
2279 #endif
2280 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
2281    USE module_dm, ONLY : ntasks_x, ntasks_y, ntasks, mytask, local_communicator
2282 #endif
2283
2284       IMPLICIT NONE
2285
2286       TYPE (domain)          :: grid
2287
2288 #include "dummy_new_decl.inc"
2289
2290       TYPE (grid_config_rec_type)              :: config_flags
2291
2292       REAL :: p_surf ,  pd_surf, p_surf_int , pb_int , ht_hold
2293       REAL :: qvf , qvf1 , qvf2, qtot
2294       REAL :: pfu, pfd, phm
2295       REAL :: z0, z1, z2, w1, w2
2296
2297       !  Local domain indices and counters.
2298
2299       INTEGER :: n_moist
2300
2301       INTEGER                             ::                       &
2302                                      ids, ide, jds, jde, kds, kde, &
2303                                      ims, ime, jms, jme, kms, kme, &
2304                                      its, ite, jts, jte, kts, kte, &
2305                                      ips, ipe, jps, jpe, kps, kpe, &
2306                                      i, j, k, kk, ispe, ktf
2307
2308       SELECT CASE ( model_data_order )
2309          CASE ( DATA_ORDER_ZXY )
2310             kds = grid%sd31 ; kde = grid%ed31 ;
2311             ids = grid%sd32 ; ide = grid%ed32 ;
2312             jds = grid%sd33 ; jde = grid%ed33 ;
2313
2314             kms = grid%sm31 ; kme = grid%em31 ;
2315             ims = grid%sm32 ; ime = grid%em32 ;
2316             jms = grid%sm33 ; jme = grid%em33 ;
2317
2318             kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
2319             its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
2320             jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
2321
2322          CASE ( DATA_ORDER_XYZ )
2323             ids = grid%sd31 ; ide = grid%ed31 ;
2324             jds = grid%sd32 ; jde = grid%ed32 ;
2325             kds = grid%sd33 ; kde = grid%ed33 ;
2326
2327             ims = grid%sm31 ; ime = grid%em31 ;
2328             jms = grid%sm32 ; jme = grid%em32 ;
2329             kms = grid%sm33 ; kme = grid%em33 ;
2330
2331             its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
2332             jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
2333             kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
2334
2335          CASE ( DATA_ORDER_XZY )
2336             ids = grid%sd31 ; ide = grid%ed31 ;
2337             kds = grid%sd32 ; kde = grid%ed32 ;
2338             jds = grid%sd33 ; jde = grid%ed33 ;
2339
2340             ims = grid%sm31 ; ime = grid%em31 ;
2341             kms = grid%sm32 ; kme = grid%em32 ;
2342             jms = grid%sm33 ; jme = grid%em33 ;
2343
2344             its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
2345             kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
2346             jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
2347
2348       END SELECT
2349
2350             ktf=MIN(kte,kde-1)
2351
2352       DO j=jts,jte
2353          DO i=its,ite
2354             grid%ph_2(i,1,j) = 0.
2355          END DO
2356       END DO
2357
2358       !  Base state potential temperature and inverse density (alpha = 1/rho) from
2359       !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
2360       !  from equation of state.  The potential temperature is a perturbation from t0.
2361
2362          n_moist = num_moist
2363
2364       print *,'n_moist,PARAM_FIRST_SCALAR',n_moist,PARAM_FIRST_SCALAR
2365
2366        DO j = jts, MIN(jte,jde-1)
2367          DO i = its, MIN(ite,ide-1)
2368
2369         IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
2370
2371                !  Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2372                !  equation) down from the top to get the pressure perturbation.  First get the pressure
2373                !  perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2374
2375                kk = kte - 1
2376                k=kk+1
2377
2378                qtot = 0.
2379                DO ispe=PARAM_FIRST_SCALAR,n_moist
2380                  qtot = qtot + 0.5*(moist(i,kk,j,ispe)+moist(i,kk,j,ispe))
2381                ENDDO
2382                qvf2 = 1./(1.+qtot)
2383                qvf1 = qtot*qvf2
2384
2385                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
2386                IF ( grid%use_theta_m == 0) THEN
2387                   qvf = 1.+rvovrd*moist(i,kk,j,P_QV)
2388                ELSE
2389                   qvf = 1.
2390                ENDIF
2391                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf*         &
2392                       (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2393                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2394                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2395
2396
2397                !  Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2398                !  inverse density fields (total and perturbation).
2399
2400         DO kk=kte-2,kts,-1
2401              k = kk + 1
2402
2403              qtot = 0.
2404              DO ispe=PARAM_FIRST_SCALAR,n_moist
2405                qtot = qtot + 0.5*(  moist(i,kk  ,j,ispe) + moist(i,kk+1,j,ispe) )
2406              ENDDO
2407                qvf2 = 1./(1.+qtot)
2408                qvf1 = qtot*qvf2
2409                grid%p(i,kk,j) = grid%p(i,kk+1,j) - ((grid%c1f(k)*grid%Mu_2(i,j)) +       &
2410                                qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/qvf2/grid%rdn(kk+1)
2411                IF ( grid%use_theta_m == 0) THEN
2412                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2413                ELSE
2414                   qvf = 1.
2415                ENDIF
2416                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
2417                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2418                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2419                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2420         ENDDO
2421
2422                !  This is the hydrostatic equation used in the model after the
2423                !  small timesteps.  In
2424                !  the model, grid%al (inverse density) is computed from the
2425                !  geopotential.
2426
2427             IF (grid%hypsometric_opt == 1) THEN
2428                   DO kk  = 2,kte
2429                      k = kk - 1
2430                      grid%ph_2(i,kk,j) = grid%ph_2(i,kk-1,j) - &
2431                                    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) &
2432                                  + (grid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,kk-1,j) )
2433                      grid%ph0(i,kk,j) = grid%ph_2(i,kk,j) + grid%phb(i,kk,j)
2434                   END DO
2435             ELSE IF (grid%hypsometric_opt == 2) THEN
2436                 ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is
2437                 ! dry pressure.
2438                 ! Note that al*p approximates Rd*T and dLOG(p) does z.
2439                 ! Here T varies mostly linear with z, the first-order
2440                 ! integration produces better result.
2441
2442                   grid%ph_2(i,1,j) = grid%phb(i,1,j)
2443                   DO k = 2,kte
2444                      pfu = grid%c3f(k  )*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k  ) + grid%p_top
2445                      pfd = grid%c3f(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k-1) + grid%p_top
2446                      phm = grid%c3h(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4h(k-1) + grid%p_top
2447                      grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
2448                   END DO
2449                   DO k = 1,kte
2450                      grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
2451                   END DO
2452                DO k = 1,kte
2453                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2454                END DO
2455             END IF ! hypsometric option
2456        ELSE  ! n_moist
2457          kk = kte - 1
2458                k = kk + 1
2459                qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk,j,P_QV))
2460                qvf2 = 1./(1.+qvf1)
2461                qvf1 = qvf1*qvf2
2462                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
2463                IF ( grid%use_theta_m == 0) THEN
2464                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2465                ELSE
2466                   qvf = 1.
2467                ENDIF
2468                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf    &
2469                            *(((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2470                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2471                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2472             !  Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2473             !  inverse density fields (total and perturbation).
2474            DO kk=kte-2,kts,-1
2475                k = kk + 1
2476                qvf1 = 0.5*(moist(i,kk,j,P_QV)+moist(i,kk+1,j,P_QV))
2477                qvf2 = 1./(1.+qvf1)
2478                qvf1 = qvf1*qvf2
2479                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)
2480                IF ( grid%use_theta_m == 0) THEN
2481                   qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
2482                ELSE
2483                   qvf = 1.
2484                ENDIF
2485                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
2486                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
2487                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
2488                grid%p_hyd(i,kk,j) = grid%p(i,kk,j) + grid%pb(i,kk,j)
2489            ENDDO
2490                !  This is the hydrostatic equation used in the model after the small timesteps.  In
2491                !  the model, grid%al (inverse density) is computed from the geopotential.
2492             IF (grid%hypsometric_opt == 1) THEN
2493                DO k  = 2,kte
2494                   grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
2495                                 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) &
2496                               + (grid%c1h(k)*grid%mu_2(i,j))*grid%alb(i,k-1,j) )
2497                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2498                END DO
2499             ELSE IF (grid%hypsometric_opt == 2) THEN
2500              ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry
2501              ! pressure.
2502              ! Note that al*p approximates Rd*T and dLOG(p) does z.
2503              ! Here T varies mostly linear with z, the first-order integration
2504              ! produces better result.
2505                grid%ph_2(i,1,j) = grid%phb(i,1,j)
2506                DO k = 2,kte
2507                   pfu = grid%c3f(k  )*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k  ) + grid%p_top
2508                   pfd = grid%c3f(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4f(k-1) + grid%p_top
2509                   phm = grid%c3h(k-1)*(grid%MUB(i,j)+grid%MU_2(i,j)) + grid%c4h(k-1) + grid%p_top
2510                   grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
2511                END DO
2512
2513                DO k = 1,kte
2514                   grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
2515                END DO
2516
2517                DO k = 1,kte
2518                   grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
2519                END DO
2520
2521             END IF ! hypsometric
2522        ENDIF ! nmoist
2523
2524 ! update surface pressure PSFC (needed for post-processing):
2525                z0 = grid%ph0(i,1,j)/g
2526                z1 = 0.5*(grid%ph0(i,1,j)+grid%ph0(i,2,j))/g
2527                z2 = 0.5*(grid%ph0(i,2,j)+grid%ph0(i,3,j))/g
2528                w1 = (z0 - z2)/(z1 - z2)
2529                w2 = 1. - w1
2530                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))
2531
2532          END DO !i
2533         ENDDO !j
2534
2535       ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2536 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
2537 # include "HALO_EM_INIT_1.inc"
2538 # include "HALO_EM_INIT_2.inc"
2539 # include "HALO_EM_INIT_3.inc"
2540 # include "HALO_EM_INIT_4.inc"
2541 # include "HALO_EM_INIT_5.inc"
2542 #endif
2543    END SUBROUTINE rebalance_cycl
2544
2545 !---------------------------------------------------------------------