dyn_em/module_initialize_scm_xy.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
  27 #ifdef DM_PARALLEL
  28    USE module_dm
  29 #endif
  30
  31
  32 CONTAINS
  33
  34
  35 !-------------------------------------------------------------------
  36 ! this is a wrapper for the solver-specific init_domain routines.
  37 ! Also dereferences the grid variables and passes them down as arguments.
  38 ! This is crucial, since the lower level routines may do message passing
  39 ! and this will get fouled up on machines that insist on passing down
  40 ! copies of assumed-shape arrays (by passing down as arguments, the
  41 ! data are treated as assumed-size -- ie. f77 -- arrays and the copying
  42 ! business is avoided).  Fie on the F90 designers.  Fie and a pox.
  43 ! NOTE:  Modified to remove all but arrays of rank 4 or more from the
  44 !        argument list.  Arrays with rank>3 are still problematic due to the
  45 !        above-noted fie- and pox-ities.  TBH 20061129.
  46
  47    SUBROUTINE init_domain ( grid )
  48
  49    IMPLICIT NONE
  50
  51    !  Input data.
  52    TYPE (domain), POINTER :: grid
  53    !  Local data.
  54    INTEGER :: idum1, idum2
  55
  56    CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
  57
  58      CALL init_domain_rk( grid &
  59 !
  60 #include "actual_new_args.inc"
  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
  73    USE module_optional_input
  74    IMPLICIT NONE
  75
  76    !  Input data.
  77    TYPE (domain), POINTER :: grid
  78
  79 # include "dummy_new_decl.inc"
  80
  81    TYPE (grid_config_rec_type)              :: config_flags
  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 ! JPH should add a read to a config file with:
  91 ! ----- check to make sure everything is initialized from the LU index, etc.
  92 ! ----- need to make a dummy category?
  93 ! cen_lat, cen_lon
  94 ! land-use category
  95 ! soil category
  96
  97    ! Local data
  98    INTEGER, PARAMETER :: nl_max = 1000
  99    REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
 100    INTEGER :: nl_in
 101
 102    INTEGER :: ii, im1, jj, jm1, loop, error, fid, lm
 103    REAL    :: B1, B2, B3, B4, B5
 104    REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
 105    REAL    :: xrad, yrad, zrad, rad, cof1, cof2
 106    REAL    :: pi, rnd
 107
 108 !  stuff from original initialization that has been dropped from the Registry
 109    REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd
 110    REAL    :: qvf1, qvf2, pd_surf
 111    INTEGER :: it
 112    real :: thtmp, ptmp, temp(3)
 113    real :: zsfc
 114
 115    LOGICAL :: moisture_init
 116    LOGICAL :: dry_sounding
 117    character (len=256) :: mminlu2
 118
 119 ! soil input
 120    INTEGER :: ns_input
 121    REAL    :: tmn_input, tsk_input
 122    REAL    :: zs_input(100),tslb_input(100),smois_input(100)
 123    LOGICAL :: real_soil = .true.
 124
 125    REAL    :: zrwa(200), zwa(200)
 126
 127
 128 #ifdef DM_PARALLEL
 129 #    include "data_calls.inc"
 130 #endif
 131
 132
 133    SELECT CASE ( model_data_order )
 134          CASE ( DATA_ORDER_ZXY )
 135    kds = grid%sd31 ; kde = grid%ed31 ;
 136    ids = grid%sd32 ; ide = grid%ed32 ;
 137    jds = grid%sd33 ; jde = grid%ed33 ;
 138
 139    kms = grid%sm31 ; kme = grid%em31 ;
 140    ims = grid%sm32 ; ime = grid%em32 ;
 141    jms = grid%sm33 ; jme = grid%em33 ;
 142
 143    kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
 144    its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
 145    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 146          CASE ( DATA_ORDER_XYZ )
 147    ids = grid%sd31 ; ide = grid%ed31 ;
 148    jds = grid%sd32 ; jde = grid%ed32 ;
 149    kds = grid%sd33 ; kde = grid%ed33 ;
 150
 151    ims = grid%sm31 ; ime = grid%em31 ;
 152    jms = grid%sm32 ; jme = grid%em32 ;
 153    kms = grid%sm33 ; kme = grid%em33 ;
 154
 155    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 156    jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
 157    kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
 158          CASE ( DATA_ORDER_XZY )
 159    ids = grid%sd31 ; ide = grid%ed31 ;
 160    kds = grid%sd32 ; kde = grid%ed32 ;
 161    jds = grid%sd33 ; jde = grid%ed33 ;
 162
 163    ims = grid%sm31 ; ime = grid%em31 ;
 164    kms = grid%sm32 ; kme = grid%em32 ;
 165    jms = grid%sm33 ; jme = grid%em33 ;
 166
 167    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 168    kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
 169    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 170
 171    END SELECT
 172
 173
 174    pi = 2.*asin(1.0)
 175    write(6,*) ' pi is ',pi
 176
 177    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 178
 179 ! here we check to see if the boundary conditions are set properly
 180
 181    CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
 182
 183    moisture_init = .true.
 184
 185     grid%itimestep=0
 186
 187     mminlu2 = ' '
 188     mminlu2(1:4) = 'USGS'
 189     CALL nl_set_mminlu(1, mminlu2)
 190 !   CALL nl_set_mminlu(1, 'USGS')
 191     CALL nl_set_iswater(1,16)
 192     CALL nl_set_isice(1,3)
 193     CALL nl_set_truelat1(1,0.)
 194     CALL nl_set_truelat2(1,0.)
 195     CALL nl_set_moad_cen_lat (1,0.)
 196     CALL nl_set_stand_lon(1,0.)
 197     CALL nl_set_pole_lon (1,0.)
 198     CALL nl_set_pole_lat (1,90.)
 199     CALL nl_set_map_proj(1,0)
 200 !   CALL model_to_grid_config_rec(1,model_config_rec,config_flags)
 201     CALL nl_get_iswater(1,grid%iswater)
 202
 203 !  here we initialize data that currently is not initialized
 204 !  in the input data
 205
 206     DO j = jts, jte
 207       DO i = its, ite
 208          grid%msft(i,j)     = 1.
 209          grid%msfu(i,j)     = 1.
 210          grid%msfv(i,j)     = 1.
 211          grid%msftx(i,j)    = 1.
 212          grid%msfty(i,j)    = 1.
 213          grid%msfux(i,j)    = 1.
 214          grid%msfuy(i,j)    = 1.
 215          grid%msfvx(i,j)    = 1.
 216          grid%msfvx_inv(i,j)    = 1.
 217          grid%msfvy(i,j)    = 1.
 218          grid%sina(i,j)     = 0.
 219          grid%cosa(i,j)     = 1.
 220          grid%e(i,j)        = 2.0*EOMEG*cos(config_flags%scm_lat*DEGRAD)
 221          grid%f(i,j)        = 2.0*EOMEG*sin(config_flags%scm_lat*DEGRAD)
 222          grid%xlat(i,j)     = config_flags%scm_lat
 223          grid%xlong(i,j)    = config_flags%scm_lon
 224          grid%xland(i,j)     = 1.
 225          grid%landmask(i,j)  = 1.
 226          grid%lu_index(i,j)  = config_flags%scm_lu_index
 227       END DO
 228    END DO
 229
 230 ! for LSM, additional variables need to be initialized
 231
 232 !  other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
 233
 234 !     CASE (SLABSCHEME)
 235
 236 !     CASE (LSMSCHEME)
 237
 238 ! JPH free of snow and ice, and only valid over land
 239         DO j = jts , MIN(jde-1,jte)
 240            DO i = its , MIN(ide-1,ite)
 241               grid%vegfra(i,j) = config_flags%scm_vegfra
 242               grid%canwat(i,j) = config_flags%scm_canwat
 243               grid%isltyp(i,j)  = config_flags%scm_isltyp
 244               grid%ivgtyp(i,j)  = config_flags%scm_lu_index
 245               grid%xice(i,j) = 0.
 246               grid%snow(i,j) = 0.
 247            END DO
 248         END DO
 249
 250 !     CASE (RUCLSMSCHEME)
 251
 252 !  END SELECT other_masked_fields
 253
 254    grid%step_number = 0
 255
 256    IF ( real_soil ) THEN ! from input file
 257
 258       IF (config_flags%sf_surface_physics .NE. 2) WRITE (6, *)   &
 259          'If using LSM option other than Noah, must edit input_soil file in test/em_scm_xy/ directory'
 260
 261       CALL read_soil(100,ns_input,tmn_input,tsk_input,zs_input,tslb_input,smois_input)
 262
 263       CALL init_module_optional_input(grid,config_flags)
 264       num_st_levels_input = ns_input
 265       num_sm_levels_input = ns_input
 266       num_sw_levels_input = ns_input
 267       DO k = 1,ns_input
 268          st_levels_input(k) = zs_input(k)*100.0 ! to cm
 269          sm_levels_input(k) = zs_input(k)*100.0 ! to cm
 270          sw_levels_input(k) = zs_input(k)*100.0 ! to cm
 271          st_input(:,k+1,:) = tslb_input(k)
 272          sm_input(:,k+1,:) = smois_input(k)
 273          sw_input(:,k+1,:) = smois_input(k)
 274       ENDDO
 275
 276       grid%tsk = tsk_input
 277       grid%sst = tsk_input
 278       grid%tmn = tmn_input
 279
 280       flag_soil_layers  = 0 ! go ahead and put skin temp in
 281       flag_soil_levels  = 0 ! go ahead and put skin moisture in
 282       flag_sst          = 0 ! don't modify for ocean
 283       flag_tavgsfc      = 0
 284       flag_soilhgt      = 0
 285
 286       CALL process_soil_real ( grid%tsk , grid%tmn , grid%tavgsfc, &
 287                    grid%landmask , grid%sst , grid%ht, grid%toposoil, &
 288                    st_input , sm_input , sw_input , &
 289                    st_levels_input , sm_levels_input , sw_levels_input , &
 290                    grid%zs , grid%dzs , grid%flag_sm_adj , grid%tslb , grid%smois , grid%sh2o , &
 291                    flag_sst , flag_tavgsfc, flag_soilhgt, flag_soil_layers, flag_soil_levels,  &
 292                    ids , ide , jds , jde , kds , kde , &
 293                    ims , ime , jms , jme , kms , kme , &
 294                    its , ite , jts , jte , kts , kte , &
 295                    model_config_rec%sf_surface_physics(grid%id) , &
 296                    model_config_rec%num_soil_layers , &
 297                    model_config_rec%real_data_init_type , &
 298                    num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 299                    num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
 300
 301    ELSE ! ideal soil
 302 ! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F
 303       CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, &
 304                      grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, &
 305                      grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, &
 306                      model_config_rec%sf_surface_physics(grid%id), &
 307                                    ids,ide, jds,jde, kds,kde,&
 308                                    ims,ime, jms,jme, kms,kme,&
 309                                    its,ite, jts,jte, kts,kte )
 310
 311     ENDIF
 312
 313     DO j = jts, jte
 314      DO k = kts, kte
 315        DO i = its, ite
 316           grid%ww(i,k,j)     = 0.
 317        END DO
 318      END DO
 319     END DO
 320
 321 ! this is adopted from Wayne Angevine's GABLS case
 322     grid%znw(1) = 1.0
 323     zrwa(kde) = exp((kde-1)/40.)
 324     zwa(kde) = grid%ztop
 325     DO k=2, kde-1
 326        zrwa(k) = exp((k-1)/40.)
 327        zwa(k) = (zrwa(k)-1.) * grid%ztop/(zrwa(kde)-1.)
 328        grid%znw(k) = 1. - (zwa(k) / zwa(kde))
 329     ENDDO
 330     grid%znw(kde) = 0.
 331
 332    DO k=1, kde-1
 333     grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
 334     grid%rdnw(k) = 1./grid%dnw(k)
 335     grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
 336    ENDDO
 337
 338    DO k=2, kde-1
 339     grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
 340     grid%rdn(k) = 1./grid%dn(k)
 341     grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
 342     grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
 343    ENDDO
 344
 345    cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
 346    cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
 347    grid%cf1  = grid%fnp(2) + cof1
 348    grid%cf2  = grid%fnm(2) - cof1 - cof2
 349    grid%cf3  = cof2
 350
 351    grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
 352    grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
 353    grid%rdx = 1./config_flags%dx
 354    grid%rdy = 1./config_flags%dy
 355
 356 !  get the sounding from the ascii sounding file, first get dry sounding and
 357 !  calculate base state
 358
 359   write(6,*) ' getting dry sounding for base state '
 360   dry_sounding = .true.
 361   CALL get_sounding( zsfc, zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 362
 363   write(6,*) ' returned from reading sounding, nl_in is ',nl_in
 364
 365 !  find ptop for the desired ztop (ztop is input from the namelist),
 366 !  and find surface pressure
 367
 368   grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
 369
 370   !  For hybrid coord
 371
 372   DO k=kts, kte
 373      IF      ( config_flags%hybrid_opt .EQ. 0 ) THEN
 374         grid%c3f(k) = grid%znw(k)
 375      ELSE IF ( config_flags%hybrid_opt .EQ. 1 ) THEN
 376         grid%c3f(k) = grid%znw(k)
 377      ELSE IF ( config_flags%hybrid_opt .EQ. 2 ) THEN
 378         B1 = 2. * grid%etac**2 * ( 1. - grid%etac )
 379         B2 = -grid%etac * ( 4. - 3. * grid%etac - grid%etac**3 )
 380         B3 = 2. * ( 1. - grid%etac**3 )
 381         B4 = - ( 1. - grid%etac**2 )
 382         B5 = (1.-grid%etac)**4
 383         grid%c3f(k) = ( B1 + B2*grid%znw(k) + B3*grid%znw(k)**2 + B4*grid%znw(k)**3 ) / B5
 384         IF ( grid%znw(k) .LT. grid%etac ) THEN
 385            grid%c3f(k) = 0.
 386         END IF
 387         IF ( k .EQ. kds ) THEN
 388            grid%c3f(k) = 1.
 389         ELSE IF ( k .EQ. kde ) THEN
 390            grid%c3f(k) = 0.
 391         END IF
 392      ELSE IF ( config_flags%hybrid_opt .EQ. 3 ) THEN
 393         grid%c3f(k) = grid%znw(k)*sin(0.5*3.14159*grid%znw(k))**2
 394         IF      ( k .EQ. kds ) THEN
 395            grid%c3f(k) = 1.
 396         ELSE IF ( k .EQ. kds ) THEN
 397            grid%c3f(kde) = 0.
 398         END IF
 399      ELSE
 400         CALL wrf_message     ( 'ERROR: --- hybrid_opt' )
 401         CALL wrf_message     ( 'ERROR: --- hybrid_opt=0    ==> Standard WRF terrain-following coordinate' )
 402         CALL wrf_message     ( 'ERROR: --- hybrid_opt=1    ==> Standard WRF terrain-following coordinate, hybrid c1, c2, c3, c4' )
 403         CALL wrf_message     ( 'ERROR: --- hybrid_opt=2    ==> Hybrid, Klemp polynomial' )
 404         CALL wrf_message     ( 'ERROR: --- hybrid_opt=3    ==> Hybrid, sin^2' )
 405         CALL wrf_error_fatal ( 'ERROR: --- Invalid option' )
 406      END IF
 407   END DO
 408
 409   !  c4 is a function of c3 and eta.
 410
 411   DO k=1, kde
 412      grid%c4f(k) = ( grid%znw(k) - grid%c3f(k) ) * ( p1000mb - grid%p_top )
 413   ENDDO
 414
 415   !  Now on half levels, just add up and divide by 2 (for c3h).  Use (eta-c3)*(p00-pt) for c4 on half levels.
 416
 417   DO k=1, kde-1
 418      grid%znu(k) = ( grid%znw(k+1) + grid%znw(k) ) * 0.5
 419      grid%c3h(k) = ( grid%c3f(k+1) + grid%c3f(k) ) * 0.5
 420      grid%c4h(k) = ( grid%znu(k) - grid%c3h(k) ) * ( p1000mb - grid%p_top )
 421   ENDDO
 422
 423   !  c1 = d(B)/d(eta).  We define c1f as c1 on FULL levels.  For a vertical difference,
 424   !  we need to use B and eta on half levels.  The k-loop ends up referring to the
 425   !  full levels, neglecting the top and bottom.
 426
 427   DO k=kds+1, kde-1
 428      grid%c1f(k) = ( grid%c3h(k) - grid%c3h(k-1) ) / ( grid%znu(k) - grid%znu(k-1) )
 429   ENDDO
 430
 431   !  The boundary conditions to get the coefficients:
 432   !  1) At k=kts: define d(B)/d(eta) = 1.  This gives us the same value of B and d(B)/d(eta)
 433   !     when doing the sigma-only B=eta.
 434   !  2) At k=kte: define d(B)/d(eta) = 0.  The curve B SMOOTHLY goes to zero, and at the very
 435   !     top, B continues to SMOOTHLY go to zero.  Note that for almost all cases of non B=eta,
 436   !     B is ALREADY=ZERO at the top, so this is a reasonable BC to assume.
 437
 438   grid%c1f(kds) = 1.
 439   IF ( ( config_flags%hybrid_opt .EQ. 0 ) .OR. ( config_flags%hybrid_opt .EQ. 1 ) ) THEN
 440      grid%c1f(kde) = 1.
 441   ELSE
 442      grid%c1f(kde) = 0.
 443   END IF
 444
 445   !  c2 = ( 1. - c1(k) ) * (p00 - pt).  There is no vertical differencing, so we can do the
 446   !  full kds to kde looping.
 447
 448   DO k=kds, kde
 449      grid%c2f(k) = ( 1. - grid%c1f(k) ) * ( p1000mb - grid%p_top )
 450   ENDDO
 451
 452   !  Now on half levels for c1 and c2.  The c1h will result from the full level c3 and full
 453   !  level eta differences.  The c2 value use the half level c1(k).
 454
 455   DO k=1, kde-1
 456      grid%c1h(k) = ( grid%c3f(k+1) - grid%c3f(k) ) / ( grid%znw(k+1) - grid%znw(k) )
 457      grid%c2h(k) = ( 1. - grid%c1h(k) ) * ( p1000mb - grid%p_top )
 458   ENDDO
 459
 460   DO j=jts,jte
 461   DO i=its,ite  ! flat surface
 462     grid%ht(i,j) = zsfc
 463     grid%phb(i,1,j) = grid%ht(i,j) * g
 464     grid%ph0(i,1,j) = grid%ht(i,j) * g
 465     grid%php(i,1,j) = 0.
 466   ENDDO
 467   ENDDO
 468
 469   DO J = jts, jte
 470   DO I = its, ite
 471
 472     p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
 473     grid%mub(i,j) = p_surf-grid%p_top
 474
 475 !  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
 476 !  interp theta (from interp) and compute 1/rho from eqn. of state
 477
 478     DO K = 1, kte-1
 479       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 480       grid%pb(i,k,j) = p_level
 481       grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
 482       grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 483     ENDDO
 484
 485 !  calc hydrostatic balance (alternatively we could interp the geopotential from the
 486 !  sounding, but this assures that the base state is in exact hydrostatic balance with
 487 !  respect to the model eqns.
 488
 489     DO k  = 2,kte
 490       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)
 491     ENDDO
 492
 493   ENDDO
 494   ENDDO
 495
 496   write(6,*) ' ptop is ',grid%p_top
 497   write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
 498
 499 !  calculate full state for each column - this includes moisture.
 500
 501   write(6,*) ' getting moist sounding for full state '
 502   dry_sounding = .false.
 503   CALL get_sounding( zsfc, zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 504
 505   DO J = jts, min(jde-1,jte)
 506   DO I = its, min(ide-1,ite)
 507
 508 !  At this point grid%p_top is already set. find the DRY mass in the column
 509 !  by interpolating the DRY pressure.
 510
 511    pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
 512
 513 !  compute the perturbation mass and the full mass
 514
 515     grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
 516     grid%mu_2(i,j) = grid%mu_1(i,j)
 517     grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
 518
 519 ! given the dry pressure and coordinate system, interp the potential
 520 ! temperature and qv
 521
 522     do k=1,kde-1
 523
 524       p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
 525
 526       moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
 527       grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
 528       grid%t_2(i,k,j)          = grid%t_1(i,k,j)
 529
 530
 531     enddo
 532
 533 !  integrate the hydrostatic equation (from the RHS of the bigstep
 534 !  vertical momentum equation) down from the top to get grid%p.
 535 !  first from the top of the model to the top pressure
 536
 537     k = kte-1  ! top level
 538
 539     qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
 540     qvf2 = 1./(1.+qvf1)
 541     qvf1 = qvf1*qvf2
 542
 543     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
 544     qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 545     grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 546                 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 547     grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 548
 549 !  down the column
 550
 551     do k=kte-2,1,-1
 552       qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
 553       qvf2 = 1./(1.+qvf1)
 554       qvf1 = qvf1*qvf2
 555       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)
 556       qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 557       grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 558                   (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 559       grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 560     enddo
 561
 562 !  this is the hydrostatic equation used in the model after the
 563 !  small timesteps.  In the model, grid%al (inverse density)
 564 !  is computed from the geopotential.
 565
 566
 567     grid%ph_1(i,1,j) = 0.
 568     DO k  = 2,kte
 569       grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
 570                    ((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)+ &
 571                     (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
 572
 573       grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
 574       grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
 575     ENDDO
 576
 577     if((i==2) .and. (j==2)) then
 578      write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
 579                               grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
 580                               grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
 581     endif
 582
 583   ENDDO
 584   ENDDO
 585
 586    write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
 587    write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
 588
 589    do k=1,kde-1
 590      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
 591                                       grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
 592                                       grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
 593    enddo
 594
 595    write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
 596    do k=1,kde-1
 597      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
 598                                       grid%p(1,k,1), grid%al(1,k,1), &
 599                                       grid%t_1(1,k,1), moist(1,k,1,P_QV)
 600    enddo
 601
 602 ! interp v
 603
 604   DO J = jts, jte
 605   DO I = its, min(ide-1,ite)
 606
 607     IF (j == jds) THEN
 608       z_at_v = grid%phb(i,1,j)/g
 609     ELSE IF (j == jde) THEN
 610       z_at_v = grid%phb(i,1,j-1)/g
 611     ELSE
 612       z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
 613     END IF
 614
 615     p_surf = interp_0( p_in, zk, z_at_v, nl_in )
 616
 617     DO K = 1, kte
 618       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 619       grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
 620       grid%v_2(i,k,j) = grid%v_1(i,k,j)
 621     ENDDO
 622
 623   ENDDO
 624   ENDDO
 625
 626 ! interp u
 627
 628   DO J = jts, min(jde-1,jte)
 629   DO I = its, ite
 630
 631     IF (i == ids) THEN
 632       z_at_u = grid%phb(i,1,j)/g
 633     ELSE IF (i == ide) THEN
 634       z_at_u = grid%phb(i-1,1,j)/g
 635     ELSE
 636       z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
 637     END IF
 638
 639     p_surf = interp_0( p_in, zk, z_at_u, nl_in )
 640
 641     DO K = 1, kte
 642       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 643       grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
 644       grid%u_2(i,k,j) = grid%u_1(i,k,j)
 645     ENDDO
 646
 647   ENDDO
 648   ENDDO
 649
 650 !  set w
 651
 652   DO J = jts, min(jde-1,jte)
 653   DO K = kts, kte
 654   DO I = its, min(ide-1,ite)
 655     grid%w_1(i,k,j) = 0.
 656     grid%w_2(i,k,j) = 0.
 657   ENDDO
 658   ENDDO
 659   ENDDO
 660
 661 !  set a few more things
 662
 663   DO J = jts, min(jde-1,jte)
 664   DO K = kts, kte-1
 665   DO I = its, min(ide-1,ite)
 666     grid%h_diabatic(i,k,j) = 0.
 667   ENDDO
 668   ENDDO
 669   ENDDO
 670
 671 ! Go ahead and initialize these from the sounding.  This will allow a run
 672 ! to actually succeed even if scm_force = 0
 673   DO k=1,kte-1
 674     grid%t_base(k) = grid%t_1(1,k,1)
 675     grid%qv_base(k) = moist(1,k,1,P_QV)
 676     grid%u_base(k) = grid%u_1(1,k,1)
 677     grid%v_base(k) = grid%v_1(1,k,1)
 678     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
 679   ENDDO
 680
 681       !  Save the dry perturbation potential temperature.
 682
 683       DO j = jts, min(jde-1,jte)
 684          DO k = kts, kte
 685             DO i = its, min(ide-1,ite)
 686                grid%th_phy_m_t0(i,k,j) = grid%t_2(i,k,j)
 687             END DO
 688          END DO
 689       END DO
 690
 691  !  Turn dry potential temperature into moist potential temperature
 692       !  at the very end of this routine
 693       !  This field will be in the model IC and and used to construct the
 694       !  BC file.
 695
 696       IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
 697       DO J  = jts, min(jde-1,jte)
 698          DO K = kts, kte-1
 699             DO I = its, min(ide-1,ite)
 700                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
 701             END DO
 702          END DO
 703       END DO
 704       ENDIF
 705
 706   RETURN
 707
 708  END SUBROUTINE init_domain_rk
 709
 710    SUBROUTINE init_module_initialize
 711    END SUBROUTINE init_module_initialize
 712
 713 !---------------------------------------------------------------------
 714
 715 !  test driver for get_sounding
 716 !
 717 !      implicit none
 718 !      integer n
 719 !      parameter(n = 1000)
 720 !      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
 721 !      logical dry
 722 !      integer nl,k
 723 !
 724 !      dry = .false.
 725 !      dry = .true.
 726 !      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
 727 !      write(6,*) ' input levels ',nl
 728 !      write(6,*) ' sounding '
 729 !      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
 730 !      do k=1,nl
 731 !        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)
 732 !      enddo
 733 !      end
 734 !
 735 !---------------------------------------------------------------------------
 736
 737       subroutine get_sounding( zsfc, zk, p, p_dry, theta, rho, &
 738                                u, v, qv, dry, nl_max, nl_in )
 739       implicit none
 740
 741       integer nl_max, nl_in
 742       real zsfc
 743       real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
 744            u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
 745       logical dry
 746
 747       integer n
 748       parameter(n=3000)
 749       logical debug
 750       parameter( debug = .true.)
 751
 752 ! input sounding data
 753
 754       real p_surf, th_surf, qv_surf
 755       real pi_surf, pi(n)
 756       real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
 757
 758 ! diagnostics
 759
 760       real rho_surf, p_input(n), rho_input(n)
 761       real pm_input(n)  !  this are for full moist sounding
 762
 763 ! local data
 764
 765       real r
 766       parameter (r = r_d)
 767       integer k, it, nl
 768       real qvf, qvf1, dz
 769
 770 !  first, read the sounding
 771
 772       call read_sounding( zsfc, p_surf, th_surf, qv_surf, &
 773                           h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
 774
 775       if(dry) then
 776        do k=1,nl
 777          qv_input(k) = 0.
 778        enddo
 779       endif
 780
 781       if(debug) write(6,*) ' number of input levels = ',nl
 782
 783         nl_in = nl
 784         if(nl_in .gt. nl_max ) then
 785           write(6,*) ' too many levels for input arrays ',nl_in,nl_max
 786           call wrf_error_fatal ( ' too many levels for input arrays ' )
 787         end if
 788
 789 !  compute diagnostics,
 790 !  first, convert qv(g/kg) to qv(g/g)
 791
 792       do k=1,nl
 793         qv_input(k) = 0.001*qv_input(k)
 794       enddo
 795
 796       p_surf = 100.*p_surf  ! convert to pascals
 797       qvf = 1. + rvovrd*qv_input(1)
 798       rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
 799       pi_surf = (p_surf/p1000mb)**(r/cp)
 800
 801       if(debug) then
 802         write(6,*) ' surface density is ',rho_surf
 803         write(6,*) ' surface pi is      ',pi_surf
 804       end if
 805
 806
 807 !  integrate moist sounding hydrostatically, starting from the
 808 !  specified surface pressure
 809 !  -> first, integrate from surface to lowest level
 810
 811           qvf = 1. + rvovrd*qv_input(1)
 812           qvf1 = 1. + qv_input(1)
 813           rho_input(1) = rho_surf
 814           dz = h_input(1)-zsfc
 815
 816 ! error check here
 817           if ( dz < 0.0 ) then
 818             write(6,*) "Your first input sounding level is below the WRF terrain elevation, aborting"
 819             stop "module_initialize_scm_xy:get_sounding"
 820           endif
 821           do it=1,10
 822             pm_input(1) = p_surf &
 823                     - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
 824             rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
 825           enddo
 826
 827 ! integrate up the column
 828
 829           do k=2,nl
 830             rho_input(k) = rho_input(k-1)
 831             dz = h_input(k)-h_input(k-1)
 832             qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
 833             qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
 834
 835             do it=1,10
 836               pm_input(k) = pm_input(k-1) &
 837                       - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
 838               rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
 839             enddo
 840           enddo
 841 !  we have the moist sounding
 842
 843 !  next, compute the dry sounding using p at the highest level from the
 844 !  moist sounding and integrating down.
 845
 846         p_input(nl) = pm_input(nl)
 847
 848           do k=nl-1,1,-1
 849             dz = h_input(k+1)-h_input(k)
 850             p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
 851           enddo
 852
 853
 854         do k=1,nl
 855           zk(k) = h_input(k)
 856           p(k) = pm_input(k)
 857           p_dry(k) = p_input(k)
 858           theta(k) = th_input(k)
 859           rho(k) = rho_input(k)
 860           u(k) = u_input(k)
 861           v(k) = v_input(k)
 862           qv(k) = qv_input(k)
 863
 864         enddo
 865
 866      if(debug) then
 867       write(6,*) ' sounding '
 868       write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
 869       do k=1,nl
 870         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)
 871       enddo
 872
 873      end if
 874
 875       end subroutine get_sounding
 876
 877 !-------------------------------------------------------
 878
 879       subroutine read_sounding( zsfc,ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
 880       implicit none
 881       integer n,nl
 882       real zsfc,ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
 883       real u10,v10,t2,q2
 884       logical end_of_file
 885       logical debug
 886
 887       integer k
 888
 889       open(unit=10,file='input_sounding',form='formatted',status='old')
 890       rewind(10)
 891
 892       read(10,*) zsfc, u10, v10, t2, q2, ps
 893       ps = ps/100.0
 894       ts = t2
 895       qvs = q2*1000
 896
 897       if(debug) then
 898         write(6,*) ' input sounding surface parameters '
 899         write(6,*) ' surface pressure (mb) ',ps
 900         write(6,*) ' surface pot. temp (K) ',ts
 901         write(6,*) ' surface mixing ratio (g/kg) ',qvs
 902       end if
 903
 904       end_of_file = .false.
 905       k = 0
 906
 907       do while (.not. end_of_file)
 908
 909         read(10,*,end=100) h(k+1), u(k+1), v(k+1), th(k+1), qv(k+1)
 910
 911         qv(k+1) = qv(k+1)*1000.0
 912         k = k+1
 913         if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
 914         go to 110
 915  100    end_of_file = .true.
 916  110    continue
 917       enddo
 918
 919       nl = k
 920
 921       close(unit=10,status = 'keep')
 922
 923       end subroutine read_sounding
 924
 925 !-------------------------------------------------------
 926
 927       subroutine read_soil( n,nl,tmn,tsk,zs,tslb,smois )
 928       implicit none
 929       integer n,nl
 930       real tmn,tsk
 931       real zs(n),tslb(n),smois(n)
 932       logical end_of_file
 933       logical debug
 934
 935       integer k
 936
 937       debug = .true.
 938
 939       open(unit=11,file='input_soil',form='formatted',status='old')
 940       rewind(11)
 941
 942       read(11,*) zs(1),tmn,tsk
 943
 944       if(debug) then
 945         write(6,*) ' input deep soil temperature (K) ',tmn
 946         write(6,*) ' input skin temperature (K) ',tsk
 947       end if
 948
 949       end_of_file = .false.
 950       k = 0
 951
 952       do while (.not. end_of_file)
 953
 954         read(11,*,end=100) zs(k+1), tslb(k+1), smois(k+1)
 955         k = k+1
 956         if(debug) write(6,'(1x,i3,3(1x,f16.7))') k, zs(k), tslb(k), smois(k)
 957         go to 110
 958  100    end_of_file = .true.
 959  110    continue
 960       enddo
 961
 962       nl = k
 963
 964       close(unit=11,status = 'keep')
 965
 966       end subroutine read_soil
 967
 968 END MODULE module_initialize_ideal