wrftladj/start_em_ad.F

   1 !-------------------------------------------------------------------
   2
   3    SUBROUTINE a_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_model_constants
  13    USE a_module_bc, ONLY : a_set_physical_bc2d
  14    USE a_module_bc_em
  15    USE module_configure, ONLY : grid_config_rec_type
  16    USE module_tiles, ONLY : set_tiles
  17 #ifdef DM_PARALLEL
  18    USE module_dm, ONLY : wrf_dm_min_real, wrf_dm_max_real, wrf_dm_maxval, &
  19         ntasks_x, ntasks_y, &
  20         local_communicator_periodic, local_communicator, mytask, ntasks
  21 #else
  22    USE module_dm, ONLY : wrf_dm_min_real,  wrf_dm_max_real
  23 #endif
  24    USE module_comm_dm
  25
  26    USE a_module_physics_init
  27    USE module_diag_pld, ONLY : pld
  28
  29    USE module_model_constants
  30
  31    IMPLICIT NONE
  32    !  Input data.
  33    TYPE (domain)          :: grid
  34
  35    LOGICAL , INTENT(IN)   :: allowed_to_read
  36
  37    !  Definitions of dummy arguments to this routine (generated from Registry).
  38 # include "dummy_new_decl.inc"
  39
  40    !  Structure that contains run-time configuration (namelist) data for domain
  41    TYPE (grid_config_rec_type)              :: config_flags
  42
  43    !  Local data
  44    INTEGER                             ::                       &
  45                                   ids, ide, jds, jde, kds, kde, &
  46                                   ims, ime, jms, jme, kms, kme, &
  47                                   ips, ipe, jps, jpe, kps, kpe, &
  48                                   its, ite, jts, jte, kts, kte, &
  49                                   ij,i,j,k,ii,jj,kk,loop,error,l
  50
  51    INTEGER ::              imsx, imex, jmsx, jmex, kmsx, kmex,    &
  52                            ipsx, ipex, jpsx, jpex, kpsx, kpex,    &
  53                            imsy, imey, jmsy, jmey, kmsy, kmey,    &
  54                            ipsy, ipey, jpsy, jpey, kpsy, kpey
  55
  56    INTEGER     :: i_m
  57
  58    REAL        :: p00, t00, a, tiso, p_surf, pd_surf, temp, tiso_tmp
  59    REAL        :: p_strat, a_strat
  60    REAL :: qvf1, qvf2, qvf, a_qvf
  61    REAL :: pfu, pfd, phm, a_pfu, a_pfd, a_phm
  62    REAL :: MPDT
  63    REAL :: spongeweight
  64    LOGICAL :: first_trip_for_this_domain, start_of_simulation, fill_w_flag
  65    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
  66
  67 #if (WRF_CHEM!=1)
  68       REAL,ALLOCATABLE,DIMENSION(:,:,:) :: cldfra_old
  69 #endif
  70
  71    REAL :: lat1 , lat2 , lat3 , lat4
  72    REAL :: lon1 , lon2 , lon3 , lon4
  73    INTEGER :: num_points_lat_lon , iloc , jloc
  74    CHARACTER (LEN=132) :: message
  75    TYPE(WRFU_TimeInterval) :: stepTime
  76    REAL, DIMENSION(:,:), ALLOCATABLE :: clat_glob
  77    logical :: f_flux  ! flag for computing averaged fluxes in cu_gd
  78
  79    INTEGER :: idex, jdex
  80    INTEGER :: im1,ip1,jm1,jp1
  81    REAL :: temp2, temp1, temp0, temp0b, tempb, tempb0
  82    CHARACTER (LEN=256) :: timestr
  83
  84    REAL :: w_max, w_min
  85    LOGICAL :: w_needs_to_be_set
  86
  87    CALL get_ijk_from_grid ( grid ,                              &
  88                            ids, ide, jds, jde, kds, kde,        &
  89                            ims, ime, jms, jme, kms, kme,        &
  90                            ips, ipe, jps, jpe, kps, kpe,        &
  91                            imsx, imex, jmsx, jmex, kmsx, kmex,  &
  92                            ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
  93                            imsy, imey, jmsy, jmey, kmsy, kmey,  &
  94                            ipsy, ipey, jpsy, jpey, kpsy, kpey )
  95
  96    kts = kps ; kte = kpe     ! note that tile is entire patch
  97    its = ips ; ite = ipe     ! note that tile is entire patch
  98    jts = jps ; jte = jpe    ! note that tile is entire patch
  99 #if (WRF_CHEM!=1)
 100          ALLOCATE(CLDFRA_OLD(IMS:IME,KMS:KME,JMS:JME),STAT=I)  ; CLDFRA_OLD = 0.
 101 #endif
 102    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 103
 104 !kludge - need to stop CG from resetting precip and phys tendencies to zero
 105 !         when we are in here due to a nest being spawned, we want to still
 106 !         recompute the base state, but that is about it
 107    !  This is temporary and will need to be changed when grid%itimestep is removed.
 108
 109    IF ( grid%itimestep .EQ. 0 ) THEN
 110       first_trip_for_this_domain = .TRUE.
 111    ELSE
 112       first_trip_for_this_domain = .FALSE.
 113    END IF
 114
 115    IF ( config_flags%use_baseparam_fr_nml ) then
 116       CALL nl_get_base_pres  ( 1 , p00 )
 117       CALL nl_get_base_temp  ( 1 , t00 )
 118       CALL nl_get_base_lapse ( 1 , a   )
 119       CALL nl_get_iso_temp   ( 1 , tiso )
 120       CALL nl_get_base_lapse_strat ( 1 , a_strat )
 121       CALL nl_get_base_pres_strat  ( 1 , p_strat )
 122       IF ( ( t00 .LT. 100. .or. p00 .LT. 10000.) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 123          WRITE(wrf_err_message,*) 'start_em: BAD BASE STATE for T00 or P00 in namelist.input file'
 124          CALL wrf_error_fatal(TRIM(wrf_err_message))
 125       END IF
 126
 127    ELSE
 128    ! get these constants from model data
 129
 130       t00  = grid%t00
 131       p00  = grid%p00
 132       a    = grid%tlp
 133       tiso = grid%tiso
 134       a_strat = grid%tlp_strat
 135       p_strat = grid%p_strat
 136       IF ( ( t00 .LT. 100. .or. p00 .LT. 10000.) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 137          WRITE(wrf_err_message,*)&
 138          'start_em: did not find base state parameters in wrfinput. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
 139          CALL wrf_error_fatal(TRIM(wrf_err_message))
 140       ENDIF
 141
 142    ENDIF
 143
 144 !  check if tiso in the input file agrees with namelist value
 145
 146    CALL nl_get_iso_temp   ( 1 , tiso_tmp )
 147    IF ( ( tiso_tmp .NE. tiso ) .AND. ( .NOT. grid%this_is_an_ideal_run ) ) THEN
 148       WRITE(wrf_err_message,*)&
 149       'start_em: namelist iso_temp is not equal to iso_temp in wrfinput. Reset nml value and rerun'
 150       CALL wrf_error_fatal(TRIM(wrf_err_message))
 151    ENDIF
 152
 153    IF ( .NOT. config_flags%restart .AND. &
 154         (( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ))) THEN
 155
 156       IF ( config_flags%map_proj .EQ. 0 ) THEN
 157          CALL wrf_error_fatal ( 'start_domain: Idealized case cannot have a separate nested input file' )
 158       END IF
 159
 160       !  Base state potential temperature and inverse density (alpha = 1/rho) from
 161       !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
 162       !  from equation of state.  The potential temperature is a perturbation from t0.
 163
 164       DO j = jts, MIN(jte,jde-1)
 165          DO i = its, MIN(ite,ide-1)
 166
 167             !  Base state pressure is a function of eta level and terrain, only, plus
 168             !  the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
 169             !  temperature, K), A (temperature difference, from 1000 mb to 300 mb, K),
 170             !  tiso (isothermal temperature at tropopause/lower stratosphere),
 171             !  p_strat (pressure at top of isothermal layer), A_strat (lapse rate in
 172             !  stratosphere above isothermal layer)
 173
 174             p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
 175
 176             DO k = 1, kte-1
 177                grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
 178                temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
 179                IF ( grid%pb(i,k,j) .LT. p_strat ) THEN
 180                   temp = tiso + A_strat * LOG ( grid%pb(i,k,j)/p_strat )
 181                ENDIF
 182                grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
 183 !              grid%t_init(i,k,j) = (t00 + A*LOG(grid%pb(i,k,j)/p00))*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
 184                grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 185             END DO
 186
 187             !  Base state mu is defined as base state surface pressure minus grid%p_top
 188
 189             grid%mub(i,j) = p_surf - grid%p_top
 190
 191             !  Integrate base geopotential, starting at terrain elevation.  This assures that
 192             !  the base state is in exact hydrostatic balance with respect to the model equations.
 193             !  This field is on full levels.
 194
 195             grid%phb(i,1,j) = grid%ht(i,j) * g
 196
 197             IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 198                DO k  = 2,kte
 199                   grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
 200                END DO
 201             ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 202                DO k  = 2,kte
 203                   pfu = grid%mub(i,j)*grid%znw(k)   + grid%p_top
 204                   pfd = grid%mub(i,j)*grid%znw(k-1) + grid%p_top
 205                   phm = grid%mub(i,j)*grid%znu(k-1) + grid%p_top
 206                   grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
 207                END DO
 208             END IF
 209          END DO
 210       END DO
 211
 212    ENDIF
 213
 214 #ifdef DM_PARALLEL
 215 #   include "HALO_EM_INIT_5_AD.inc"
 216 #   include "HALO_EM_INIT_4_AD.inc"
 217 #   include "HALO_EM_INIT_3_AD.inc"
 218 #   include "HALO_EM_INIT_2_AD.inc"
 219 #   include "HALO_EM_INIT_1_AD.inc"
 220 #endif
 221
 222         ! initialize advective tendency diagnostics for non-chem
 223         if ( grid%itimestep .eq. 0 .and. config_flags%tenddiag .eq. USETENDDIAG ) then
 224           a_advh_t(:,:,:,:) = 0.
 225           a_advz_t(:,:,:,:) = 0.
 226         endif
 227
 228    IF (num_scalar > 0) THEN
 229
 230 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
 231
 232       loop_3d_s   : DO loop = 1 , num_scalar
 233          CALL a_set_physical_bc3d( a_scalar(:,:,:,loop) , 'r' , config_flags ,               &
 234                                  ids , ide , jds , jde , kds , kde ,        &
 235                                  ims , ime , jms , jme , kms , kme ,        &
 236                                  its , ite , jts , jte , kts , kte ,        &
 237                                  its , ite , jts , jte , kts , kte )
 238       END DO loop_3d_s
 239
 240    ENDIF
 241
 242    IF (num_moist > 0) THEN
 243
 244 ! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
 245
 246       loop_3d_m   : DO loop = 1 , num_moist
 247          CALL a_set_physical_bc3d( a_moist(:,:,:,loop) , 'r' , config_flags ,                 &
 248                                  ids , ide , jds , jde , kds , kde ,        &
 249                                  ims , ime , jms , jme , kms , kme ,        &
 250                                  its , ite , jts , jte , kts , kte ,        &
 251                                  its , ite , jts , jte , kts , kte )
 252       END DO loop_3d_m
 253
 254    ENDIF
 255
 256    CALL a_set_physical_bc3d( grid%a_w_1 , 'W' , config_flags ,                 &
 257                          ids , ide , jds , jde , kds , kde ,        &
 258                          ims , ime , jms , jme , kms , kme ,        &
 259                          its , ite , jts , jte , kts , kte ,        &
 260                          its , ite , jts , jte , kts , kte )
 261    CALL a_set_physical_bc3d( grid%a_w_2 , 'W' , config_flags ,                 &
 262                          ids , ide , jds , jde , kds , kde ,        &
 263                          ims , ime , jms , jme , kms , kme ,        &
 264                          its , ite , jts , jte , kts , kte ,        &
 265                          its , ite , jts , jte , kts , kte )
 266
 267    CALL a_set_physical_bc3d( grid%a_ph_1 , 'W' , config_flags ,                 &
 268                          ids , ide , jds , jde , kds , kde ,        &
 269                          ims , ime , jms , jme , kms , kme ,        &
 270                          its , ite , jts , jte , kts , kte ,        &
 271                          its , ite , jts , jte , kts , kte )
 272
 273    CALL a_set_physical_bc3d( grid%a_ph_2 , 'W' , config_flags ,                 &
 274                          ids , ide , jds , jde , kds , kde ,        &
 275                          ims , ime , jms , jme , kms , kme ,        &
 276                          its , ite , jts , jte , kts , kte ,        &
 277                          its , ite , jts , jte , kts , kte )
 278
 279    CALL a_set_physical_bc3d( grid%a_t_1 , 't' , config_flags ,                 &
 280                          ids , ide , jds , jde , kds , kde ,        &
 281                          ims , ime , jms , jme , kms , kme ,        &
 282                          its , ite , jts , jte , kts , kte ,        &
 283                          its , ite , jts , jte , kts , kte )
 284
 285    CALL a_set_physical_bc3d( grid%a_t_2 , 't' , config_flags ,                 &
 286                          ids , ide , jds , jde , kds , kde ,        &
 287                          ims , ime , jms , jme , kms , kme ,        &
 288                          its , ite , jts , jte , kts , kte ,        &
 289                          its , ite , jts , jte , kts , kte )
 290
 291    CALL a_set_physical_bc2d( grid%a_mu_1, 't' , config_flags ,   &
 292                            ids , ide , jds , jde ,  &
 293                            ims , ime , jms , jme ,  &
 294                            its , ite , jts , jte ,  &
 295                            its , ite , jts , jte   )
 296    CALL a_set_physical_bc2d( grid%a_mu_2, 't' , config_flags ,   &
 297                            ids , ide , jds , jde ,  &
 298                            ims , ime , jms , jme ,  &
 299                            its , ite , jts , jte ,  &
 300                            its , ite , jts , jte   )
 301    CALL a_set_physical_bc3d( grid%a_php , 'W' , config_flags ,                 &
 302                          ids , ide , jds , jde , kds , kde ,        &
 303                          ims , ime , jms , jme , kms , kme ,        &
 304                          its , ite , jts , jte , kts , kte ,        &
 305                          its , ite , jts , jte , kts , kte )
 306
 307    CALL a_set_physical_bc3d( grid%a_al , 't' , config_flags ,                 &
 308                          ids , ide , jds , jde , kds , kde ,        &
 309                          ims , ime , jms , jme , kms , kme ,        &
 310                          its , ite , jts , jte , kts , kte ,        &
 311                          its , ite , jts , jte , kts , kte )
 312    CALL a_set_physical_bc3d( grid%a_alt , 't' , config_flags ,                 &
 313                          ids , ide , jds , jde , kds , kde ,        &
 314                          ims , ime , jms , jme , kms , kme ,        &
 315                          its , ite , jts , jte , kts , kte ,        &
 316                          its , ite , jts , jte , kts , kte )
 317    CALL a_set_physical_bc3d( grid%a_tke_2, 't' , config_flags ,                 &
 318                          ids , ide , jds , jde , kds , kde ,        &
 319                          ims , ime , jms , jme , kms , kme ,        &
 320                          its , ite , jts , jte , kts , kte ,        &
 321                          its , ite , jts , jte , kts , kte )
 322
 323    IF ( first_trip_for_this_domain ) THEN
 324
 325 ! set GMT outside of phy_init because phy_init may not be called on this
 326 ! process if, for example, it is a moving nest and if this part of the domain is not
 327 ! being initialized (not the leading edge).
 328    CALL domain_setgmtetc( grid, start_of_simulation )
 329 !tgs
 330    IF(config_flags%cycling) start_of_simulation = .true.
 331 !  print *,'cycling, start_of_simulation -->',config_flags%cycling, start_of_simulation
 332
 333 ! Set the w at the surface.  If this is the start of a forecast, or if this is a cycled run
 334 ! and the start of that forecast, we define the w field.  However, if this a restart, then
 335 ! we leave w alone.  Initial value is V dot grad(topo) at surface, then smoothly decaying
 336 ! above that.
 337    END IF
 338
 339    IF ( ( .NOT. start_of_simulation .OR. config_flags%cycling ) .AND. ( .NOT. config_flags%restart ) ) THEN
 340
 341 ! If W already exists (not zero), then we leave it alone.  How to do this?  We find the
 342 ! max/min on this node at the surface.  If parallel, we collect the max/min from all procs.
 343 ! If the max/min throughout the entire domain at the surface is identically 0, then we say
 344 ! that the W field is NOT initialized, and we run the set_w_surface routines for the
 345 ! two time levels of W.  If the field is already initialized, we do NOT run those two
 346 ! routines.
 347
 348       w_max = grid%w_2(its,1,jts)
 349       w_min = grid%w_2(its,1,jts)
 350       DO j = jts, MIN(jte,jde-1)
 351          DO i = its, MIN(ite,ide-1)
 352          w_max = MAX ( w_max , grid%w_2(i,1,j) )
 353          w_min = MIN ( w_min , grid%w_2(i,1,j) )
 354          END DO
 355       END DO
 356 #ifdef DM_PARALLEL
 357       w_max = wrf_dm_max_real ( w_max )
 358       w_min = wrf_dm_min_real ( w_min )
 359 #endif
 360
 361       IF ( ( ABS(w_max) .LT. 1.E-6 ) .AND. &
 362            ( ABS(w_min) .LT. 1.E-6 ) ) THEN
 363          w_needs_to_be_set = .TRUE.
 364       ELSE
 365          IF ( config_flags%use_input_w ) THEN
 366             w_needs_to_be_set = .FALSE.
 367          ELSE
 368             w_needs_to_be_set = .TRUE.
 369          END IF
 370       END IF
 371
 372       IF ( w_needs_to_be_set ) THEN
 373
 374          fill_w_flag = .true.
 375          CALL a_set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
 376                               grid%w_1, grid%a_w_1, grid%ht, &
 377                               grid%u_1, grid%a_u_1, grid%v_1, grid%a_v_1, grid%cf1, &
 378                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
 379                               ids, ide, jds, jde, kds, kde,                    &
 380                               ims, ime, jms, jme, kms, kme,                    &
 381                               its, ite, jts, jte, kts, kte                     )
 382          CALL a_set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
 383                               grid%w_2, grid%a_w_2, grid%ht,  &
 384                               grid%u_2, grid%a_u_2, grid%v_2, grid%a_v_2, grid%cf1, &
 385                               grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
 386                               ids, ide, jds, jde, kds, kde,                    &
 387                               ims, ime, jms, jme, kms, kme,                    &
 388                               its, ite, jts, jte, kts, kte                     )
 389       ENDIF
 390
 391    ENDIF
 392
 393 ! finished setting kinematic condition for w at the surface
 394
 395    CALL a_set_physical_bc3d( grid%a_u_1 , 'U' , config_flags ,                  &
 396                          ids , ide , jds , jde , kds , kde ,        &
 397                          ims , ime , jms , jme , kms , kme ,        &
 398                          its , ite , jts , jte , kts , kte ,        &
 399                          its , ite , jts , jte , kts , kte )
 400    CALL a_set_physical_bc3d( grid%a_u_2 , 'U' , config_flags ,                  &
 401                          ids , ide , jds , jde , kds , kde ,        &
 402                          ims , ime , jms , jme , kms , kme ,        &
 403                          its , ite , jts , jte , kts , kte ,        &
 404                          its , ite , jts , jte , kts , kte )
 405
 406    CALL a_set_physical_bc3d( grid%a_v_1 , 'V' , config_flags ,                  &
 407                          ids , ide , jds , jde , kds , kde ,        &
 408                          ims , ime , jms , jme , kms , kme ,        &
 409                          its , ite , jts , jte , kts , kte ,        &
 410                          its , ite , jts , jte , kts , kte )
 411    CALL a_set_physical_bc3d( grid%a_v_2 , 'V' , config_flags ,                  &
 412                          ids , ide , jds , jde , kds , kde ,        &
 413                          ims , ime , jms , jme , kms , kme ,        &
 414                          its , ite , jts , jte , kts , kte ,        &
 415                          its , ite , jts , jte , kts , kte )
 416
 417  !  set physical boundary conditions for all initialized variables
 418
 419 !-----------------------------------------------------------------------
 420 !  Stencils for patch communications  (WCS, 29 June 2001)
 421 !  Note:  the size of this halo exchange reflects the
 422 !         fact that we are carrying the uncoupled variables
 423 !         as state variables in the mass coordinate model, as
 424 !         opposed to the coupled variables as in the height
 425 !         coordinate model.
 426 !
 427 !                           * * * * *
 428 !         *        * * *    * * * * *
 429 !       * + *      * + *    * * + * *
 430 !         *        * * *    * * * * *
 431 !                           * * * * *
 432 !
 433 !j  grid%u_1                          x
 434 !j  grid%u_2                          x
 435 !j  grid%v_1                          x
 436 !j  grid%v_2                          x
 437 !j  grid%w_1                          x
 438 !j  grid%w_2                          x
 439 !j  grid%t_1                          x
 440 !j  grid%t_2                          x
 441 !j grid%ph_1                          x
 442 !j grid%ph_2                          x
 443 !
 444 !j  grid%t_init                       x
 445 !
 446 !j  grid%phb   x
 447 !j  grid%ph0   x
 448 !j  grid%php   x
 449 !j   grid%pb   x
 450 !j   grid%al   x
 451 !j  grid%alt   x
 452 !j  grid%alb   x
 453 !
 454 !  the following are 2D (xy) variables
 455 !
 456 !j grid%mu_1                          x
 457 !j grid%mu_2                          x
 458 !j grid%mub   x
 459 !j grid%mu0   x
 460 !j grid%ht    x
 461 !j grid%msftx x
 462 !j grid%msfty x
 463 !j grid%msfux x
 464 !j grid%msfuy x
 465 !j grid%msfvx x
 466 !j grid%msfvy x
 467 !j grid%sina  x
 468 !j grid%cosa  x
 469 !j grid%e     x
 470 !j grid%f     x
 471 !
 472 !  4D variables
 473 !
 474 ! moist                        x
 475 !  chem                        x
 476 !scalar                        x
 477
 478 !--------------------------------------------------------------
 479
 480 #ifdef DM_PARALLEL
 481 #  include "HALO_EM_INIT_5_AD.inc"
 482 #  include "HALO_EM_INIT_4_AD.inc"
 483 #  include "HALO_EM_INIT_3_AD.inc"
 484 #  include "HALO_EM_INIT_2_AD.inc"
 485 #  include "HALO_EM_INIT_1_AD.inc"
 486 #endif
 487
 488    IF ( first_trip_for_this_domain ) THEN
 489
 490    CALL wrf_debug ( 100 , 'start_domain_em_ad: Before call to a_phy_init' )
 491
 492 ! namelist MPDT does not exist yet, so set it here
 493 ! MPDT is the call frequency for microphysics in minutes (0 means every step)
 494    MPDT = 0.
 495
 496    CALL set_tiles ( grid , grid%imask_nostag, ims, ime, jms, jme, ips, ipe, jps, jpe )
 497 !
 498 ! Phy init can do reads and broadcasts when initializing physics -- landuse for example. However, if
 499 ! we're running on a reduced mesh (that is, some tasks don't have any work) we have to at least let them
 500 ! pass through this code so the broadcasts don't hang on the other, active tasks.  Set the number of
 501 ! tiles to a minimum of 1 and assume that the backwards patch ranges (ips=0, ipe=-1) will prevent
 502 ! anything else from happening on the blank tasks.  JM 20080605
 503 !
 504    if ( allowed_to_read ) grid%num_tiles = max(1,grid%num_tiles)
 505 !
 506 ! Phy_init is not necessarily thread-safe; do not multi-thread this loop.
 507 ! The tiling is to handle the fact that we may be masking off part of the computation.
 508 !
 509    DO ij = 1, grid%num_tiles
 510
 511 !tgs do not need physics initialization for backward DFI integration
 512     IF ( ( grid%dfi_opt .EQ. DFI_NODFI ) .or. &
 513           ( ( grid%dfi_stage .NE. DFI_BCK ) .and. &
 514           ( grid%dfi_stage .NE. DFI_STARTBCK ) ) ) THEN    !tgs, mods by tah
 515
 516      CALL a_phy_init( grid%id , config_flags, grid%DT, grid%RESTART, grid%znw, grid%znu,   &
 517                       grid%p_top, grid%tsk, grid%RADT,grid%BLDT,grid%CUDT, MPDT, &
 518                       grid%rucuten, grid%rvcuten, grid%rthcuten, grid%a_rthcuten, &
 519                       grid%rqvcuten, grid%a_rqvcuten, grid%rqrcuten, grid%rqccuten, &
 520                       grid%a_rqccuten, &
 521                       grid%rqscuten, grid%rqicuten,                          &
 522                       grid%rushten, grid%rvshten, grid%rthshten,             &
 523                       grid%rqvshten, grid%rqrshten, grid%rqcshten,           &
 524                       grid%rqsshten, grid%rqishten, grid%rqgshten,           &
 525                       grid%rublten,grid%a_rublten,grid%rvblten,grid%a_rvblten, &
 526                       grid%rthblten, grid%a_rthblten,              &
 527                       grid%rqvblten,grid%a_rqvblten,grid%rqcblten,grid%rqiblten, &
 528                       grid%rthraten,grid%rthratenlw,grid%rthratensw,         &
 529                       grid%stepbl,grid%stepra,grid%stepcu,                   &
 530                       grid%w0avg, grid%a_w0avg, grid%rainnc, grid%a_rainnc, &
 531                       grid%rainc, grid%a_rainc, grid%raincv, grid%a_raincv, &
 532                       grid%rainncv, grid%a_rainncv, &
 533                       grid%snownc, grid%a_snownc, grid%snowncv, grid%a_snowncv, &
 534                       grid%graupelnc, grid%a_graupelnc, grid%graupelncv, grid%a_graupelncv, &
 535                       grid%nca,grid%swrad_scat,                    &
 536                       grid%cldefi,grid%lowlyr,                          &
 537                       grid%mass_flux,                              &
 538                       grid%rthften, grid%rqvften,                       &
 539                       grid%cldfra,                                      &
 540 #if (WRF_CHEM!=1)
 541                       cldfra_old,                                  &
 542 #endif
 543                       grid%glw,grid%gsw,grid%emiss,grid%embck,            &
 544                       grid%lu_index,                                      &
 545                       grid%landuse_ISICE, grid%landuse_LUCATS,            &
 546                       grid%landuse_LUSEAS, grid%landuse_ISN,              &
 547                       grid%lu_state,                                      &
 548                       grid%xlat,grid%xlong,grid%albedo,grid%albbck,grid%GMT,grid%JULYR,grid%JULDAY,     &
 549                       grid%levsiz, num_ozmixm, num_aerosolc, grid%paerlev,  &
 550                       grid%alevsiz, grid%no_src_types,                      &
 551                       grid%tmn,grid%xland,grid%znt,grid%z0,grid%ust,grid%mol,grid%pblh,grid%tke_pbl,    &
 552                       grid%exch_h,grid%thc,grid%snowc,grid%mavail,grid%hfx,grid%qfx,grid%rainbl, &
 553                       grid%tslb,grid%zs,grid%dzs,config_flags%num_soil_layers,grid%warm_rain,  &
 554                       grid%adv_moist_cond,                         &
 555                       grid%apr_gr,grid%apr_w,grid%apr_mc,grid%apr_st,grid%apr_as,      &
 556                       grid%apr_capma,grid%apr_capme,grid%apr_capmi,          &
 557                       grid%xice,grid%xicem,grid%vegfra,grid%snow,grid%canwat,grid%smstav,         &
 558                       grid%smstot, grid%sfcrunoff,grid%udrunoff,grid%grdflx,grid%acsnow,      &
 559                       grid%acsnom,grid%ivgtyp,grid%isltyp, grid%sfcevp,grid%smois,     &
 560                       grid%sh2o, grid%snowh, grid%smfr3d,                    &
 561                       grid%snoalb,                 &
 562                       grid%DX,grid%DY,grid%f_ice_phy,grid%f_rain_phy,grid%f_rimef_phy, &
 563                       grid%mp_restart_state,grid%tbpvs_state,grid%tbpvs0_state,&
 564                       allowed_to_read, grid%moved, start_of_simulation,               &
 565                       grid%LAGDAY, &
 566                       ids, ide, jds, jde, kds, kde,           &
 567                       ims, ime, jms, jme, kms, kme,           &
 568                       grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, &
 569                       config_flags%num_urban_ndm,                           & !multi-layer urban
 570                       config_flags%urban_map_zrd,                           & !multi-layer urban
 571                       config_flags%urban_map_zwd,                           & !multi-layer urban
 572                       config_flags%urban_map_gd,                            & !multi-layer urban
 573                       config_flags%urban_map_zd,                            & !multi-layer urban
 574                       config_flags%urban_map_zdf,                           & !multi-layer urban
 575                       config_flags%urban_map_bd,                            & !multi-layer urban
 576                       config_flags%urban_map_wd,                            & !multi-layer urban
 577                       config_flags%urban_map_gbd,                           & !multi-layer urban
 578                       config_flags%urban_map_fbd,                           & !multi-layer urban
 579                       config_flags%num_urban_hi,                            & !multi-layer urban
 580                       grid%raincv_a,grid%raincv_b,                                    &
 581                       grid%gd_cloud, grid%gd_cloud2,                                  & ! Optional
 582                       grid%gd_cloud_a, grid%gd_cloud2_a,                              & ! Optional
 583                       grid%QC_CU, grid%QI_CU,                                         & ! Optional
 584                       ozmixm,grid%pin,                             &     ! Optional
 585                       grid%m_ps_1,grid%m_ps_2,grid%m_hybi,aerosolc_1,aerosolc_2,&  ! Optional
 586                       grid%rundgdten,grid%rvndgdten,grid%rthndgdten,         &     ! Optional
 587                       grid%rphndgdten,grid%rqvndgdten,grid%rmundgdten,       &     ! Optional
 588                       grid%FGDT,grid%stepfg,                        &     ! Optional
 589                       grid%cugd_tten,grid%cugd_ttens,grid%cugd_qvten,        &   ! Optional
 590                       grid%cugd_qvtens,grid%cugd_qcten,                 &   ! Optional
 591                       grid%ISNOWXY, grid%ZSNSOXY, grid%TSNOXY,                                  &   ! Optional Noah-MP
 592                       grid%SNICEXY, grid%SNLIQXY, grid%TVXY, grid%TGXY, grid%CANICEXY,          &   ! Optional Noah-MP
 593                       grid%CANLIQXY, grid%EAHXY, grid%TAHXY, grid%CMXY,                         &   ! Optional Noah-MP
 594                       grid%CHXY, grid%FWETXY, grid%SNEQVOXY,grid%ALBOLDXY,grid%QSNOWXY,grid%QRAINXY, &   ! Optional Noah-MP
 595                       grid%WSLAKEXY, grid%ZWTXY, grid%WAXY, grid%WTXY, grid%LFMASSXY, grid%RTMASSXY, & ! Optional Noah-MP
 596                       grid%STMASSXY, grid%WOODXY, grid%STBLCPXY, grid%FASTCPXY,                 &   ! Optional Noah-MP
 597                       grid%XSAIXY,                                                              &   ! Optional Noah-MP
 598                       grid%T2MVXY, grid%T2MBXY, grid%CHSTARXY,                                  &   ! Optional Noah-MP
 599                       grid%DZR, grid%DZB, grid%DZG,                          & !Optional urban
 600                       grid%TR_URB2D,grid%TB_URB2D,grid%TG_URB2D,grid%TC_URB2D,    & !Optional urban
 601                       grid%QC_URB2D, grid%XXXR_URB2D,grid%XXXB_URB2D,        & !Optional urban
 602                       grid%XXXG_URB2D, grid%XXXC_URB2D,                 & !Optional urban
 603                       grid%TRL_URB3D, grid%TBL_URB3D, grid%TGL_URB3D,        & !Optional urban
 604                       grid%SH_URB2D, grid%LH_URB2D, grid%G_URB2D, grid%RN_URB2D,  & !Optional urban
 605                       grid%TS_URB2D, grid%FRC_URB2D, grid%UTYPE_URB2D,      & !Optional urban
 606                       grid%TRB_URB4D,grid%TW1_URB4D,grid%TW2_URB4D,grid%TGB_URB4D,grid%TLEV_URB3D,  & !multi-layer urban
 607                       grid%QLEV_URB3D,grid%TW1LEV_URB3D,grid%TW2LEV_URB3D,        & !multi-layer urban
 608                       grid%TGLEV_URB3D,grid%TFLEV_URB3D,grid%SF_AC_URB3D,         & !multi-layer urban
 609                       grid%LF_AC_URB3D,grid%CM_AC_URB3D,grid%SFVENT_URB3D,grid%LFVENT_URB3D,   & !multi-layer urban
 610                       grid%SFWIN1_URB3D,grid%SFWIN2_URB3D, & !multi-layer urban
 611                       grid%SFW1_URB3D,grid%SFW2_URB3D,grid%SFR_URB3D,grid%SFG_URB3D,     & !multi-layer urban
 612                       grid%A_U_BEP,grid%A_V_BEP,grid%A_T_BEP,grid%A_Q_BEP,             & !multi-layer urban
 613                       grid%A_E_BEP,grid%B_U_BEP,grid%B_V_BEP,grid%B_T_BEP,             & !multi-layer urban
 614                       grid%B_Q_BEP,grid%B_E_BEP,grid%DLG_BEP,                          & !multi-layer urban
 615                       grid%DL_U_BEP,grid%SF_BEP,grid%VL_BEP,                           & !multi-layer urban
 616                       grid%TML,grid%T0ML,grid%HML,grid%H0ML,grid%HUML,grid%HVML,grid%TMOML,     & !Optional oml
 617                       grid%itimestep, grid%fdob,            &
 618                       t00, p00, a,                      & ! for obs_nudge base state
 619                       grid%TYR, grid%TYRA, grid%TDLY, grid%TLAG, grid%NYEAR, grid%NDAY,grid%tmn_update, &
 620                       grid%achfx, grid%aclhf, grid%acgrdflx                 &
 621                       ,grid%te_temf,grid%cf3d_temf,grid%wm_temf        & ! WA
 622                       ,grid%massflux_EDKF, grid%entr_EDKF, grid%detr_EDKF                &
 623                       ,grid%thl_up,grid%thv_up,grid%rt_up                                &
 624                       ,grid%rv_up,grid%rc_up,grid%u_up,grid%v_up,grid%frac_up            &
 625                       )
 626        ENDIF   !tgs
 627
 628    ENDDO
 629
 630    CALL wrf_debug ( 100 , 'start_domain_em_ad: After call to a_phy_init' )
 631
 632    END IF
 633
 634 !
 635 !-------------
 636
 637    IF(.not.config_flags%restart)THEN
 638
 639     DO j=jts,min(jte,jde-1)
 640        DO i=its,min(ite,ide-1)
 641           DO k=kte, kts+1, -1
 642              grid%a_p(i,k-1,j) = grid%a_p(i,k-1,j) + 2.*grid%a_p_hyd_w(i,k,j)
 643              grid%a_p_hyd_w(i,k-1,j) = grid%a_p_hyd_w(i,k-1,j) - grid%a_p_hyd_w(i,k,j)
 644              grid%a_p_hyd_w(i,k,j) = 0.0
 645           ENDDO
 646           grid%a_p(i,1,j) = grid%a_p(i,1,j) + grid%a_p_hyd_w(i,1,j)
 647           grid%a_p_hyd_w(i,1,j) = 0.0
 648        ENDDO
 649     ENDDO
 650
 651 !  reconstitute base-state fields
 652
 653      IF ( first_trip_for_this_domain ) THEN
 654        DO j = jts,min(jte,jde-1)
 655        DO k = kts,kte-1
 656        DO i = its, min(ite,ide-1)
 657            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 658              grid%t_1(i,k,j)=grid%t_2(i,k,j)
 659            ENDIF
 660        ENDDO
 661        ENDDO
 662        ENDDO
 663
 664        DO j = jts,min(jte,jde-1)
 665        DO k = kts,kte
 666        DO i = its, min(ite,ide-1)
 667           grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
 668        ENDDO
 669        ENDDO
 670        ENDDO
 671
 672        DO j = jts,min(jte,jde-1)
 673        DO i = its, min(ite,ide-1)
 674            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 675              grid%mu_1(i,j)=grid%mu_2(i,j)
 676            ENDIF
 677        ENDDO
 678        ENDDO
 679      END IF
 680
 681     IF(config_flags%max_dom .EQ. 1)THEN
 682 ! with single domain, grid%t_init from wrfinput is OK to use
 683      DO j = jts,min(jte,jde-1)
 684      DO k = kts,kte-1
 685      DO i = its, min(ite,ide-1)
 686        IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 687          grid%pb(i,k,j) = grid%znu(k)*grid%mub(i,j)+grid%p_top
 688          grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 689        ENDIF
 690      ENDDO
 691      ENDDO
 692      ENDDO
 693     ENDIF
 694
 695 ! Use equations from calc_p_rho_phi to derive p and al from ph
 696
 697     IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 698        DO j=jts,min(jte,jde-1)
 699        DO k=kts,kte-1
 700        DO i=its,min(ite,ide-1)
 701           grid%al(i,k,j)=-1./(grid%mub(i,j)+grid%mu_1(i,j))*(grid%alb(i,k,j)*grid%mu_1(i,j)  &
 702                      +grid%rdnw(k)*(grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)))
 703        ENDDO
 704        ENDDO
 705        ENDDO
 706     ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 707        DO j=jts,min(jte,jde-1)
 708        DO k=kts,kte-1
 709        DO i=its,min(ite,ide-1)
 710           pfu = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znw(k+1)+grid%p_top
 711           pfd = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znw(k)  +grid%p_top
 712           phm = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znu(k)  +grid%p_top
 713           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)) &
 714                             /phm/LOG(pfd/pfu)-grid%alb(i,k,j)
 715        ENDDO
 716        ENDDO
 717        ENDDO
 718     END IF
 719
 720     DO j=min(jte,jde-1),jts,-1
 721     DO k=kte-1,kts,-1
 722     DO i=min(ite,ide-1),its,-1
 723         grid%a_p(i,k,j) = grid%a_p(i,k,j) + grid%a_p_hyd(i,k,j)
 724         IF ( .NOT. config_flags%var4d_run ) THEN
 725            qvf = 1. + rvovrd*grid%moist(i,k,j,P_QV)
 726            temp2 = p1000mb*(grid%alb(i,k,j)+grid%al(i,k,j))
 727            temp1 = t0 + grid%t_1(i,k,j)
 728            temp0 = temp1*qvf/temp2
 729            IF (r_d*temp0 .LE. 0.0 .AND. (cpovcv .EQ. 0.0 .OR. cpovcv .NE. &
 730                 INT(cpovcv))) THEN
 731              temp0b = 0.0
 732            ELSE
 733              temp0b = r_d*cpovcv*(r_d*temp0)**(cpovcv-1)*p1000mb*grid%a_p(i,k,j)/temp2
 734            END IF
 735            grid%a_t_1(i,k,j) = grid%a_t_1(i,k,j) + qvf*temp0b
 736            a_qvf = temp1*temp0b
 737            grid%a_p(i,k,j) = 0.0
 738            grid%a_moist(i,k,j,P_QV) = grid%a_moist(i,k,j,P_QV) + rvovrd*a_qvf
 739         END IF
 740         IF ( .NOT. config_flags%var4d_run ) THEN
 741            grid%a_al(i,k,j) = grid%a_al(i,k,j) + grid%a_alt(i,k,j) - temp0*p1000mb*temp0b
 742         ELSE
 743            grid%a_al(i,k,j) = grid%a_al(i,k,j) + grid%a_alt(i,k,j)
 744         END IF
 745         grid%a_alt(i,k,j) = 0.0
 746         grid%a_p_hyd(i,k,j) = 0.0
 747     ENDDO
 748     ENDDO
 749     ENDDO
 750
 751     IF ( config_flags%hypsometric_opt .EQ. 1 ) THEN
 752        DO j=min(jte,jde-1),jts,-1
 753        DO k=kte-1,kts,-1
 754        DO i=min(ite,ide-1),its,-1
 755           temp = grid%mub(i,j) + grid%mu_1(i,j)
 756           tempb = -(grid%a_al(i,k,j)/temp)
 757           grid%a_mu_1(i,j) = grid%a_mu_1(i,j) + (grid%alb(i,k,j)-(grid%alb(i,k,j)*grid%mu_1(i,j) &
 758              +grid%rdnw(k)*(grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)))/temp)*tempb
 759           grid%a_ph_1(i,k+1,j) = grid%a_ph_1(i,k+1,j) + grid%rdnw(k)*tempb
 760           grid%a_ph_1(i,k,j) = grid%a_ph_1(i,k,j) - grid%rdnw(k)*tempb
 761           grid%a_al(i,k,j) = 0.0
 762        ENDDO
 763        ENDDO
 764        ENDDO
 765     ELSE IF ( config_flags%hypsometric_opt .EQ. 2 ) THEN
 766        DO j=min(jte,jde-1),jts,-1
 767        DO k=kte-1,kts,-1
 768        DO i=min(ite,ide-1),its,-1
 769           pfu = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znw(k+1) + grid%p_top
 770           phm = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znu(k) + grid%p_top
 771           pfd = (grid%mub(i,j)+grid%mu_1(i,j))*grid%znw(k) + grid%p_top
 772           temp0 = pfd/pfu
 773           temp1 = LOG(temp0)
 774           temp = phm*temp1
 775           tempb = grid%a_al(i,k,j)/temp
 776           tempb0 = -((grid%phb(i,k+1,j)-grid%phb(i,k,j)+grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j))*tempb/temp)
 777           temp0b = phm*tempb0/(temp0*pfu)
 778           grid%a_ph_1(i,k+1,j) = grid%a_ph_1(i,k+1,j) + tempb
 779           grid%a_ph_1(i,k,j) = grid%a_ph_1(i,k,j) - tempb
 780           a_phm = temp1*tempb0
 781           a_pfd = temp0b
 782           a_pfu = -(temp0*temp0b)
 783           grid%a_al(i,k,j) = 0.0
 784           grid%a_mu_1(i,j) = grid%a_mu_1(i,j) + grid%znw(k)*a_pfd + grid%znw(k+1)*a_pfu + grid%znu(k)*a_phm
 785        ENDDO
 786        ENDDO
 787        ENDDO
 788     END IF
 789
 790 !  if this is for a nested domain, the defined/interpolated fields are the _2
 791
 792      IF ( first_trip_for_this_domain ) THEN
 793
 794 ! data that is expected to be zero must be explicitly initialized as such
 795 !    grid%h_diabatic = 0.
 796
 797        DO j = jts,min(jte,jde-1)
 798        DO k = kts,kte-1
 799        DO i = its, min(ite,ide-1)
 800            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 801              grid%a_t_2(i,k,j) = grid%a_t_2(i,k,j) + grid%a_t_1(i,k,j)
 802              grid%a_t_1(i,k,j) = 0.0
 803            ENDIF
 804        ENDDO
 805        ENDDO
 806        ENDDO
 807
 808        DO j = jts,min(jte,jde-1)
 809        DO k = kts,kte
 810        DO i = its, min(ite,ide-1)
 811           grid%a_ph_2(i,k,j) = grid%a_ph_2(i,k,j) + grid%a_ph_1(i,k,j)
 812           grid%a_ph_1(i,k,j) = 0.0
 813        ENDDO
 814        ENDDO
 815        ENDDO
 816
 817        DO j = jts,min(jte,jde-1)
 818        DO i = its, min(ite,ide-1)
 819            IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
 820              grid%a_mu_2(i,j) = grid%a_mu_2(i,j) + grid%a_mu_1(i,j)
 821              grid%a_mu_1(i,j) = 0.0
 822            ENDIF
 823        ENDDO
 824        ENDDO
 825      END IF
 826
 827    ENDIF
 828
 829
 830 #if (WRF_CHEM!=1)
 831       DEALLOCATE(CLDFRA_OLD)
 832 #endif
 833
 834      CALL wrf_debug ( 100 , 'start_domain_em_ad: Returning' )
 835
 836      RETURN
 837
 838    END SUBROUTINE a_start_domain_em
 839