dyn_em/module_initialize_tropical_cyclone.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    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    !  Local data
  82    INTEGER                             ::                       &
  83                                   ids, ide, jds, jde, kds, kde, &
  84                                   ims, ime, jms, jme, kms, kme, &
  85                                   its, ite, jts, jte, kts, kte, &
  86                                   i, j, k
  87
  88    INTEGER, PARAMETER :: nl_max = 1000
  89    REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
  90    INTEGER :: nl_in
  91
  92    INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc, lm
  93    REAL    :: B1, B2, B3, B4, B5
  94    REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
  95    REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
  96 !   REAL, EXTERNAL :: interp_0
  97    REAL    :: pi, rnd
  98
  99 ! variables/arrays for analytic vortex:
 100     integer :: nref,kref,nloop,i1,i2
 101     real :: r0,zdd,dd1,dd2,xref,vr,fcor,qvs,e1,tx,px,qx,ric,rjc,rr,diff,sst
 102     real*8 :: rmax,vmax,frac,angle
 103     real, dimension(:), allocatable :: rref,zref,th0,qv0,thv0,prs0,pi0,rh0
 104     real, dimension(:,:), allocatable :: vref,piref,pref,thref,thvref,qvref
 105     real :: pi_in,dz
 106
 107 !  stuff from original initialization that has been dropped from the Registry
 108    REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
 109    REAL    :: qvf1, qvf2, pd_surf
 110    INTEGER :: it
 111    real :: thtmp, ptmp, temp(3)
 112
 113    LOGICAL :: moisture_init
 114    LOGICAL :: stretch_grid, dry_sounding
 115    character (len=256) :: mminlu2
 116
 117 #ifdef DM_PARALLEL
 118 #    include "data_calls.inc"
 119 #endif
 120
 121
 122    SELECT CASE ( model_data_order )
 123          CASE ( DATA_ORDER_ZXY )
 124    kds = grid%sd31 ; kde = grid%ed31 ;
 125    ids = grid%sd32 ; ide = grid%ed32 ;
 126    jds = grid%sd33 ; jde = grid%ed33 ;
 127
 128    kms = grid%sm31 ; kme = grid%em31 ;
 129    ims = grid%sm32 ; ime = grid%em32 ;
 130    jms = grid%sm33 ; jme = grid%em33 ;
 131
 132    kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
 133    its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
 134    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 135          CASE ( DATA_ORDER_XYZ )
 136    ids = grid%sd31 ; ide = grid%ed31 ;
 137    jds = grid%sd32 ; jde = grid%ed32 ;
 138    kds = grid%sd33 ; kde = grid%ed33 ;
 139
 140    ims = grid%sm31 ; ime = grid%em31 ;
 141    jms = grid%sm32 ; jme = grid%em32 ;
 142    kms = grid%sm33 ; kme = grid%em33 ;
 143
 144    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 145    jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
 146    kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
 147          CASE ( DATA_ORDER_XZY )
 148    ids = grid%sd31 ; ide = grid%ed31 ;
 149    kds = grid%sd32 ; kde = grid%ed32 ;
 150    jds = grid%sd33 ; jde = grid%ed33 ;
 151
 152    ims = grid%sm31 ; ime = grid%em31 ;
 153    kms = grid%sm32 ; kme = grid%em32 ;
 154    jms = grid%sm33 ; jme = grid%em33 ;
 155
 156    its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
 157    kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
 158    jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
 159
 160    END SELECT
 161
 162 !-----------------------------------------------------------------------
 163 !               USER SETTINGS
 164
 165 !  Parameters for analytic vortex:
 166 !  Reference:  Rotunno and Emanuel, 1987, JAS, p. 549
 167
 168     r0     =   412500.0     ! outer radius (m)
 169     rmax   =    82500.0     ! approximate radius of max winds (m)
 170     vmax   =       15.0     ! approximate value of max wind speed (m/s)
 171     zdd    =    20000.0     ! depth of vortex (m)
 172
 173
 174 ! other settings:
 175
 176     fcor   =  5.0e-5        ! Coriolis parameter (1/s)
 177     sst    =  28.0          ! sea-surface temperature (Celsius)
 178
 179 !-----------------------------------------------------------------------
 180
 181    stretch_grid = .true.
 182    delt = 6.
 183 !   z_scale = .50
 184    z_scale = .40
 185    pi = 2.*asin(1.0)
 186    write(6,*) ' pi is ',pi
 187    nxc = (ide-ids)/2
 188    nyc = jde/2
 189    icm = ide/2
 190 ! lm is the half width of the land in terms of grid points
 191    lm = 25
 192    write(6,*) 'lm,icm-lm,icm+lm = ', lm,icm-lm,icm+lm
 193
 194    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
 195
 196 ! here we check to see if the boundary conditions are set properly
 197
 198    CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
 199
 200    moisture_init = .true.
 201
 202     grid%itimestep=0
 203
 204 #ifdef DM_PARALLEL
 205    CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
 206    CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
 207 #endif
 208
 209     mminlu2 = ' '
 210     mminlu2(1:4) = 'USGS'
 211     CALL nl_set_mminlu(1, mminlu2)
 212 !   CALL nl_set_mminlu(1, 'USGS')
 213     CALL nl_set_iswater(1,16)
 214     CALL nl_set_isice(1,3)
 215     CALL nl_set_cen_lat(1,20.)
 216     CALL nl_set_cen_lon(1,-105.)
 217     CALL nl_set_truelat1(1,0.)
 218     CALL nl_set_truelat2(1,0.)
 219     CALL nl_set_moad_cen_lat (1,0.)
 220     CALL nl_set_stand_lon (1,0.)
 221     CALL nl_set_pole_lon (1,0.)
 222     CALL nl_set_pole_lat (1,90.)
 223     CALL nl_set_map_proj(1,0)
 224 !   CALL model_to_grid_config_rec(1,model_config_rec,config_flags)
 225     CALL nl_get_iswater(1,grid%iswater)
 226
 227 !  here we initialize data that currently is not initialized
 228 !  in the input data
 229
 230
 231     DO j = jts, jte
 232       DO i = its, ite
 233          grid%ht(i,j)       = 0.
 234          grid%msft(i,j)     = 1.
 235          grid%msfu(i,j)     = 1.
 236          grid%msfv(i,j)     = 1.
 237          grid%msftx(i,j)    = 1.
 238          grid%msfty(i,j)    = 1.
 239          grid%msfux(i,j)    = 1.
 240          grid%msfuy(i,j)    = 1.
 241          grid%msfvx(i,j)    = 1.
 242          grid%msfvy(i,j)    = 1.
 243          grid%msfvx_inv(i,j)= 1.
 244          grid%sina(i,j)     = 0.
 245          grid%cosa(i,j)     = 1.
 246          grid%xlong(i,j)    = 0.0
 247          grid%e(i,j)        = 0.0
 248          grid%f(i,j)        = fcor
 249          grid%xlat(i,j)     = asin(0.5*fcor/EOMEG)/DEGRAD
 250 ! Hard-wire the ocean configuration
 251          grid%xland(i,j)     = 2.
 252          grid%lu_index(i,j)  = 16
 253          grid%tsk(i,j) = 273.15 + sst
 254          ! I think tmn is not used for ocean points, but set a value anyway:
 255          grid%tmn(i,j) = grid%tsk(i,j) - 10.0
 256       END DO
 257    END DO
 258
 259     print *,'   f = ',grid%f(its,jts)
 260     print *,'   lat = ',grid%xlat(its,jts)
 261
 262 ! for Noah LSM, additional variables need to be initialized
 263
 264    other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
 265
 266       CASE (SLABSCHEME)
 267
 268       CASE (LSMSCHEME)
 269
 270         DO j = jts , MIN(jde-1,jte)
 271            DO i = its , MIN(ide-1,ite)
 272               IF (grid%xland(i,j) .lt. 1.5) THEN
 273                  grid%vegfra(i,j) = 0.5
 274                  grid%canwat(i,j) = 0.
 275                  grid%ivgtyp(i,j) = 18
 276                  grid%isltyp(i,j) = 8
 277                  grid%xice(i,j) = 0.
 278                  grid%snow(i,j) = 0.
 279               ELSE
 280                  grid%vegfra(i,j) = 0.
 281                  grid%canwat(i,j) = 0.
 282                  grid%ivgtyp(i,j) = 16
 283                  grid%isltyp(i,j) = 14
 284                  grid%xice(i,j) = 0.
 285                  grid%snow(i,j) = 0.
 286               ENDIF
 287            END DO
 288         END DO
 289
 290       CASE (RUCLSMSCHEME)
 291
 292    END SELECT other_masked_fields
 293
 294    DO j = jts, jte
 295     DO k = kts, kte
 296       DO i = its, ite
 297          grid%ww(i,k,j)     = 0.
 298       END DO
 299     END DO
 300    END DO
 301
 302    grid%step_number = 0
 303
 304 ! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F
 305       CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, &
 306                      grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, &
 307                      grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, &
 308                      model_config_rec%sf_surface_physics(grid%id), &
 309                                    ids,ide, jds,jde, kds,kde,&
 310                                    ims,ime, jms,jme, kms,kme,&
 311                                    its,ite, jts,jte, kts,kte )
 312
 313 ! set up the grid
 314
 315    IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
 316      DO k=1, kde
 317       grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
 318                                 (1.-exp(-1./z_scale))
 319      ENDDO
 320    ELSE
 321      DO k=1, kde
 322       grid%znw(k) = 1. - float(k-1)/float(kde-1)
 323      ENDDO
 324    ENDIF
 325
 326    DO k=1, kde-1
 327     grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
 328     grid%rdnw(k) = 1./grid%dnw(k)
 329     grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
 330    ENDDO
 331
 332    DO k=2, kde-1
 333     grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
 334     grid%rdn(k) = 1./grid%dn(k)
 335     grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
 336     grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
 337    ENDDO
 338
 339    cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
 340    cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
 341    grid%cf1  = grid%fnp(2) + cof1
 342    grid%cf2  = grid%fnm(2) - cof1 - cof2
 343    grid%cf3  = cof2
 344
 345    grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
 346    grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
 347    grid%rdx = 1./config_flags%dx
 348    grid%rdy = 1./config_flags%dy
 349
 350 !  get the sounding from the ascii sounding file, first get dry sounding and
 351 !  calculate base state
 352
 353   write(6,*) ' getting dry sounding for base state '
 354   dry_sounding = .true.
 355   CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 356
 357   write(6,*) ' returned from reading sounding, nl_in is ',nl_in
 358
 359 !  find ptop for the desired ztop (ztop is input from the namelist),
 360 !  and find surface pressure
 361
 362   grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
 363
 364   !  For hybrid coord
 365
 366   DO k=kts, kte
 367      IF      ( config_flags%hybrid_opt .EQ. 0 ) THEN
 368         grid%c3f(k) = grid%znw(k)
 369      ELSE IF ( config_flags%hybrid_opt .EQ. 1 ) THEN
 370         grid%c3f(k) = grid%znw(k)
 371      ELSE IF ( config_flags%hybrid_opt .EQ. 2 ) THEN
 372         B1 = 2. * grid%etac**2 * ( 1. - grid%etac )
 373         B2 = -grid%etac * ( 4. - 3. * grid%etac - grid%etac**3 )
 374         B3 = 2. * ( 1. - grid%etac**3 )
 375         B4 = - ( 1. - grid%etac**2 )
 376         B5 = (1.-grid%etac)**4
 377         grid%c3f(k) = ( B1 + B2*grid%znw(k) + B3*grid%znw(k)**2 + B4*grid%znw(k)**3 ) / B5
 378         IF ( grid%znw(k) .LT. grid%etac ) THEN
 379            grid%c3f(k) = 0.
 380         END IF
 381         IF ( k .EQ. kds ) THEN
 382            grid%c3f(k) = 1.
 383         ELSE IF ( k .EQ. kde ) THEN
 384            grid%c3f(k) = 0.
 385         END IF
 386      ELSE IF ( config_flags%hybrid_opt .EQ. 3 ) THEN
 387         grid%c3f(k) = grid%znw(k)*sin(0.5*3.14159*grid%znw(k))**2
 388         IF      ( k .EQ. kds ) THEN
 389            grid%c3f(k) = 1.
 390         ELSE IF ( k .EQ. kds ) THEN
 391            grid%c3f(kde) = 0.
 392         END IF
 393      ELSE
 394         CALL wrf_message     ( 'ERROR: --- hybrid_opt' )
 395         CALL wrf_message     ( 'ERROR: --- hybrid_opt=0    ==> Standard WRF terrain-following coordinate' )
 396         CALL wrf_message     ( 'ERROR: --- hybrid_opt=1    ==> Standard WRF terrain-following coordinate, hybrid c1, c2, c3, c4' )
 397         CALL wrf_message     ( 'ERROR: --- hybrid_opt=2    ==> Hybrid, Klemp polynomial' )
 398         CALL wrf_message     ( 'ERROR: --- hybrid_opt=3    ==> Hybrid, sin^2' )
 399         CALL wrf_error_fatal ( 'ERROR: --- Invalid option' )
 400      END IF
 401   END DO
 402
 403   !  c4 is a function of c3 and eta.
 404
 405   DO k=1, kde
 406      grid%c4f(k) = ( grid%znw(k) - grid%c3f(k) ) * ( p1000mb - grid%p_top )
 407   ENDDO
 408
 409   !  Now on half levels, just add up and divide by 2 (for c3h).  Use (eta-c3)*(p00-pt) for c4 on half levels.
 410
 411   DO k=1, kde-1
 412      grid%znu(k) = ( grid%znw(k+1) + grid%znw(k) ) * 0.5
 413      grid%c3h(k) = ( grid%c3f(k+1) + grid%c3f(k) ) * 0.5
 414      grid%c4h(k) = ( grid%znu(k) - grid%c3h(k) ) * ( p1000mb - grid%p_top )
 415   ENDDO
 416
 417   !  c1 = d(B)/d(eta).  We define c1f as c1 on FULL levels.  For a vertical difference,
 418   !  we need to use B and eta on half levels.  The k-loop ends up referring to the
 419   !  full levels, neglecting the top and bottom.
 420
 421   DO k=kds+1, kde-1
 422      grid%c1f(k) = ( grid%c3h(k) - grid%c3h(k-1) ) / ( grid%znu(k) - grid%znu(k-1) )
 423   ENDDO
 424
 425   !  The boundary conditions to get the coefficients:
 426   !  1) At k=kts: define d(B)/d(eta) = 1.  This gives us the same value of B and d(B)/d(eta)
 427   !     when doing the sigma-only B=eta.
 428   !  2) At k=kte: define d(B)/d(eta) = 0.  The curve B SMOOTHLY goes to zero, and at the very
 429   !     top, B continues to SMOOTHLY go to zero.  Note that for almost all cases of non B=eta,
 430   !     B is ALREADY=ZERO at the top, so this is a reasonable BC to assume.
 431
 432   grid%c1f(kds) = 1.
 433   IF ( ( config_flags%hybrid_opt .EQ. 0 ) .OR. ( config_flags%hybrid_opt .EQ. 1 ) ) THEN
 434      grid%c1f(kde) = 1.
 435   ELSE
 436      grid%c1f(kde) = 0.
 437   END IF
 438
 439   !  c2 = ( 1. - c1(k) ) * (p00 - pt).  There is no vertical differencing, so we can do the
 440   !  full kds to kde looping.
 441
 442   DO k=kds, kde
 443      grid%c2f(k) = ( 1. - grid%c1f(k) ) * ( p1000mb - grid%p_top )
 444   ENDDO
 445
 446   !  Now on half levels for c1 and c2.  The c1h will result from the full level c3 and full
 447   !  level eta differences.  The c2 value use the half level c1(k).
 448
 449   DO k=1, kde-1
 450      grid%c1h(k) = ( grid%c3f(k+1) - grid%c3f(k) ) / ( grid%znw(k+1) - grid%znw(k) )
 451      grid%c2h(k) = ( 1. - grid%c1h(k) ) * ( p1000mb - grid%p_top )
 452   ENDDO
 453
 454   DO j=jts,jte
 455   DO i=its,ite  ! flat surface
 456     grid%phb(i,1,j) = 0.
 457     grid%php(i,1,j) = 0.
 458     grid%ph0(i,1,j) = 0.
 459     grid%ht(i,j) = 0.
 460   ENDDO
 461   ENDDO
 462
 463   DO J = jts, jte
 464   DO I = its, ite
 465
 466     p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
 467     grid%mub(i,j) = p_surf-grid%p_top
 468
 469 !  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
 470 !  interp theta (from interp) and compute 1/rho from eqn. of state
 471
 472     DO K = 1, kte-1
 473       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 474       grid%pb(i,k,j) = p_level
 475       grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
 476       grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
 477     ENDDO
 478
 479 !  calc hydrostatic balance (alternatively we could interp the geopotential from the
 480 !  sounding, but this assures that the base state is in exact hydrostatic balance with
 481 !  respect to the model eqns.
 482
 483     DO k  = 2,kte
 484       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)
 485     ENDDO
 486
 487   ENDDO
 488   ENDDO
 489
 490   write(6,*) ' ptop is ',grid%p_top
 491   write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
 492
 493 !  calculate full state for each column - this includes moisture.
 494
 495   write(6,*) ' getting moist sounding for full state '
 496   dry_sounding = .false.
 497   CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
 498
 499   DO J = jts, min(jde-1,jte)
 500   DO I = its, min(ide-1,ite)
 501
 502 !  At this point grid%p_top is already set. find the DRY mass in the column
 503 !  by interpolating the DRY pressure.
 504
 505    pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
 506
 507 !  compute the perturbation mass and the full mass
 508
 509     grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
 510     grid%mu_2(i,j) = grid%mu_1(i,j)
 511     grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
 512
 513 ! given the dry pressure and coordinate system, interp the potential
 514 ! temperature and qv
 515
 516     do k=1,kde-1
 517
 518       p_level = grid%c3h(k)*(pd_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 519
 520       moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
 521       grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
 522       grid%t_2(i,k,j)          = grid%t_1(i,k,j)
 523
 524
 525     enddo
 526
 527 !  integrate the hydrostatic equation (from the RHS of the bigstep
 528 !  vertical momentum equation) down from the top to get grid%p.
 529 !  first from the top of the model to the top pressure
 530
 531     k = kte-1  ! top level
 532
 533     qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
 534     qvf2 = 1./(1.+qvf1)
 535     qvf1 = qvf1*qvf2
 536
 537     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
 538     qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 539     grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 540                 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 541     grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 542
 543 !  down the column
 544
 545     do k=kte-2,1,-1
 546       qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
 547       qvf2 = 1./(1.+qvf1)
 548       qvf1 = qvf1*qvf2
 549       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)
 550       qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 551       grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 552                   (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 553       grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 554     enddo
 555
 556 !  this is the hydrostatic equation used in the model after the
 557 !  small timesteps.  In the model, grid%al (inverse density)
 558 !  is computed from the geopotential.
 559
 560
 561     grid%ph_1(i,1,j) = 0.
 562     DO k  = 2,kte
 563       grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
 564                    ((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)+ &
 565                     (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
 566
 567       grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
 568       grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
 569     ENDDO
 570
 571     if((i==2) .and. (j==2)) then
 572      write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
 573                               grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
 574                               grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
 575     endif
 576
 577   ENDDO
 578   ENDDO
 579
 580 !-----------------------------------------------------------------------
 581 !  Analytic vortex.
 582 !  Reference:  Rotunno and Emanuel, 1987, JAS, p. 549
 583
 584     dd2 = 2.0 * rmax / ( r0 + rmax )
 585
 586     nref = 1 + int( float(ide-ids+1)/2.0 )
 587     kref = kte-1
 588
 589     print *,'  ids,ide,kds,kds  = ',ids,ide,kds,kde
 590     print *,'  its,ite,kts,kts  = ',its,ite,kts,kte
 591     print *,'  nref,fcor        = ',nref,fcor
 592     print *,'  r0,rmax,vmax,zdd = ',r0,rmax,vmax,zdd
 593
 594     allocate(  rref(nref)         )
 595     allocate(  zref(0:kref+1)     )
 596     allocate(   th0(0:kref+1)     )
 597     allocate(   qv0(0:kref+1)     )
 598     allocate(  thv0(0:kref+1)     )
 599     allocate(  prs0(0:kref+1)     )
 600     allocate(   pi0(0:kref+1)     )
 601     allocate(   rh0(0:kref+1)     )
 602     allocate(  vref(nref,0:kref+1))
 603     allocate( piref(nref,0:kref+1))
 604     allocate(  pref(nref,0:kref+1))
 605     allocate( thref(nref,0:kref+1))
 606     allocate(thvref(nref,0:kref+1))
 607     allocate( qvref(nref,0:kref+1))
 608
 609     ! get base state:
 610     print *,'  zref,th0,qv0,thv0:'
 611     do k=1,kref
 612       th0(k) = t0+grid%t_1(1,k,1)
 613       qv0(k) = moist(1,k,1,P_QV)
 614       thv0(k) = th0(k)*(1.0+(r_v/r_d)*qv0(k))/(1.0+qv0(k))
 615       zref(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
 616       print *,k,zref(k),th0(k),qv0(k),thv0(k)
 617     enddo
 618
 619     print *,'  prs0,pi0,rh0:'
 620     do k=1,kref
 621       prs0(k) = grid%p(1,k,1)+grid%pb(1,k,1)
 622       pi0(k) = (prs0(k)/p0)**(r_d/cp)
 623       E1=1000.0*SVP1*EXP(SVP2*(th0(k)*pi0(k)-SVPT0)/(th0(k)*pi0(k)-SVP3))
 624       qvs = EP_2*E1/(prs0(k)-E1)
 625       rh0(k) = qv0(k)/qvs
 626       print *,k,prs0(k),pi0(k),rh0(k)
 627     enddo
 628
 629     zref(0) = -zref(1)
 630     zref(kref+1) = zref(kref)+(zref(kref)-zref(kref-1))
 631
 632       rref=0.0
 633       vref=0.0
 634      piref=0.0
 635       pref=0.0
 636      thref=0.0
 637     thvref=0.0
 638      qvref=0.0
 639
 640     do i=1,nref
 641       rref(i) = config_flags%dx*(float(i-1)+0.5)
 642     enddo
 643
 644     print *,'  zref,dz:'
 645     do k=0,kref+1
 646       if( k.ge.2 .and. k.le.kref )then
 647         print *,k,zref(k),zref(k)-zref(k-1)
 648       else
 649         print *,k,zref(k)
 650       endif
 651     enddo
 652
 653     print *,'  vref:'
 654     do k=1,kref
 655       do i=1,nref
 656         if(rref(i).lt.r0)then
 657           dd1 = 2.0 * rmax / ( rref(i) + rmax )
 658           vr = sqrt( vmax**2 * (rref(i)/rmax)**2     &
 659           * ( dd1 ** 3 - dd2 ** 3 ) + 0.25*fcor*fcor*rref(i)*rref(i) )   &
 660                   - 0.5 * fcor * rref(i)
 661         else
 662           vr = 0.0
 663         endif
 664         if(zref(k).lt.zdd)then
 665           vref(i,k) = vr * (zdd-zref(k))/(zdd-0.0)
 666         else
 667           vref(i,k) = 0.0
 668         endif
 669         if(k.eq.1) print *,i,rref(i),vref(i,k)
 670       enddo
 671     enddo
 672
 673     print *,'  Iterate:'
 674     DO nloop=1,20
 675
 676       ! get qv and thv from rh and th:
 677       do k=1,kref
 678       do i=1,nref
 679         tx = (pi0(k)+piref(i,k))*(th0(k)+thref(i,k))
 680         px = p0*((pi0(k)+piref(i,k))**(cp/r_d))
 681         E1 = 1000.0*SVP1*EXP(SVP2*(tx-SVPT0)/(tx-SVP3))
 682         qvs = EP_2*E1/(px-E1)
 683         qvref(i,k) = rh0(k)*qvs
 684         thvref(i,k)=(th0(k)+thref(i,k))*(1.0+(r_v/r_d)*qvref(i,k))   &
 685                                            /(1.0+qvref(i,k))
 686       enddo
 687       enddo
 688
 689       ! get nondimensional pressure perturbation (piref):
 690       do k=1,kref
 691         piref(nref,k)=0.0
 692         do i=nref,2,-1
 693           piref(i-1,k) = piref(i,k)                                       &
 694        + (rref(i-1)-rref(i))/(cp*0.5*(thvref(i-1,k)+thvref(i,k))) * 0.5 * &
 695            ( vref(i  ,k)*vref(i  ,k)/rref(i)                              &
 696             +vref(i-1,k)*vref(i-1,k)/rref(i-1)                            &
 697              + fcor * ( vref(i,k) + vref(i-1,k) ) )
 698         enddo
 699       enddo
 700
 701       do i=1,nref
 702         piref(i,   0) = piref(i, 1)
 703         piref(i,kref+1) = piref(i,kref)
 704       enddo
 705
 706       ! get potential temperature perturbation (thref):
 707       do k=2,kref
 708       do i=1,nref
 709         thref(i,k) = 0.5*( cp*0.5*(thvref(i,k)+thvref(i,k+1))*(piref(i,k+1)-piref(i,k))/(zref(k+1)-zref(k))     &
 710                           +cp*0.5*(thvref(i,k)+thvref(i,k-1))*(piref(i,k)-piref(i,k-1))/(zref(k)-zref(k-1)) )   &
 711                         *thv0(k)/g
 712         thref(i,k)=(thv0(k)+thref(i,k))*(1.0+qvref(i,k))/(1.0+(r_v/r_d)*qvref(i,k))-th0(k)
 713       enddo
 714       enddo
 715
 716       k=1
 717       do i=1,nref
 718         thref(i,k) = ( cp*0.5*(thvref(i,k)+thvref(i,k+1))*(piref(i,k+1)-piref(i,k))/(zref(k+1)-zref(k)) )   &
 719                         *thv0(k)/g
 720         thref(i,k)=(thv0(k)+thref(i,k))*(1.0+qvref(i,k))/(1.0+(r_v/r_d)*qvref(i,k))-th0(k)
 721       enddo
 722
 723       print *,'  th,qv,pi = ',nloop,thref(1,1),qvref(1,1),piref(1,1)
 724
 725     ENDDO   ! enddo for iteration
 726
 727     ! reference (total) pressure:
 728     do k=1,kref
 729     do i=1,nref
 730       pref(i,k) = p0*( ( pi0(k)+piref(i,k) )**(cp/r_d) )
 731     enddo
 732     enddo
 733
 734     ! analytic axisymmetric vortex is ready ... now interpolate to 3D grid:
 735     ! (note:  vortex is placed in center of domain)
 736
 737     ric = float(ide-ids+1)/2.0
 738     rjc = float(jde-jds+1)/2.0
 739
 740     print *,'  ids,ide,jds,jde = ',ids,ide,jds,jde
 741     print *,'  ric,rjc = ',ric,rjc
 742
 743     print *,'  zk:'
 744     do k=1,kte
 745       zk(k) = zref(k)
 746       print *,k,zk(k)
 747     enddo
 748
 749     nl_in = kte-1
 750     print *,'  nl_in = ',nl_in
 751
 752     DO J = jts, min(jde-1,jte)
 753     DO I = its, min(ide-1,ite)
 754       rr = sqrt( ( (float(i)-ric)*config_flags%dx )**2 + ( (float(j)-rjc)*config_flags%dy )**2 )
 755       rr = min( rr , rref(nref) )
 756             diff = -1.0e20
 757             ii = 0
 758             do while( diff.lt.0.0 )
 759               ii = ii + 1
 760               diff = rref(ii)-rr
 761             enddo
 762             i2 = max( ii , 2 )
 763             i1 = i2-1
 764             frac = (      rr-rref(i1))   &
 765                   /(rref(i2)-rref(i1))
 766             do k=1,nl_in
 767               pi_in    = pi0(k)+piref(i1,k)+(piref(i2,k)-piref(i1,k))*frac
 768               qv(k)    =        qvref(i1,k)+(qvref(i2,k)-qvref(i1,k))*frac
 769               theta(k) = th0(k)+thref(i1,k)+(thref(i2,k)-thref(i1,k))*frac
 770               p_in(k) = p1000mb*(pi_in**(cp/r_d))
 771               qvf = 1. + rvovrd*qv(k)
 772               rho(k) = 1./((r_d/p1000mb)*theta(k)*qvf*((p_in(k)/p1000mb)**cvpm))
 773             enddo
 774
 775             pd_in(nl_in) = p_in(nl_in)
 776
 777             do k=nl_in-1,1,-1
 778               dz = zk(k+1)-zk(k)
 779               pd_in(k) = pd_in(k+1) + 0.5*dz*(rho(k)+rho(k+1))*g
 780             enddo
 781
 782 !  At this point grid%p_top is already set. find the DRY mass in the column
 783 !  by interpolating the DRY pressure.
 784
 785    pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
 786
 787 !  compute the perturbation mass and the full mass
 788
 789     grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
 790     grid%mu_2(i,j) = grid%mu_1(i,j)
 791     grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
 792
 793 ! given the dry pressure and coordinate system, interp the potential
 794 ! temperature and qv
 795
 796     do k=1,kde-1
 797
 798       p_level = grid%c3h(k)*(pd_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 799
 800       moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
 801       grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
 802       grid%t_2(i,k,j)          = grid%t_1(i,k,j)
 803
 804
 805     enddo
 806
 807
 808
 809 !  integrate the hydrostatic equation (from the RHS of the bigstep
 810 !  vertical momentum equation) down from the top to get grid%p.
 811 !  first from the top of the model to the top pressure
 812
 813     k = kte-1  ! top level
 814
 815     qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
 816     qvf2 = 1./(1.+qvf1)
 817     qvf1 = qvf1*qvf2
 818
 819     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
 820     qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 821     grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 822                 (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 823     grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 824
 825 !  down the column
 826
 827     do k=kte-2,1,-1
 828       qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
 829       qvf2 = 1./(1.+qvf1)
 830       qvf1 = qvf1*qvf2
 831       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)
 832       qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 833       grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 834                   (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 835       grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
 836     enddo
 837
 838 !  this is the hydrostatic equation used in the model after the
 839 !  small timesteps.  In the model, grid%al (inverse density)
 840 !  is computed from the geopotential.
 841
 842
 843     grid%ph_1(i,1,j) = 0.
 844     DO k  = 2,kte
 845       grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
 846                    (grid%c1h(k-1)*(grid%mub(i,j)+grid%mu_1(i,j))+grid%c2h(k-1))*grid%al(i,k-1,j)+ &
 847                     (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
 848
 849       grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
 850       grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
 851     ENDDO
 852
 853     ENDDO   ! do loop for i
 854     ENDDO   ! do loop for j
 855
 856 !-------------------------------------------
 857 ! Done with mass fields, now get winds:
 858
 859 ! interp v
 860
 861   DO J = jts, jte
 862   DO I = its, min(ide-1,ite)
 863
 864       rr = sqrt( ( (float(i)-ric)*config_flags%dx )**2 + ( (float(j)-0.5-rjc)*config_flags%dy )**2 )
 865       rr = min( rr , rref(nref) )
 866             diff = -1.0e20
 867             ii = 0
 868             do while( diff.lt.0.0 )
 869               ii = ii + 1
 870               diff = rref(ii)-rr
 871             enddo
 872             i2 = max( ii , 2 )
 873             i1 = i2-1
 874             frac = (      rr-rref(i1))   &
 875                   /(rref(i2)-rref(i1))
 876             angle = datan2(dble( (float(j)-0.5-rjc)*config_flags%dy ),   &
 877                            dble( (float(i)-ric)*config_flags%dx ) )
 878             do k=1,kte
 879               v(k) = (vref(i1,k)+( vref(i2,k)- vref(i1,k))*frac )*cos(angle)
 880               p_in(k) = pref(i1,k)+(pref(i2,k)-pref(i1,k))*frac
 881             enddo
 882
 883     IF (j == jds) THEN
 884       z_at_v = grid%phb(i,1,j)/g
 885     ELSE IF (j == jde) THEN
 886       z_at_v = grid%phb(i,1,j-1)/g
 887     ELSE
 888       z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
 889     END IF
 890
 891     p_surf = interp_0( p_in, zk, z_at_v, nl_in )
 892
 893     DO K = 1, kte
 894       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 895       grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
 896       grid%v_2(i,k,j) = grid%v_1(i,k,j)
 897     ENDDO
 898
 899   ENDDO
 900   ENDDO
 901
 902 ! interp u
 903
 904   DO J = jts, min(jde-1,jte)
 905   DO I = its, ite
 906
 907       rr = sqrt( ( (float(i)-ric-0.5)*config_flags%dx )**2 + ( (float(j)-rjc)*config_flags%dy )**2 )
 908       rr = min( rr , rref(nref) )
 909             diff = -1.0e20
 910             ii = 0
 911             do while( diff.lt.0.0 )
 912               ii = ii + 1
 913               diff = rref(ii)-rr
 914             enddo
 915             i2 = max( ii , 2 )
 916             i1 = i2-1
 917             frac = (      rr-rref(i1))   &
 918                   /(rref(i2)-rref(i1))
 919             angle = datan2(dble( (float(j)-rjc)*config_flags%dy ),   &
 920                            dble( (float(i)-0.5-ric)*config_flags%dx ) )
 921             do k=1,kte
 922               u(k) = -(vref(i1,k)+( vref(i2,k)- vref(i1,k))*frac )*sin(angle)
 923               p_in(k) = pref(i1,k)+(pref(i2,k)-pref(i1,k))*frac
 924             enddo
 925
 926     IF (i == ids) THEN
 927       z_at_u = grid%phb(i,1,j)/g
 928     ELSE IF (i == ide) THEN
 929       z_at_u = grid%phb(i-1,1,j)/g
 930     ELSE
 931       z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
 932     END IF
 933
 934     p_surf = interp_0( p_in, zk, z_at_u, nl_in )
 935
 936     DO K = 1, kte
 937       p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
 938       grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
 939       grid%u_2(i,k,j) = grid%u_1(i,k,j)
 940     ENDDO
 941
 942   ENDDO
 943   ENDDO
 944
 945     ! All done ... deallocate arrays:
 946
 947     deallocate(  rref )
 948     deallocate(  zref )
 949     deallocate(   th0 )
 950     deallocate(   qv0 )
 951     deallocate(  thv0 )
 952     deallocate(  prs0 )
 953     deallocate(   pi0 )
 954     deallocate(   rh0 )
 955     deallocate(  vref )
 956     deallocate( piref )
 957     deallocate(  pref )
 958     deallocate( thref )
 959     deallocate(thvref )
 960     deallocate( qvref )
 961
 962     print *,'  completed vortex init successfully '
 963
 964 !-----------------------------------------------------------------------
 965
 966    if (0.gt.1) then
 967 !#if 0
 968 !  The tropical_cyclone case is adapted from the squall line case
 969 !  so we just turn off the thermal perturbation
 970
 971 !  thermal perturbation to kick off convection
 972         call random_seed
 973   write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
 974   write(6,*) ' delt for perturbation ',delt
 975
 976   DO J = jts, min(jde-1,jte)
 977 !    yrad = config_flags%dy*float(j-nyc)/4000.
 978      yrad = 0.
 979     DO I = its, min(ide-1,ite)
 980       xrad = config_flags%dx*float(i-nxc)/10000.
 981 !     xrad = 0.
 982       DO K = 1, 35
 983
 984 !  put in preturbation theta (bubble) and recalc density.  note,
 985 !  the mass in the column is not changing, so when theta changes,
 986 !  we recompute density and geopotential
 987         zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j)  &
 988                    +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
 989         zrad = (zrad-1500.)/1500.
 990         RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
 991 !        IF(RAD <= 1.) THEN
 992         call RANDOM_NUMBER(rnd)
 993           grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*(rnd-0.5)
 994          !  grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
 995            grid%t_2(i,k,j)=grid%t_1(i,k,j)
 996            qvf = 1. + rvovrd*moist(i,k,j,P_QV)
 997            grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
 998                         (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
 999            grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
1000 !        ENDIF
1001       ENDDO
1002
1003 !  rebalance hydrostatically
1004
1005       DO k  = 2,kte
1006         grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (grid%dnw(k-1))*(       &
1007                      (grid%c1h(k-1)*(grid%mub(i,j)+grid%mu_1(i,j))+grid%c2h(k-1))*grid%al(i,k-1,j)+ &
1008                       (grid%c1h(k-1)*grid%mu_1(i,j))*grid%alb(i,k-1,j)  )
1009
1010         grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
1011         grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
1012       ENDDO
1013
1014     ENDDO
1015   ENDDO
1016         endif
1017 !#endif
1018
1019    write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
1020    write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
1021    do k=1,kde-1
1022      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
1023                                       grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
1024                                       grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
1025    enddo
1026
1027    write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
1028    do k=1,kde-1
1029      write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
1030                                       grid%p(1,k,1), grid%al(1,k,1), &
1031                                       grid%t_1(1,k,1), moist(1,k,1,P_QV)
1032    enddo
1033
1034 !! interp v
1035 !
1036 !  DO J = jts, jte
1037 !  DO I = its, min(ide-1,ite)
1038 !
1039 !    IF (j == jds) THEN
1040 !      z_at_v = grid%phb(i,1,j)/g
1041 !    ELSE IF (j == jde) THEN
1042 !      z_at_v = grid%phb(i,1,j-1)/g
1043 !    ELSE
1044 !      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
1045 !    END IF
1046 !
1047 !    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
1048 !
1049 !    DO K = 1, kte
1050 !      p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
1051 !      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
1052 !      grid%v_2(i,k,j) = grid%v_1(i,k,j)
1053 !    ENDDO
1054 !
1055 !  ENDDO
1056 !  ENDDO
1057 !
1058 !! interp u
1059 !
1060 !  DO J = jts, min(jde-1,jte)
1061 !  DO I = its, ite
1062 !
1063 !    IF (i == ids) THEN
1064 !      z_at_u = grid%phb(i,1,j)/g
1065 !    ELSE IF (i == ide) THEN
1066 !      z_at_u = grid%phb(i-1,1,j)/g
1067 !    ELSE
1068 !      z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
1069 !    END IF
1070 !
1071 !    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
1072 !
1073 !    DO K = 1, kte
1074 !      p_level = grid%c3h(k)*(p_surf - grid%p_top)+grid%c4h(k) + grid%p_top
1075 !      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
1076 !      grid%u_2(i,k,j) = grid%u_1(i,k,j)
1077 !    ENDDO
1078 !
1079 !  ENDDO
1080 !  ENDDO
1081
1082 !  set w
1083
1084   DO J = jts, min(jde-1,jte)
1085   DO K = kts, kte
1086   DO I = its, min(ide-1,ite)
1087     grid%w_1(i,k,j) = 0.
1088     grid%w_2(i,k,j) = 0.
1089   ENDDO
1090   ENDDO
1091   ENDDO
1092
1093 !  set a few more things
1094
1095   DO J = jts, min(jde-1,jte)
1096   DO K = kts, kte-1
1097   DO I = its, min(ide-1,ite)
1098     grid%h_diabatic(i,k,j) = 0.
1099   ENDDO
1100   ENDDO
1101   ENDDO
1102
1103   DO k=1,kte-1
1104     grid%t_base(k) = grid%t_1(1,k,1)
1105     grid%qv_base(k) = moist(1,k,1,P_QV)
1106     grid%u_base(k) = grid%u_1(1,k,1)
1107     grid%v_base(k) = grid%v_1(1,k,1)
1108     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
1109   ENDDO
1110
1111   DO J = jts, min(jde-1,jte)
1112   DO I = its, min(ide-1,ite)
1113      thtmp   = grid%t_2(i,1,j)+t0
1114      ptmp    = grid%p(i,1,j)+grid%pb(i,1,j)
1115      temp(1) = thtmp * (ptmp/p1000mb)**rcp
1116      thtmp   = grid%t_2(i,2,j)+t0
1117      ptmp    = grid%p(i,2,j)+grid%pb(i,2,j)
1118      temp(2) = thtmp * (ptmp/p1000mb)**rcp
1119      thtmp   = grid%t_2(i,3,j)+t0
1120      ptmp    = grid%p(i,3,j)+grid%pb(i,3,j)
1121      temp(3) = thtmp * (ptmp/p1000mb)**rcp
1122
1123 !     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
1124      ! I don't know why this is here, so I have commented it out:
1125 !!!     grid%tmn(I,J)=grid%tsk(I,J)-0.5
1126   ENDDO
1127   ENDDO
1128
1129       !  Save the dry perturbation potential temperature.
1130
1131       DO j = jts, min(jde-1,jte)
1132          DO k = kts, kte
1133             DO i = its, min(ide-1,ite)
1134                grid%th_phy_m_t0(i,k,j) = grid%t_2(i,k,j)
1135             END DO
1136          END DO
1137       END DO
1138
1139   !  Turn dry potential temperature into moist potential temperature
1140       !  at the very end of this routine
1141       !  This field will be in the model IC and and used to construct the
1142       !  BC file.
1143
1144       IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
1145       DO J  = jts, min(jde-1,jte)
1146          DO K = kts, kte-1
1147             DO I = its, min(ide-1,ite)
1148                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
1149             END DO
1150          END DO
1151       END DO
1152       ENDIF
1153
1154   RETURN
1155
1156  END SUBROUTINE init_domain_rk
1157
1158    SUBROUTINE init_module_initialize
1159    END SUBROUTINE init_module_initialize
1160
1161 !---------------------------------------------------------------------
1162
1163 !  test driver for get_sounding
1164 !
1165 !      implicit none
1166 !      integer n
1167 !      parameter(n = 1000)
1168 !      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
1169 !      logical dry
1170 !      integer nl,k
1171 !
1172 !      dry = .false.
1173 !      dry = .true.
1174 !      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
1175 !      write(6,*) ' input levels ',nl
1176 !      write(6,*) ' sounding '
1177 !      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
1178 !      do k=1,nl
1179 !        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)
1180 !      enddo
1181 !      end
1182 !
1183 !---------------------------------------------------------------------------
1184
1185       subroutine get_sounding( zk, p, p_dry, theta, rho, &
1186                                u, v, qv, dry, nl_max, nl_in )
1187       implicit none
1188
1189       integer nl_max, nl_in
1190       real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
1191            u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
1192       logical dry
1193
1194       integer n
1195       parameter(n=3000)
1196       logical debug
1197       parameter( debug = .true.)
1198
1199 ! input sounding data
1200
1201       real p_surf, th_surf, qv_surf
1202       real pi_surf, pi(n)
1203       real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
1204
1205 ! diagnostics
1206
1207       real rho_surf, p_input(n), rho_input(n)
1208       real pm_input(n)  !  this are for full moist sounding
1209
1210 ! local data
1211
1212       real r
1213       parameter (r = r_d)
1214       integer k, it, nl
1215       real qvf, qvf1, dz
1216
1217 !  first, read the sounding
1218
1219       call read_sounding( p_surf, th_surf, qv_surf, &
1220                           h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
1221
1222       if(dry) then
1223        do k=1,nl
1224          qv_input(k) = 0.
1225        enddo
1226       endif
1227
1228       if(debug) write(6,*) ' number of input levels = ',nl
1229
1230         nl_in = nl
1231         if(nl_in .gt. nl_max ) then
1232           write(6,*) ' too many levels for input arrays ',nl_in,nl_max
1233           call wrf_error_fatal ( ' too many levels for input arrays ' )
1234         end if
1235
1236 !  compute diagnostics,
1237 !  first, convert qv(g/kg) to qv(g/g)
1238
1239       do k=1,nl
1240         qv_input(k) = 0.001*qv_input(k)
1241       enddo
1242
1243       p_surf = 100.*p_surf  ! convert to pascals
1244       qvf = 1. + rvovrd*qv_input(1)
1245       rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
1246       pi_surf = (p_surf/p1000mb)**(r/cp)
1247
1248       if(debug) then
1249         write(6,*) ' surface density is ',rho_surf
1250         write(6,*) ' surface pi is      ',pi_surf
1251       end if
1252
1253
1254 !  integrate moist sounding hydrostatically, starting from the
1255 !  specified surface pressure
1256 !  -> first, integrate from surface to lowest level
1257
1258           qvf = 1. + rvovrd*qv_input(1)
1259           qvf1 = 1. + qv_input(1)
1260           rho_input(1) = rho_surf
1261           dz = h_input(1)
1262           do it=1,10
1263             pm_input(1) = p_surf &
1264                     - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
1265             rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
1266           enddo
1267
1268 ! integrate up the column
1269
1270           do k=2,nl
1271             rho_input(k) = rho_input(k-1)
1272             dz = h_input(k)-h_input(k-1)
1273             qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
1274             qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
1275
1276             do it=1,10
1277               pm_input(k) = pm_input(k-1) &
1278                       - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
1279               rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
1280             enddo
1281           enddo
1282
1283 !  we have the moist sounding
1284
1285 !  next, compute the dry sounding using p at the highest level from the
1286 !  moist sounding and integrating down.
1287
1288         p_input(nl) = pm_input(nl)
1289
1290           do k=nl-1,1,-1
1291             dz = h_input(k+1)-h_input(k)
1292             p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
1293           enddo
1294
1295
1296         do k=1,nl
1297           zk(k) = h_input(k)
1298           p(k) = pm_input(k)
1299           p_dry(k) = p_input(k)
1300           theta(k) = th_input(k)
1301           rho(k) = rho_input(k)
1302           u(k) = u_input(k)
1303           v(k) = v_input(k)
1304           qv(k) = qv_input(k)
1305
1306         enddo
1307
1308      if(debug) then
1309       write(6,*) ' sounding '
1310       write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
1311       do k=1,nl
1312         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)
1313       enddo
1314
1315      end if
1316
1317       end subroutine get_sounding
1318
1319 !-------------------------------------------------------
1320
1321       subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
1322       implicit none
1323       integer n,nl
1324       real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
1325       logical end_of_file
1326       logical debug
1327
1328       integer k
1329
1330       open(unit=10,file='input_sounding',form='formatted',status='old')
1331       rewind(10)
1332       read(10,*) ps, ts, qvs
1333       if(debug) then
1334         write(6,*) ' input sounding surface parameters '
1335         write(6,*) ' surface pressure (mb) ',ps
1336         write(6,*) ' surface pot. temp (K) ',ts
1337         write(6,*) ' surface mixing ratio (g/kg) ',qvs
1338       end if
1339
1340       end_of_file = .false.
1341       k = 0
1342
1343       do while (.not. end_of_file)
1344
1345         read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
1346         k = k+1
1347         if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
1348         go to 110
1349  100    end_of_file = .true.
1350  110    continue
1351       enddo
1352
1353       nl = k
1354
1355       close(unit=10,status = 'keep')
1356
1357       end subroutine read_sounding
1358
1359 END MODULE module_initialize_ideal