dyn_em/module_initialize_fire.F

   1 !IDEAL:MODEL_LAYER:INITIALIZATION
   2 !
   3
   4 !  This MODULE holds the routines which are used to perform various initializations
   5 !  for the individual domains.
   6
   7 !  This MODULE CONTAINS the following routines:
   8
   9 !  initialize_field_test - 1. Set different fields to different constant
  10 !                             values.  This is only a test.  If the correct
  11 !                             domain is not found (based upon the "id")
  12 !                             then a fatal error is issued.
  13
  14 !-----------------------------------------------------------------------
  15
  16 MODULE module_initialize_ideal
  17
  18    USE module_domain
  19    USE module_io_domain
  20    USE module_state_description
  21    USE module_model_constants
  22    USE module_bc
  23    USE module_timing
  24    USE module_configure
  25    USE module_init_utilities
  26    USE module_soil_pre        !AK/ak for full surface initialization
  27 #ifdef DM_PARALLEL
  28    USE module_dm
  29 #endif
  30    USE module_fr_fire_util, ONLY: continue_at_boundary,crash,read_array_2d_real, &
  31      read_array_2d_integer,interpolate_2d,set_ideal_coord,print_2d_stats
  32    USE module_fr_fire_phys, ONLY: fuelmc_g,read_namelist_fire
  33
  34 CONTAINS
  35
  36
  37 !-------------------------------------------------------------------
  38 ! this is a wrapper for the solver-specific init_domain routines.
  39 ! Also dereferences the grid variables and passes them down as arguments.
  40 ! This is crucial, since the lower level routines may do message passing
  41 ! and this will get fouled up on machines that insist on passing down
  42 ! copies of assumed-shape arrays (by passing down as arguments, the
  43 ! data are treated as assumed-size -- ie. f77 -- arrays and the copying
  44 ! business is avoided).  Fie on the F90 designers.  Fie and a pox.
  45
  46    SUBROUTINE init_domain ( grid )
  47
  48    IMPLICIT NONE
  49
  50    !  Input data.
  51    TYPE (domain), POINTER :: grid
  52    !  Local data.
  53    INTEGER :: idum1, idum2
  54
  55    CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
  56
  57      CALL init_domain_rk( grid &
  58 !
  59 #include "actual_new_args.inc"
  60 !
  61                         )
  62
  63    END SUBROUTINE init_domain
  64
  65 !-------------------------------------------------------------------
  66
  67    SUBROUTINE init_domain_rk ( grid &
  68 !
  69 # include "dummy_new_args.inc"
  70 !
  71 )
  72    IMPLICIT NONE
  73
  74    !  Input data.
  75    TYPE (domain), POINTER :: grid
  76
  77 # include "dummy_new_decl.inc"
  78
  79    TYPE (grid_config_rec_type)              :: config_flags
  80
  81    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
  82
  83    !  Local data
  84    INTEGER                             ::                       &
  85                                   ids, ide, jds, jde, kds, kde, &
  86                                   ims, ime, jms, jme, kms, kme, &
  87                                   its, ite, jts, jte, kts, kte, &
  88                                   i, j, k
  89
  90    INTEGER, PARAMETER :: nl_max = 1000
  91    REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
  92    INTEGER :: nl_in
  93
  94    INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
  95    REAL    :: B1, B2, B3, B4, B5
  96    REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
  97    REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
  98    REAL    :: x_rad, y_rad, z_rad, hght_pert   !Ak/ak
  99    character (len=256) :: mminlu2              !AK/ak land use scheme (USGS)
 100 !   REAL, EXTERNAL :: interp_0
 101    REAL    :: hm
 102    REAL    :: pi
 103
 104 !  stuff from original initialization that has been dropped from the Registry
 105    REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
 106    REAL    :: qvf1, qvf2, pd_surf
 107    INTEGER :: it
 108    real :: thtmp, ptmp, temp(3)
 109
 110    LOGICAL :: moisture_init
 111    LOGICAL :: stretch_grd, dry_sounding
 112    LOGICAL :: stretch_hyp, sfc_init           !AK/ak switches for hyperbolic grid streching and surface initialization
 113
 114    INTEGER :: xs , xe , ys , ye
 115    INTEGER :: mtn_type
 116    INTEGER :: ks, ke, id
 117    INTEGER :: &  ! fire mesh sizes
 118               iots,iote,jots,jote, &            ! tile dims out
 119              ifds,ifde, kfds,kfde, jfds,jfde,                              &
 120              ifms,ifme, kfms,kfme, jfms,jfme,                              &
 121              ifts,ifte, kfts,kfte, jfts,jfte
 122
 123    REAL :: mtn_ht, mtn_max, mtn_x, mtn_y, mtn_z, grad_max
 124    REAL :: tign_max,tign_min
 125    REAL :: mtn_axs, mtn_ays, mtn_axe, mtn_aye
 126    REAL :: mtn_fxs, mtn_fys, mtn_fxe, mtn_fye
 127    REAL :: mtn_xs, mtn_ys, mtn_xe, mtn_ye
 128    REAL :: fdx,fdy ! fire mesh step
 129    INTEGER:: ir,jr ! refinement factors
 130    REAL :: minhfx,maxhfx,totheat
 131
 132    logical have_fire_ht,have_fire_grad,have_atm_grad
 133
 134 !*** executable
 135
 136    SELECT CASE ( model_data_order )
 137          CASE ( DATA_ORDER_ZXY )
 138    kds = grid%sd31 ; kde = grid%ed31 ;
 139    ids = grid%sd32 ; ide = grid%ed32 ;
 140    jds = grid%sd33 ; jde = grid%ed33 ;
 141
 142    kms = grid%sm31 ; kme = grid%em31 ;
 143    ims = grid%sm32 ; ime = grid%em32 ;
 144    jms = grid%sm33 ; jme = grid%em33 ;
 145
 146    kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
 147    its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
 148    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 149          CASE ( DATA_ORDER_XYZ )
 150    ids = grid%sd31 ; ide = grid%ed31 ;
 151    jds = grid%sd32 ; jde = grid%ed32 ;
 152    kds = grid%sd33 ; kde = grid%ed33 ;
 153
 154    ims = grid%sm31 ; ime = grid%em31 ;
 155    jms = grid%sm32 ; jme = grid%em32 ;
 156    kms = grid%sm33 ; kme = grid%em33 ;
 157
 158    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 159    jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
 160    kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
 161          CASE ( DATA_ORDER_XZY )
 162    ids = grid%sd31 ; ide = grid%ed31 ;
 163    kds = grid%sd32 ; kde = grid%ed32 ;
 164    jds = grid%sd33 ; jde = grid%ed33 ;
 165
 166    ims = grid%sm31 ; ime = grid%em31 ;
 167    kms = grid%sm32 ; kme = grid%em32 ;
 168    jms = grid%sm33 ; jme = grid%em33 ;
 169
 170    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 171    kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
 172    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 173
 174    END SELECT
 175
 176 !   z_scale = .40
 177    pi = 2.*asin(1.0)
 178    write(6,*) ' pi is ',pi
 179    nxc = (ide-ids)/2
 180    nyc = (jde-jds)/2
 181
 182    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 183
 184 ! here we check to see if the boundary conditions are set properly
 185
 186    CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
 187
 188     delt = config_flags%delt_perturbation
 189     x_rad = config_flags%xrad_perturbation
 190     y_rad = config_flags%yrad_perturbation
 191     z_rad = config_flags%zrad_perturbation
 192     hght_pert = config_flags%hght_perturbation
 193
 194     stretch_grd = config_flags%stretch_grd
 195     stretch_hyp = config_flags%stretch_hyp
 196         z_scale = config_flags%z_grd_scale
 197        sfc_init = config_flags%sfc_full_init
 198
 199  moisture_init = .true.   !AK/ak
 200
 201     grid%itimestep=0
 202
 203 #ifdef DM_PARALLEL
 204    CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
 205    CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
 206 #endif
 207
 208 !AK/ak land use initialization (USGS)
 209    IF (sfc_init) THEN
 210     mminlu2=' '
 211     mminlu2(1:4) = 'USGS'                !Ak/ak
 212     CALL nl_set_mminlu(1, mminlu2)       !Ak/ak
 213     CALL nl_get_iswater(1,grid%iswater) ! Ak/ak
 214    ENDIF
 215
 216     CALL nl_set_iswater(1,0)
 217     CALL nl_set_cen_lat(1,40.)
 218     CALL nl_set_cen_lon(1,-105.)
 219     CALL nl_set_truelat1(1,0.)
 220     CALL nl_set_truelat2(1,0.)
 221     CALL nl_set_moad_cen_lat (1,0.)
 222     CALL nl_set_stand_lon (1,0.)
 223     CALL nl_set_pole_lon (1,0.)
 224     CALL nl_set_pole_lat (1,90.)
 225     CALL nl_set_map_proj(1,0)
 226
 227 !  here we initialize data we currently is not initialized
 228 !  in the input data
 229
 230     DO j = jts, jte
 231       DO i = its, ite
 232          grid%msftx(i,j)    = 1.
 233          grid%msfty(i,j)    = 1.
 234          grid%msfux(i,j)    = 1.
 235          grid%msfuy(i,j)    = 1.
 236          grid%msfvx(i,j)    = 1.
 237          grid%msfvx_inv(i,j)= 1.
 238          grid%msfvy(i,j)    = 1.
 239          grid%sina(i,j)     = 0.
 240          grid%cosa(i,j)     = 1.
 241          grid%e(i,j)        = 0.
 242          grid%f(i,j)        = 0.
 243       END DO
 244     END DO
 245
 246 ! ***** fire
 247       write(6,*) '*************************************'
 248
 249 !AK/ak surface initialization latitude, longitude, landuse index from from LANDUSE.TBL skin temperature and soil temperature
 250    IF (sfc_init) THEN
 251     DO j = jts, jte
 252       DO i = its, ite
 253           grid%xlat(i,j) = config_flags%fire_lat_init     !Ak/sk (35)
 254          grid%xlong(i,j) = config_flags%fire_lon_init     !Ak/ak (-111)
 255          grid%xland(i,j) = 1.                             !Ak/ak
 256       grid%lu_index(i,j) = config_flags%sfc_lu_index      !AK/ak land use index (28)
 257            grid%tsk(i,j) = config_flags%sfc_tsk           !AK/ak  surface skin temperature [K] (285)
 258            grid%tmn(i,j) = config_flags%sfc_tmn           !AK/ak  soil temperature at lower boundary [K] (285)
 259       END DO
 260     END DO
 261       ! read land use data from files, overwriting the constant
 262       if(config_flags%fire_read_lu) &
 263           call read_array_2d_real('input_lu',grid%lu_index,ids,ide,jds,jde,ims,ime,jms,jme)
 264       if(config_flags%fire_read_tsk) &
 265           call read_array_2d_real   ('input_tsk',grid%tsk,    ids,ide,jds,jde,ims,ime,jms,jme)
 266       if(config_flags%fire_read_tmn) &
 267           call read_array_2d_real   ('input_tmn',grid%tmn,    ids,ide,jds,jde,ims,ime,jms,jme)
 268
 269 ! for Noah LSM, additional variables need to be initializedi  !AK/ak |
 270
 271   other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
 272
 273       CASE (SLABSCHEME)
 274       write(6,*) ' SLAB surface scheme activated'
 275
 276       CASE (LSMSCHEME)
 277       write(6,*) ' Noah unified LSM scheme activated with:'
 278       write(6,*) '    vegetation fraction=',config_flags%sfc_vegfra
 279       write(6,*) '          canopy  water=',config_flags%sfc_canwat
 280       write(6,*) '     dominant veg. type=',config_flags%sfc_ivgtyp
 281       write(6,*) '     dominant soil type=',config_flags%sfc_isltyp
 282
 283         DO j = jts , MIN(jde-1,jte)
 284            DO i = its , MIN(ide-1,ite)
 285                grid%vegfra(i,j) = config_flags%sfc_vegfra  !0.5
 286                grid%canwat(i,j) = config_flags%sfc_canwat  !0.
 287                grid%ivgtyp(i,j) = config_flags%sfc_ivgtyp  !18
 288                grid%isltyp(i,j) = config_flags%sfc_isltyp  !7
 289                grid%xice(i,j) = 0.
 290                grid%snow(i,j) = 0.
 291            END DO
 292         END DO
 293
 294       CASE (RUCLSMSCHEME)
 295        write(6,*) ' RUS surface scheme activated'
 296     END SELECT other_masked_fields                         !AK/ak |
 297
 298     ENDIF
 299
 300     DO j = jts, jte
 301     DO k = kts, kte
 302       DO i = its, ite
 303          grid%ww(i,k,j)     = 0.
 304       END DO
 305    END DO
 306    END DO
 307
 308    grid%step_number = 0
 309
 310    IF (sfc_init) THEN
 311
 312       write(6,*) ' full surface initialization activated '
 313       ! write(6,*) ' land use index =', config_flags%sfc_lu_index
 314       ! write(6,*) ' skin temperature=',grid%tsk(10,10),&
 315       !         '[K] soil temperature=', grid%tmn(10,10),'[K]'
 316 ! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F
 317       CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, &
 318                      grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, &
 319                      grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, &
 320                      model_config_rec%sf_surface_physics(grid%id), &
 321                                    ids,ide, jds,jde, kds,kde,&
 322                                    ims,ime, jms,jme, kms,kme,&
 323                                    its,ite, jts,jte, kts,kte )
 324
 325
 326      ELSE
 327       write(6,*) 'full surface initialization is turned off!! '
 328     ENDIF    !end of surface initialization
 329
 330 ! set up the grid
 331    write(6,*) '*************************************'
 332
 333    IF (model_config_rec%eta_levels(1) .EQ. -1) THEN !we do not have eta_levels from namelist
 334      IF (stretch_grd) THEN ! exponential or hyperbolic tangential stretch for eta
 335
 336       IF (stretch_hyp) THEN ! hyperbolic tangential stretch (more levels at the surface)
 337        write(6,*) ' hyperbolic tangential stretching activated with z_scale =',z_scale
 338        DO k=1, kde
 339         grid%znw(k) = -1.* (tanh(z_scale*(float(k-1) / float(kde-1) -1.)))/ &
 340                                    (tanh(z_scale))
 341        ENDDO
 342       ELSE                 ! exponential stretch for eta (nearly constant dz)
 343       write(6,*) ' exponential grid stretching activated with z_scale =',z_scale
 344        DO k=1, kde
 345        grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
 346                               (1.-exp(-1./z_scale))
 347        ENDDO
 348       ENDIF
 349      ELSE
 350      write(6,*) ' no grid stretching'
 351        DO k=1, kde
 352         grid%znw(k) = 1. - float(k-1)/float(kde-1)
 353        ENDDO
 354      ENDIF
 355    ELSE
 356      CALL wrf_debug(0,"module_initialize_les: vertical nesting is enabled, using eta_levels specified in namelist.input")
 357      ks = 0
 358      DO id=1,grid%id
 359         ks = ks+model_config_rec%e_vert(id)
 360      ENDDO
 361      IF (ks .GT. max_eta) THEN
 362         CALL wrf_error_fatal("too many vertical levels, increase max_eta in frame/module_driver_constants.F")
 363      ENDIF
 364         !Now set the eta_levels to what we specified in the namelist. We've
 365         !packed all the domains' eta_levels into a 'vector' and now we need
 366         !to pull only the section of the vector associated with our domain
 367         !of interest, which is between indicies ks and ke.
 368      IF (grid%id .EQ. 1) THEN
 369        ks = 1
 370        ke = model_config_rec%e_vert(1)
 371      ELSE
 372        id = 1
 373        ks = 1
 374        ke = 0
 375        DO WHILE (grid%id .GT. id)
 376          id = id+1
 377          ks = ks+model_config_rec%e_vert(id-1)
 378          ke = ks+model_config_rec%e_vert(id)
 379        ENDDO
 380      ENDIF
 381      DO k=1,kde
 382        grid%znw(k) = model_config_rec%eta_levels(ks+k-1)
 383      ENDDO
 384      !Check the value of the first and last eta level for our domain,
 385      !then check that the vector of eta levels is only decreasing
 386      IF (grid%znw(1) .NE. 1.0) THEN
 387        CALL wrf_error_fatal("error with specified eta_levels, first level is not 1.0")
 388      ENDIF
 389      IF (grid%znw(kde) .NE. 0.0) THEN
 390        CALL wrf_error_fatal("error with specified eta_levels, last level is not 0.0")
 391      ENDIF
 392      DO k=2,kde
 393        IF (grid%znw(k) .GT. grid%znw(k-1)) THEN
 394          CALL wrf_error_fatal("eta_levels are not uniformly decreasing from 1.0 to 0.0")
 395        ENDIF
 396      ENDDO
 397    ENDIF
 398    write(6,*) '*************************************'
 399
 400    DO k=1, kde-1
 401     grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
 402     grid%rdnw(k) = 1./grid%dnw(k)
 403     grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
 404    ENDDO
 405
 406    DO k=2, kde-1
 407     grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
 408     grid%rdn(k) = 1./grid%dn(k)
 409     grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
 410     grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
 411    ENDDO
 412
 413    cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
 414    cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
 415    grid%cf1  = grid%fnp(2) + cof1
 416    grid%cf2  = grid%fnm(2) - cof1 - cof2
 417    grid%cf3  = cof2
 418
 419    grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
 420    grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
 421    grid%rdx = 1./config_flags%dx
 422    grid%rdy = 1./config_flags%dy
 423
 424 !  get the sounding from the ascii sounding file, first get dry sounding and
 425 !  calculate base state
 426
 427   dry_sounding = .true.
 428   IF ( wrf_dm_on_monitor() ) THEN
 429   write(6,*) ' getting dry sounding for base state '
 430
 431   CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 432   ENDIF
 433   CALL wrf_dm_bcast_real( zk , nl_max )
 434   CALL wrf_dm_bcast_real( p_in , nl_max )
 435   CALL wrf_dm_bcast_real( pd_in , nl_max )
 436   CALL wrf_dm_bcast_real( theta , nl_max )
 437   CALL wrf_dm_bcast_real( rho , nl_max )
 438   CALL wrf_dm_bcast_real( u , nl_max )
 439   CALL wrf_dm_bcast_real( v , nl_max )
 440   CALL wrf_dm_bcast_real( qv , nl_max )
 441   CALL wrf_dm_bcast_integer ( nl_in , 1 )
 442
 443   write(6,*) ' returned from reading sounding, nl_in is ',nl_in
 444
 445 !  find ptop for the desired ztop (ztop is input from the namelist),
 446 !  and find surface pressure
 447
 448   grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
 449
 450   !  For hybrid coord
 451
 452   DO k=kts, kte
 453      IF      ( config_flags%hybrid_opt .EQ. 0 ) THEN
 454         grid%c3f(k) = grid%znw(k)
 455      ELSE IF ( config_flags%hybrid_opt .EQ. 1 ) THEN
 456         grid%c3f(k) = grid%znw(k)
 457      ELSE IF ( config_flags%hybrid_opt .EQ. 2 ) THEN
 458         B1 = 2. * grid%etac**2 * ( 1. - grid%etac )
 459         B2 = -grid%etac * ( 4. - 3. * grid%etac - grid%etac**3 )
 460         B3 = 2. * ( 1. - grid%etac**3 )
 461         B4 = - ( 1. - grid%etac**2 )
 462         B5 = (1.-grid%etac)**4
 463         grid%c3f(k) = ( B1 + B2*grid%znw(k) + B3*grid%znw(k)**2 + B4*grid%znw(k)**3 ) / B5
 464         IF ( grid%znw(k) .LT. grid%etac ) THEN
 465            grid%c3f(k) = 0.
 466         END IF
 467         IF ( k .EQ. kds ) THEN
 468            grid%c3f(k) = 1.
 469         ELSE IF ( k .EQ. kde ) THEN
 470            grid%c3f(k) = 0.
 471         END IF
 472      ELSE IF ( config_flags%hybrid_opt .EQ. 3 ) THEN
 473         grid%c3f(k) = grid%znw(k)*sin(0.5*3.14159*grid%znw(k))**2
 474         IF      ( k .EQ. kds ) THEN
 475            grid%c3f(k) = 1.
 476         ELSE IF ( k .EQ. kds ) THEN
 477            grid%c3f(kde) = 0.
 478         END IF
 479      ELSE
 480         CALL wrf_message     ( 'ERROR: --- hybrid_opt' )
 481         CALL wrf_message     ( 'ERROR: --- hybrid_opt=0    ==> Standard WRF terrain-following coordinate' )
 482         CALL wrf_message     ( 'ERROR: --- hybrid_opt=1    ==> Standard WRF terrain-following coordinate, hybrid c1, c2, c3, c4' )
 483         CALL wrf_message     ( 'ERROR: --- hybrid_opt=2    ==> Hybrid, Klemp polynomial' )
 484         CALL wrf_message     ( 'ERROR: --- hybrid_opt=3    ==> Hybrid, sin^2' )
 485         CALL wrf_error_fatal ( 'ERROR: --- Invalid option' )
 486      END IF
 487   END DO
 488
 489   !  c4 is a function of c3 and eta.
 490
 491   DO k=1, kde
 492      grid%c4f(k) = ( grid%znw(k) - grid%c3f(k) ) * ( p1000mb - grid%p_top )
 493   ENDDO
 494
 495   !  Now on half levels, just add up and divide by 2 (for c3h).  Use (eta-c3)*(p00-pt) for c4 on half levels.
 496
 497   DO k=1, kde-1
 498      grid%znu(k) = ( grid%znw(k+1) + grid%znw(k) ) * 0.5
 499      grid%c3h(k) = ( grid%c3f(k+1) + grid%c3f(k) ) * 0.5
 500      grid%c4h(k) = ( grid%znu(k) - grid%c3h(k) ) * ( p1000mb - grid%p_top )
 501   ENDDO
 502
 503   !  c1 = d(B)/d(eta).  We define c1f as c1 on FULL levels.  For a vertical difference,
 504   !  we need to use B and eta on half levels.  The k-loop ends up referring to the
 505   !  full levels, neglecting the top and bottom.
 506
 507   DO k=kds+1, kde-1
 508      grid%c1f(k) = ( grid%c3h(k) - grid%c3h(k-1) ) / ( grid%znu(k) - grid%znu(k-1) )
 509   ENDDO
 510
 511   !  The boundary conditions to get the coefficients:
 512   !  1) At k=kts: define d(B)/d(eta) = 1.  This gives us the same value of B and d(B)/d(eta)
 513   !     when doing the sigma-only B=eta.
 514   !  2) At k=kte: define d(B)/d(eta) = 0.  The curve B SMOOTHLY goes to zero, and at the very
 515   !     top, B continues to SMOOTHLY go to zero.  Note that for almost all cases of non B=eta,
 516   !     B is ALREADY=ZERO at the top, so this is a reasonable BC to assume.
 517
 518   grid%c1f(kds) = 1.
 519   IF ( ( config_flags%hybrid_opt .EQ. 0 ) .OR. ( config_flags%hybrid_opt .EQ. 1 ) ) THEN
 520      grid%c1f(kde) = 1.
 521   ELSE
 522      grid%c1f(kde) = 0.
 523   END IF
 524
 525   !  c2 = ( 1. - c1(k) ) * (p00 - pt).  There is no vertical differencing, so we can do the
 526   !  full kds to kde looping.
 527
 528   DO k=kds, kde
 529      grid%c2f(k) = ( 1. - grid%c1f(k) ) * ( p1000mb - grid%p_top )
 530   ENDDO
 531
 532   !  Now on half levels for c1 and c2.  The c1h will result from the full level c3 and full
 533   !  level eta differences.  The c2 value use the half level c1(k).
 534
 535   DO k=1, kde-1
 536      grid%c1h(k) = ( grid%c3f(k+1) - grid%c3f(k) ) / ( grid%znw(k+1) - grid%znw(k) )
 537      grid%c2h(k) = ( 1. - grid%c1h(k) ) * ( p1000mb - grid%p_top )
 538   ENDDO
 539
 540 ! get fire mesh dimensions
 541     CALL get_ijk_from_subgrid (  grid ,                   &
 542        ifds,ifde, jfds,jfde,kfds,kfde, &
 543        ifms,ifme, jfms,jfme,kfms,kfme,  &
 544        ifts,ifte, jfts,jfte,kfts,kfte)
 545
 546   write (6,*)' ******** SFIRE ideal initialization ********'
 547
 548   ! fire grid step size
 549   fdx = grid%dx/grid%sr_x
 550   fdy = grid%dy/grid%sr_y
 551   ! refinement ratios
 552   ir = grid%sr_x
 553   jr = grid%sr_y
 554
 555   write (6,*)' atm  mesh step ',grid%dx,grid%dy
 556   write (6,*)' fire mesh step ',fdx,fdy
 557   write (6,*)' refinement ratio ',grid%sr_x,grid%sr_y
 558   write (6,*)' atm  domain bounds ',ids,ide, jds,jde,kds,kde
 559   write (6,*)' atm  memory bounds ',ims,ime, jms,jme,kms,kme
 560   write (6,*)' atm  tile   bounds ',its,ite, jts,jte,kts,kte
 561   write (6,*)' fire domain bounds ',ifds,ifde, jfds,jfde,kfds,kfde
 562   write (6,*)' fire memory bounds ',ifms,ifme, jfms,jfme,kfms,kfme
 563   write (6,*)' fire tile   bounds ',ifts,ifte, jfts,jfte,kfts,kfte
 564   write (6,*)' Note that atm mesh and fire mesh are cell-centered'
 565
 566 ! set ideal coordinates
 567   call set_ideal_coord( fdx,fdy, &
 568                 ifds,ifde,jfds,jfde,  &
 569                 ifms,ifme,jfms,jfme,  &
 570                 ifts,ifte,jfts,jfte,  &
 571                 grid%fxlong,grid%fxlat          )
 572 ! Avoid setting atmospheric coordinates when using nested domains ME
 573   if (config_flags%max_dom.eq.1) then
 574       call set_ideal_coord( grid%dx,grid%dy, &
 575                   ids,ide,jds,jde,  &
 576                   ims,ime,jms,jme,  &
 577                   its,ite,jts,jte,  &
 578                   grid%xlong,grid%xlat          )
 579   endif
 580
 581 ! set terrain height
 582
 583   DO j=jts,jte
 584   DO i=its,ite
 585     grid%ht(i,j) = 0.
 586   ENDDO
 587   ENDDO
 588
 589   if(config_flags%fire_fmc_read.eq.2) then
 590           write(6,*)'Reading fuel moisture from file input_fmc_g'
 591           call read_array_2d_real   ('input_fmc_g',grid%fmc_g,    ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
 592   endif
 593
 594   call read_namelist_fire(.false.)  ! read fuel coefficienrs
 595
 596   if(config_flags%fire_fmc_read.eq.0) then
 597           write(6,*)'Setting fuel moisture in wrfinput to constant ', fuelmc_g
 598           do j = jfds,jfde
 599              do i= ifds,ifde
 600                 grid%fmc_g(i,j)=fuelmc_g
 601              enddo
 602           enddo
 603   endif
 604
 605   call print_2d_stats(ifds,ifde,jfds,jfde, &
 606                          ifms,ifme,jfms,jfme, &
 607                          grid%fmc_g, 'fmc_g')
 608
 609
 610   if(config_flags%fire_fuel_read.eq.2) then
 611       write(6,*)'Reading fuel map from file input_fc'
 612       call read_array_2d_real('input_fc',grid%nfuel_cat,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
 613   endif
 614
 615   have_fire_grad=.false.
 616   have_atm_grad=.false.
 617   have_fire_ht=.false.
 618
 619   !******* set terrain height
 620
 621   ! copy params from the namelist
 622   mtn_type = config_flags%fire_mountain_type
 623   mtn_xs   = config_flags%fire_mountain_start_x
 624   mtn_ys   = config_flags%fire_mountain_start_y
 625   mtn_xe   = config_flags%fire_mountain_end_x
 626   mtn_ye   = config_flags%fire_mountain_end_y
 627   mtn_ht   = config_flags%fire_mountain_height
 628
 629   IF(mtn_type .ne. 0)THEN
 630
 631     ! idealized mountain
 632
 633     ! atmospheric grid coordinates of the mountain
 634     mtn_axs = mtn_xs/grid%dx + ids - 0.5
 635     mtn_axe = mtn_xe/grid%dx + ids - 0.5
 636     mtn_ays = mtn_ys/grid%dy + jds - 0.5
 637     mtn_aye = mtn_ye/grid%dy + jds - 0.5
 638
 639     ! fire grid coordinates of the mountain
 640     mtn_fxs = mtn_xs/fdx + ifds - 0.5
 641     mtn_fxe = mtn_xe/fdx + ifds - 0.5
 642     mtn_fys = mtn_ys/fdy + jfds - 0.5
 643     mtn_fye = mtn_ye/fdy + jfds - 0.5
 644
 645     write(6,*)' Mountain height ',mtn_ht,' type',mtn_type
 646     write(6,*)' Mountain (m) LL=(0,0) ',mtn_xs,':',mtn_xe,' ',mtn_ys,':',mtn_ye
 647     write(6,*)' Mountain on atm grid  ',mtn_axs,':',mtn_axe,' ',mtn_ays,':',mtn_aye
 648     write(6,*)' Mountain on fire grid ',mtn_fxs,':',mtn_fxe,' ',mtn_fys,':',mtn_fye
 649
 650     mtn_max = 0.
 651     DO j=jts,jte
 652     DO i=its,ite
 653       mtn_x = pi + 2*pi* max(0. , min( (i - mtn_axs)/(mtn_axe - mtn_axs), 1. ))
 654       mtn_y = pi + 2*pi* max(0. , min( (j - mtn_ays)/(mtn_aye - mtn_ays), 1. ))
 655       SELECT CASE(mtn_type)
 656         CASE (1) ! circ/elliptic mountain
 657           mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y))
 658         CASE (2)  ! EW ridge
 659           mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y))
 660         CASE (3)  ! NS ridge
 661           mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x))
 662         CASE DEFAULT
 663           call wrf_error_fatal ( ' bad fire_mountain_type ' )
 664       END SELECT
 665       mtn_max = max(mtn_max, mtn_z)
 666       grid%ht(i,j) = mtn_z
 667     ENDDO
 668     ENDDO
 669
 670     write(6, *)' Atm  tile ',its,':',ite,' ',jts,':',jte,' max terrain height ',mtn_max
 671
 672     DO j=jfts,jfte
 673         DO i=ifts,ifte
 674           mtn_x = pi + 2*pi* max(0. , min( (i - mtn_fxs)/(mtn_fxe - mtn_fxs), 1. ))
 675           mtn_y = pi + 2*pi* max(0. , min( (j - mtn_fys)/(mtn_fye - mtn_fys), 1. ))
 676           SELECT CASE(mtn_type)
 677             CASE (1) ! circ/elliptic mountain
 678               mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y))
 679             CASE (2)  ! EW ridge
 680               mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y))
 681             CASE (3)  ! NS ridge
 682               mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x))
 683             CASE DEFAULT
 684               call wrf_error_fatal ( ' bad fire_mountain_type ' )
 685           END SELECT
 686           grid%zsf(i,j) = mtn_z
 687         ENDDO
 688     ENDDO
 689
 690     have_fire_ht=.true.
 691
 692   ELSE ! mtn_type
 693
 694     if(config_flags%fire_read_atm_ht)then !
 695       call read_array_2d_real('input_ht',grid%ht,ids,ide,jds,jde,ims,ime,jms,jme)
 696       ! no flag - we always have the terrain height on atm mesh, zero if not set
 697     endif
 698
 699     if(config_flags%fire_read_fire_ht)then !
 700       call read_array_2d_real('input_zsf',grid%zsf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
 701       have_fire_ht=.true.
 702     endif
 703
 704     if(config_flags%fire_read_atm_grad)then !
 705       call crash('Reading terrain gradient on atm mesh from file not supported.')
 706       have_atm_grad=.true.
 707     endif
 708
 709     if(config_flags%fire_read_fire_grad)then !
 710       call read_array_2d_real('input_dzdxf',grid%dzdxf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
 711       call read_array_2d_real('input_dzdyf',grid%dzdyf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
 712       have_fire_grad=.true.
 713     endif
 714   ENDIF  ! mtn_type
 715
 716   if(have_fire_ht)then
 717       write(6, *)'Fine-resolution terrain height on the fire mesh used.'
 718   else
 719       write(6,*)'Interpolating the terrain height from the atm mesh to the fire mesh'
 720       call interpolate_2d(  &
 721         ims,ime,jms,jme, & ! memory dims atm grid tile
 722         its,ite,jts,jte, & ! where atm grid values set
 723         ifms,ifme,jfms,jfme,    & ! array dims fire grid
 724         ifts,ifte,jfts,jfte,  & ! dimensions fire grid tile
 725         ir,jr,                & ! refinement ratio
 726         real(ids),real(jds),ifds+(ir-1)*0.5,jfds+(jr-1)*0.5, & ! line up by lower left corner of domain
 727         grid%ht,                                      & ! atm grid arrays in
 728         grid%zsf)                                      ! fire grid arrays out
 729       have_fire_ht=.true.
 730   endif
 731
 732
 733   if(have_fire_grad)then
 734      write(6, *)'Fine-resolution terrain gradient on the fire mesh used.'
 735   else
 736
 737      write(6,*)'Computing the terrain gradient from fire mesh height'
 738      if(.not.have_fire_ht)then
 739         write(6,*)'WARNING: Fire mesh terrain height not given, setting to zero'
 740         do j=jfts,jfte
 741            do i=ifts,ifte
 742               grid%zsf(i,j) = 0.
 743            enddo
 744         enddo
 745       endif
 746
 747       ! extend the terrain height one beyond the domain
 748       call continue_at_boundary(1,1,0., & ! do x direction or y direction
 749               ifms,ifme,jfms,jfme, &            ! memory dims
 750               ifds,ifde,jfds,jfde, &            ! domain dims
 751               ifds,ifde,jfds,jfde, &            ! patch dims = domain, not parallel!
 752               ifts,ifte,jfts,jfte, &            ! tile dims
 753               iots,iote,jots,jote, &            ! tile dims out
 754               grid%zsf)                         ! array
 755
 756       ! compute the terrain gradient by differencing
 757       do j=jfts,jfte
 758          do i=ifts,ifte
 759             grid%dzdxf(i,j) = (grid%zsf(i+1,j)-grid%zsf(i-1,j))/(2.*fdx)
 760             grid%dzdyf(i,j) = (grid%zsf(i,j+1)-grid%zsf(i,j-1))/(2.*fdy)
 761          enddo
 762       enddo
 763       have_fire_grad=.true.
 764    endif ! have_fire_grad
 765
 766    if(.not.have_fire_grad)call crash('Fire mesh terrain gradient not set')
 767
 768     mtn_max = 0.
 769     DO j=jts,jte
 770       DO i=its,ite
 771         mtn_max = max(mtn_max, grid%ht(i,j))
 772       ENDDO
 773     ENDDO
 774     write(6, *)' Max terrain height on the atmosphere mesh ',mtn_max
 775
 776     mtn_max = 0.
 777     grad_max =0.
 778     DO j=jfts,jfte
 779       DO i=ifts,ifte
 780         mtn_max = max(mtn_max, grid%zsf(i,j))
 781         grad_max = max( grad_max, sqrt(grid%dzdxf(i,j)**2+grid%dzdyf(i,j)**2) )
 782       ENDDO
 783     ENDDO
 784     write(6, *)' Max terrain height on the fire mesh       ',mtn_max
 785     write(6, *)' Max terrain gradient on the fire mesh     ',grad_max
 786
 787 ! the rest of initialization dependent on the atmosphere grid terrain height set
 788
 789   DO j=jts,jte
 790   DO i=its,ite
 791     grid%phb(i,1,j) = g * grid%ht(i,j)
 792     grid%ph0(i,1,j) = g * grid%ht(i,j)
 793   ENDDO
 794   ENDDO
 795
 796   DO J = jts, jte
 797   DO I = its, ite
 798
 799     p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
 800     grid%mub(i,j) = p_surf-grid%p_top
 801
 802 !  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
 803 !  interp theta (from interp) and compute 1/rho from eqn. of state
 804
 805     DO K = 1, kte-1
 806       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 807       grid%pb(i,k,j) = p_level
 808       grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
 809       grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 810     ENDDO
 811
 812 !  calc hydrostatic balance (alternatively we could interp the geopotential from the
 813 !  sounding, but this assures that the base state is in exact hydrostatic balance with
 814 !  respect to the model eqns.
 815
 816     DO k  = 2,kte
 817       grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*(grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))*grid%alb(i,k-1,j)
 818     ENDDO
 819
 820   ENDDO
 821   ENDDO
 822
 823   IF ( wrf_dm_on_monitor() ) THEN
 824     write(6,*) ' ptop is ',grid%p_top
 825     write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
 826   ENDIF
 827
 828 !  calculate full state for each column - this includes moisture.
 829
 830   write(6,*) ' getting moist sounding for full state '
 831   dry_sounding = .false.
 832   CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 833
 834   DO J = jts, min(jde-1,jte)
 835   DO I = its, min(ide-1,ite)
 836
 837 !  At this point grid%p_top is already set. find the DRY mass in the column
 838 !  by interpolating the DRY pressure.
 839
 840    pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
 841
 842 !  compute the perturbation mass and the full mass
 843
 844     grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
 845     grid%mu_2(i,j) = grid%mu_1(i,j)
 846     grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
 847
 848 ! given the dry pressure and coordinate system, interp the potential
 849 ! temperature and qv
 850
 851     do k=1,kde-1
 852
 853       p_level = grid%c3h(k)*(pd_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 854
 855       moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
 856       grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
 857       grid%t_2(i,k,j)          = grid%t_1(i,k,j)
 858
 859
 860     enddo
 861
 862 !  integrate the hydrostatic equation (from the RHS of the bigstep
 863 !  vertical momentum equation) down from the top to get grid%p.
 864 !  first from the top of the model to the top pressure
 865
 866     k = kte-1  ! top level
 867
 868     qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
 869     qvf2 = 1./(1.+qvf1)
 870     qvf1 = qvf1*qvf2
 871
 872     grid%p(i,k,j) = - 0.5*((grid%c1f(k+1)*grid%mu_1(i,j))+qvf1*(grid%c1f(k+1)*grid%mub(i,j)+grid%c2f(k+1)))/grid%rdnw(k)/qvf2
 873     qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 874     grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 875                 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 876     grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 877
 878 !  down the column
 879
 880     do k=kte-2,1,-1
 881       qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
 882       qvf2 = 1./(1.+qvf1)
 883       qvf1 = qvf1*qvf2
 884       grid%p(i,k,j) = grid%p(i,k+1,j) - ((grid%c1f(k+1)*grid%mu_1(i,j)) + qvf1*(grid%c1f(k+1)*grid%mub(i,j)+grid%c2f(k+1)))/qvf2/grid%rdn(k+1)
 885       qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 886       grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 887                   (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 888       grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 889     enddo
 890
 891 !  this is the hydrostatic equation used in the model after the
 892 !  small timesteps.  In the model, grid%al (inverse density)
 893 !  is computed from the geopotential.
 894
 895
 896     grid%ph_1(i,1,j) = 0.
 897     DO k  = 2,kte
 898       grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
 899                    ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ &
 900                     (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
 901
 902       grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
 903       grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
 904     ENDDO
 905
 906     IF ( wrf_dm_on_monitor() ) THEN
 907     if((i==2) .and. (j==2)) then
 908      write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
 909                               grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
 910                               grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
 911     endif
 912     ENDIF
 913
 914   ENDDO
 915   ENDDO
 916
 917 ! checking if the perturbation (bubble) is to be applied
 918   IF ((delt/=0.) .and. (x_rad > 0.) &
 919                  .and. (y_rad > 0.) &
 920                  .and. (z_rad > 0.)) THEN
 921 !  thermal perturbation to kick off convection
 922
 923   write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
 924   write(6,'(A23,f18.16)') ' delt for perturbation ',delt
 925   write(6,'(A30,f18.12)') ' x radius of the perturbation ' ,x_rad
 926   write(6,'(A30,f18.12)') ' y radius of the perturbation ' ,y_rad
 927   write(6,'(A30,f18.12)') ' z radius of the perturbation ' ,z_rad
 928   write(6,'(A30,f18.12)') ' height of the perturbation   ' ,hght_pert
 929
 930   DO J = jts, min(jde-1,jte)
 931     yrad = config_flags%dy*float(j-nyc)/y_rad
 932 !   yrad = 0.
 933     DO I = its, min(ide-1,ite)
 934       xrad = config_flags%dx*float(i-nxc)/x_rad
 935 !     xrad = 0.
 936       DO K = 1, kte-1
 937
 938 !  put in preturbation theta (bubble) and recalc density.  note,
 939 !  the mass in the column is not changing, so when theta changes,
 940 !  we recompute density and geopotential
 941
 942         zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j)  &
 943                    +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
 944         zrad = (zrad-hght_pert)/z_rad
 945         RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
 946         IF(RAD <= 1.) THEN
 947            grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
 948            grid%t_2(i,k,j)=grid%t_1(i,k,j)
 949            qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 950            grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 951                         (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 952            grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 953         ENDIF
 954       ENDDO
 955
 956 !  rebalance hydrostatically
 957
 958       DO k  = 2,kte
 959         grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
 960                      ((grid%c1h(k-1)*grid%mub(i,j)+grid%c2h(k-1))+(grid%c1h(k-1)*grid%mu_1(i,j)))*grid%al(i,k-1,j)+ &
 961                       (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
 962
 963         grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
 964         grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
 965       ENDDO
 966
 967     ENDDO
 968   ENDDO
 969
 970 !End of setting up the perturbation (bubble)
 971 ENDIF
 972
 973    IF ( wrf_dm_on_monitor() ) THEN
 974    write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
 975    write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
 976    do k=1,kde-1
 977      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
 978                                       grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
 979                                       grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
 980    enddo
 981
 982    write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
 983    do k=1,kde-1
 984      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
 985                                       grid%p(1,k,1), grid%al(1,k,1), &
 986                                       grid%t_1(1,k,1), moist(1,k,1,P_QV)
 987    enddo
 988    ENDIF
 989
 990 ! interp v
 991
 992   DO J = jts, jte
 993   DO I = its, min(ide-1,ite)
 994
 995     IF (j == jds) THEN
 996       z_at_v = grid%phb(i,1,j)/g
 997     ELSE IF (j == jde) THEN
 998       z_at_v = grid%phb(i,1,j-1)/g
 999     ELSE
1000       z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
1001     END IF
1002     p_surf = interp_0( p_in, zk, z_at_v, nl_in )
1003
1004     DO K = 1, kte-1
1005       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
1006       grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
1007       grid%v_2(i,k,j) = grid%v_1(i,k,j)
1008     ENDDO
1009
1010   ENDDO
1011   ENDDO
1012
1013 ! interp u
1014
1015   DO J = jts, min(jde-1,jte)
1016   DO I = its, ite
1017
1018     IF (i == ids) THEN
1019       z_at_u = grid%phb(i,1,j)/g
1020     ELSE IF (i == ide) THEN
1021       z_at_u = grid%phb(i-1,1,j)/g
1022     ELSE
1023       z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
1024     END IF
1025
1026     p_surf = interp_0( p_in, zk, z_at_u, nl_in )
1027
1028     DO K = 1, kte-1
1029       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
1030       grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
1031       grid%u_2(i,k,j) = grid%u_1(i,k,j)
1032     ENDDO
1033
1034   ENDDO
1035   ENDDO
1036
1037 !  set w
1038
1039   DO J = jts, min(jde-1,jte)
1040   DO K = kts, kte
1041   DO I = its, min(ide-1,ite)
1042     grid%w_1(i,k,j) = 0.
1043     grid%w_2(i,k,j) = 0.
1044   ENDDO
1045   ENDDO
1046   ENDDO
1047
1048 !  set a few more things
1049
1050   DO J = jts, min(jde-1,jte)
1051   DO K = kts, kte-1
1052   DO I = its, min(ide-1,ite)
1053     grid%h_diabatic(i,k,j) = 0.
1054   ENDDO
1055   ENDDO
1056   ENDDO
1057
1058   IF ( wrf_dm_on_monitor() ) THEN
1059   DO k=1,kte-1
1060     grid%t_base(k) = grid%t_1(1,k,1)
1061     grid%qv_base(k) = moist(1,k,1,P_QV)
1062     grid%u_base(k) = grid%u_1(1,k,1)
1063     grid%v_base(k) = grid%v_1(1,k,1)
1064     grid%z_base(k) = 0.5*(grid%phb(1,k,1)+grid%phb(1,k+1,1)+grid%ph_1(1,k,1)+grid%ph_1(1,k+1,1))/g
1065   ENDDO
1066   ENDIF
1067   CALL wrf_dm_bcast_real( grid%t_base , kte )
1068   CALL wrf_dm_bcast_real( grid%qv_base , kte )
1069   CALL wrf_dm_bcast_real( grid%u_base , kte )
1070   CALL wrf_dm_bcast_real( grid%v_base , kte )
1071   CALL wrf_dm_bcast_real( grid%z_base , kte )
1072
1073   DO J = jts, min(jde-1,jte)
1074   DO I = its, min(ide-1,ite)
1075      thtmp   = grid%t_2(i,1,j)+t0
1076      ptmp    = grid%p(i,1,j)+grid%pb(i,1,j)
1077      temp(1) = thtmp * (ptmp/p1000mb)**rcp
1078      thtmp   = grid%t_2(i,2,j)+t0
1079      ptmp    = grid%p(i,2,j)+grid%pb(i,2,j)
1080      temp(2) = thtmp * (ptmp/p1000mb)**rcp
1081      thtmp   = grid%t_2(i,3,j)+t0
1082      ptmp    = grid%p(i,3,j)+grid%pb(i,3,j)
1083      temp(3) = thtmp * (ptmp/p1000mb)**rcp
1084
1085 !     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3) !AK/AK it is already declared via namelist.input if sfc_init=.true.
1086 !     grid%tmn(I,J)=grid%tsk(I,J)-0.5                                  !AK/AK it is already declared via namelist.input if sfc_init=.true.
1087   ENDDO
1088   ENDDO
1089
1090   IF (.NOT.sfc_init) THEN
1091   write(6,*) ' setting tsk and tmn default'
1092     DO J = jts, min(jde-1,jte)
1093     DO I = its, min(ide-1,ite)
1094      grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
1095      grid%tmn(I,J)=grid%tsk(I,J)-0.5
1096    ENDDO
1097    ENDDO
1098  ENDIF
1099
1100       !  Save the dry perturbation potential temperature.
1101
1102       DO j = jts, min(jde-1,jte)
1103          DO k = kts, kte
1104             DO i = its, min(ide-1,ite)
1105                grid%th_phy_m_t0(i,k,j) = grid%t_2(i,k,j)
1106             END DO
1107          END DO
1108       END DO
1109
1110  !  Turn dry potential temperature into moist potential temperature
1111       !  at the very end of this routine
1112       !  This field will be in the model IC and and used to construct the
1113       !  BC file.
1114
1115       IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
1116       DO J  = jts, min(jde-1,jte)
1117          DO K = kts, kte-1
1118             DO I = its, min(ide-1,ite)
1119                grid%t_2(i,k,j) = ( grid%t_2(i,k,j) + T0 ) * (1. + (R_v/R_d) * moist(i,k,j,p_qv)) - T0
1120             END DO
1121          END DO
1122       END DO
1123       ENDIF
1124
1125  END SUBROUTINE init_domain_rk
1126
1127    SUBROUTINE init_module_initialize
1128    END SUBROUTINE init_module_initialize
1129
1130 !---------------------------------------------------------------------
1131
1132 !  test driver for get_sounding
1133 !
1134 !      implicit none
1135 !      integer n
1136 !      parameter(n = 1000)
1137 !      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
1138 !      logical dry
1139 !      integer nl,k
1140 !
1141 !      dry = .false.
1142 !      dry = .true.
1143 !      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
1144 !      write(6,*) ' input levels ',nl
1145 !      write(6,*) ' sounding '
1146 !      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
1147 !      do k=1,nl
1148 !        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k)
1149 !      enddo
1150 !      end
1151 !
1152 !---------------------------------------------------------------------------
1153
1154       subroutine get_sounding( zk, p, p_dry, theta, rho, &
1155                                u, v, qv, dry, nl_max, nl_in )
1156       implicit none
1157
1158       integer nl_max, nl_in
1159       real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
1160            u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
1161       logical dry
1162
1163       integer n
1164       parameter(n=1000)
1165       logical debug
1166       parameter( debug = .true.)
1167
1168 ! input sounding data
1169
1170       real p_surf, th_surf, qv_surf
1171       real pi_surf, pi(n)
1172       real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
1173
1174 ! diagnostics
1175
1176       real rho_surf, p_input(n), rho_input(n)
1177       real pm_input(n)  !  this are for full moist sounding
1178
1179 ! local data
1180
1181       real r
1182       parameter (r = r_d)
1183       integer k, it, nl
1184       real qvf, qvf1, dz
1185
1186 !  first, read the sounding
1187
1188       call read_sounding( p_surf, th_surf, qv_surf, &
1189                           h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
1190
1191       if(dry) then
1192        do k=1,nl
1193          qv_input(k) = 0.
1194        enddo
1195       endif
1196
1197       if(debug) write(6,*) ' number of input levels = ',nl
1198
1199         nl_in = nl
1200         if(nl_in .gt. nl_max ) then
1201           write(6,*) ' too many levels for input arrays ',nl_in,nl_max
1202           call wrf_error_fatal ( ' too many levels for input arrays ' )
1203         end if
1204
1205 !  compute diagnostics,
1206 !  first, convert qv(g/kg) to qv(g/g)
1207
1208       do k=1,nl
1209         qv_input(k) = 0.001*qv_input(k)
1210       enddo
1211
1212       p_surf = 100.*p_surf  ! convert to pascals
1213       qvf = 1. + rvovrd*qv_input(1)
1214       rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
1215       pi_surf = (p_surf/p1000mb)**(r/cp)
1216
1217       if(debug) then
1218         write(6,*) ' surface density is ',rho_surf
1219         write(6,*) ' surface pi is      ',pi_surf
1220       end if
1221
1222
1223 !  integrate moist sounding hydrostatically, starting from the
1224 !  specified surface pressure
1225 !  -> first, integrate from surface to lowest level
1226
1227           qvf = 1. + rvovrd*qv_input(1)
1228           qvf1 = 1. + qv_input(1)
1229           rho_input(1) = rho_surf
1230           dz = h_input(1)
1231           do it=1,10
1232             pm_input(1) = p_surf &
1233                     - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
1234             rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
1235           enddo
1236
1237 ! integrate up the column
1238
1239           do k=2,nl
1240             rho_input(k) = rho_input(k-1)
1241             dz = h_input(k)-h_input(k-1)
1242             qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
1243             qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
1244
1245             do it=1,10
1246               pm_input(k) = pm_input(k-1) &
1247                       - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
1248               rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
1249             enddo
1250           enddo
1251
1252 !  we have the moist sounding
1253
1254 !  next, compute the dry sounding using p at the highest level from the
1255 !  moist sounding and integrating down.
1256
1257         p_input(nl) = pm_input(nl)
1258
1259           do k=nl-1,1,-1
1260             dz = h_input(k+1)-h_input(k)
1261             p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
1262           enddo
1263
1264
1265         do k=1,nl
1266
1267           zk(k) = h_input(k)
1268           p(k) = pm_input(k)
1269           p_dry(k) = p_input(k)
1270           theta(k) = th_input(k)
1271           rho(k) = rho_input(k)
1272           u(k) = u_input(k)
1273           v(k) = v_input(k)
1274           qv(k) = qv_input(k)
1275
1276         enddo
1277
1278      if(debug) then
1279       write(6,*) ' sounding '
1280       write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
1281       do k=1,nl
1282         write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k)
1283       enddo
1284
1285      end if
1286
1287       end subroutine get_sounding
1288
1289 !-------------------------------------------------------
1290
1291       subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
1292       implicit none
1293       integer n,nl
1294       real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
1295       logical end_of_file
1296       logical debug
1297
1298       integer k
1299
1300       open(unit=10,file='input_sounding',form='formatted',status='old')
1301       rewind(10)
1302       read(10,*) ps, ts, qvs
1303       if(debug) then
1304         write(6,*) ' input sounding surface parameters '
1305         write(6,*) ' surface pressure (mb) ',ps
1306         write(6,*) ' surface pot. temp (K) ',ts
1307         write(6,*) ' surface mixing ratio (g/kg) ',qvs
1308       end if
1309
1310       end_of_file = .false.
1311       k = 0
1312
1313       do while (.not. end_of_file)
1314
1315         read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
1316         k = k+1
1317         if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
1318         go to 110
1319  100    end_of_file = .true.
1320  110    continue
1321       enddo
1322
1323       nl = k
1324
1325       close(unit=10,status = 'keep')
1326
1327       end subroutine read_sounding
1328
1329 END MODULE module_initialize_ideal