dyn_em/nest_init_utils.F

   1 ! careful adding any 1 stuff in here, adjust_tempqv has
   2 ! c3 and c4 to compute pressure with two reference pressures
   3
   4 SUBROUTINE init_domain_constants_em ( parent , nest )
   5    USE module_domain, ONLY : domain
   6    IMPLICIT NONE
   7    TYPE(domain)  :: parent , nest
   8
   9    INTEGER iswater, islake, isice, isurban, isoilwater, map_proj, julyr, julday
  10    REAL    truelat1 , truelat2 , gmt , moad_cen_lat , stand_lon, pole_lat, pole_lon
  11    CHARACTER (LEN=256) :: char_junk
  12
  13 ! single-value constants
  14
  15    nest%p_top   = parent%p_top
  16    nest%save_topo_from_real   = parent%save_topo_from_real
  17    nest%cfn     = parent%cfn
  18    nest%cfn1    = parent%cfn1
  19    nest%rdx     = 1./nest%dx
  20    nest%rdy     = 1./nest%dy
  21 !  nest%dts     = nest%dt/float(nest%time_step_sound)
  22    nest%dtseps  = parent%dtseps  ! used in height model only?
  23    nest%resm    = parent%resm    ! used in height model only?
  24    nest%zetatop = parent%zetatop ! used in height model only?
  25    nest%cf1     = parent%cf1
  26    nest%cf2     = parent%cf2
  27    nest%cf3     = parent%cf3
  28    nest%gmt     = parent%gmt
  29    nest%julyr   = parent%julyr
  30    nest%julday  = parent%julday
  31    nest%iswater = parent%iswater
  32    nest%isice   = parent%isice
  33    nest%isurban = parent%isurban
  34    nest%islake  = parent%islake
  35    nest%isoilwater = parent%isoilwater
  36    nest%mminlu  = trim(parent%mminlu)
  37    nest%tiso    = parent%tiso
  38    nest%tlp     = parent%tlp
  39    nest%p00     = parent%p00
  40    nest%t00     = parent%t00
  41    nest%tlp_strat= parent%tlp_strat
  42    nest%p_strat = parent%p_strat
  43 !cyl: variables for trajectory /float
  44    nest%traj_k = parent%traj_k
  45    nest%traj_long = parent%traj_long
  46    nest%traj_lat = parent%traj_lat
  47    nest%this_is_an_ideal_run = parent%this_is_an_ideal_run
  48    nest%lake_depth_flag = parent%lake_depth_flag
  49    nest%bathymetry_flag = parent%bathymetry_flag
  50
  51    CALL nl_get_mminlu ( 1, char_junk )
  52    CALL nl_get_iswater( 1, iswater )
  53    CALL nl_get_islake ( 1, islake )
  54    CALL nl_get_isice  ( 1, isice )
  55    CALL nl_get_isurban( 1, isurban )
  56    CALL nl_get_isoilwater(1, isoilwater )
  57    CALL nl_get_truelat1 ( 1 , truelat1 )
  58    CALL nl_get_truelat2 ( 1 , truelat2 )
  59    CALL nl_get_moad_cen_lat ( 1 , moad_cen_lat )
  60    CALL nl_get_stand_lon ( 1 , stand_lon )
  61    CALL nl_get_pole_lat ( 1 , pole_lat )
  62    CALL nl_get_pole_lon ( 1 , pole_lon )
  63    CALL nl_get_map_proj ( 1 , map_proj )
  64    CALL nl_get_gmt ( 1 , gmt)
  65    CALL nl_get_julyr ( 1 , julyr)
  66    CALL nl_get_julday ( 1 , julday)
  67    IF ( nest%id .NE. 1 ) THEN
  68      CALL nl_set_gmt (nest%id, gmt)
  69      CALL nl_set_julyr (nest%id, julyr)
  70      CALL nl_set_julday (nest%id, julday)
  71      CALL nl_set_iswater ( nest%id, iswater )
  72      CALL nl_set_islake  ( nest%id, islake )
  73      CALL nl_set_isice   ( nest%id, isice )
  74      CALL nl_set_isurban ( nest%id, isurban )
  75      CALL nl_set_isoilwater ( nest%id, isoilwater )
  76      CALL nl_set_mminlu ( nest%id, char_junk )
  77      CALL nl_set_truelat1 ( nest%id , truelat1 )
  78      CALL nl_set_truelat2 ( nest%id , truelat2 )
  79      CALL nl_set_moad_cen_lat ( nest%id , moad_cen_lat )
  80      CALL nl_set_stand_lon ( nest%id , stand_lon )
  81      CALL nl_set_pole_lat ( nest%id , pole_lat )
  82      CALL nl_set_pole_lon ( nest%id , pole_lon )
  83      CALL nl_set_map_proj ( nest%id , map_proj )
  84    END IF
  85    nest%gmt     = gmt
  86    nest%julday  = julday
  87    nest%julyr   = julyr
  88    nest%iswater = iswater
  89    nest%islake  = islake
  90    nest%isice   = isice
  91    nest%isoilwater = isoilwater
  92    nest%mminlu  = trim(char_junk)
  93    nest%truelat1= truelat1
  94    nest%truelat2= truelat2
  95    nest%moad_cen_lat= moad_cen_lat
  96    nest%stand_lon= stand_lon
  97    nest%pole_lat= pole_lat
  98    nest%pole_lon= pole_lon
  99    nest%map_proj= map_proj
 100
 101    nest%step_number  = parent%step_number
 102
 103 ! 1D constants (Z)
 104
 105 !DAVE - maybe test if vert_nest instead
 106 !  IF  ( parent%e_vert .EQ. nest%e_vert ) THEN
 107       nest%fnm(1:parent%e_vert)       = parent%fnm(1:parent%e_vert)
 108       nest%fnp(1:parent%e_vert)       = parent%fnp(1:parent%e_vert)
 109       nest%rdnw(1:parent%e_vert)      = parent%rdnw(1:parent%e_vert)
 110       nest%rdn(1:parent%e_vert)       = parent%rdn(1:parent%e_vert)
 111       nest%dnw(1:parent%e_vert)       = parent%dnw(1:parent%e_vert)
 112       nest%dn(1:parent%e_vert)        = parent%dn(1:parent%e_vert)
 113       nest%znu(1:parent%e_vert)       = parent%znu(1:parent%e_vert)
 114       nest%znw(1:parent%e_vert)       = parent%znw(1:parent%e_vert)
 115       nest%t_base(1:parent%e_vert)    = parent%t_base(1:parent%e_vert)
 116       nest%u_base(1:parent%e_vert)    = parent%u_base(1:parent%e_vert)
 117       nest%v_base(1:parent%e_vert)    = parent%v_base(1:parent%e_vert)
 118       nest%qv_base(1:parent%e_vert)   = parent%qv_base(1:parent%e_vert)
 119       nest%z_base(1:parent%e_vert)    = parent%z_base(1:parent%e_vert)
 120       nest%c1h(1:parent%e_vert)       = parent%c1h(1:parent%e_vert)
 121       nest%c2h(1:parent%e_vert)       = parent%c2h(1:parent%e_vert)
 122       nest%c3h(1:parent%e_vert)       = parent%c3h(1:parent%e_vert)
 123       nest%c4h(1:parent%e_vert)       = parent%c4h(1:parent%e_vert)
 124       nest%c1f(1:parent%e_vert)       = parent%c1f(1:parent%e_vert)
 125       nest%c2f(1:parent%e_vert)       = parent%c2f(1:parent%e_vert)
 126       nest%c3f(1:parent%e_vert)       = parent%c3f(1:parent%e_vert)
 127       nest%c4f(1:parent%e_vert)       = parent%c4f(1:parent%e_vert)
 128 !  END IF
 129
 130    nest%dzs       = parent%dzs
 131    nest%zs        = parent%zs
 132
 133 END SUBROUTINE init_domain_constants_em
 134
 135
 136 !---------------------------------------------------------------------------------------------------
 137
 138 SUBROUTINE init_domain_vert_nesting ( parent, nest, use_baseparam_fr_nml )
 139
 140 !KAL this is a driver to initialize the vertical coordinates for the nest when vertical nesting is used
 141    USE module_domain
 142    IMPLICIT NONE
 143    TYPE(domain), POINTER                      :: parent, nest
 144    LOGICAL                                    :: use_baseparam_fr_nml
 145    !local
 146    REAL, DIMENSION(parent%e_vert)             :: znw_c
 147
 148    INTERFACE
 149
 150       SUBROUTINE vert_cor_vertical_nesting_integer(nest,znw_c,k_dim_c)
 151          USE module_domain
 152          TYPE(domain), POINTER :: nest
 153          integer , intent(in) :: k_dim_c
 154          real , dimension(k_dim_c), INTENT(IN) :: znw_c
 155       END SUBROUTINE vert_cor_vertical_nesting_integer
 156
 157       SUBROUTINE vert_cor_vertical_nesting_arbitrary(nest,znw_c,kde_c,use_baseparam_fr_nml)
 158           USE module_domain
 159           TYPE(domain), POINTER :: nest
 160           INTEGER, INTENT(IN   ) :: kde_c
 161           REAL, DIMENSION(kde_c), INTENT(IN   ) :: znw_c
 162           LOGICAL, INTENT(IN   ) :: use_baseparam_fr_nml
 163       END SUBROUTINE vert_cor_vertical_nesting_arbitrary
 164
 165    END INTERFACE
 166
 167
 168    ! save the coarse grid values here
 169
 170    znw_c = nest%znw(1:parent%e_vert)
 171
 172    ! calculate the nest (fine) grid values here
 173    ! one of these calls goes to integer refinement in the vertical direction, and one goes to arbitrary refinement.  Eventually the call to integer refinement will be obsolete.
 174    if (nest%vert_refine_method .EQ. 1) then  !if you are in this subroutine there is vertical nesting- (i.e. nest%e_vert /= parent%e_vert to enter this subroutine)
 175       CALL vert_cor_vertical_nesting_integer(nest,znw_c,parent%e_vert)
 176    elseif (nest%vert_refine_method .EQ. 2) then
 177       CALL vert_cor_vertical_nesting_arbitrary(nest,znw_c,parent%e_vert,use_baseparam_fr_nml)
 178    endif
 179
 180 END SUBROUTINE init_domain_vert_nesting
 181
 182 !-----------------------------------------------------------------------------------------
 183
 184
 185 !this is a direct copy of a subroutine that is in ndown, but I couldn't link to the subroutine in ndown because it is compiled after this file
 186 !so a dependecy on ndown will not work. Additionally, ndown is not compiled for ideal cases. The variable is named parent in ndown, but it is actually operating on the nest. It has been renamed to nest here.
 187       SUBROUTINE vert_cor_vertical_nesting_integer(nest,znw_c,k_dim_c)
 188          USE module_domain
 189          USE module_model_constants
 190    IMPLICIT NONE
 191          TYPE(domain), POINTER :: nest
 192          integer , intent(in) :: k_dim_c
 193          real , dimension(k_dim_c), INTENT(IN) :: znw_c
 194
 195        integer :: kde_c , kde_n ,n_refine,ii,kkk,k
 196        real :: dznw_m,cof1,cof2
 197
 198        INTEGER :: ids, ide, jds, jde, kds, kde, &
 199                   ims, ime, jms, jme, kms, kme, &
 200                   its, ite, jts, jte, kts, kte
 201
 202 !KAL this subroutine recalculates the vertical coordinates for the nest when vertical nesting is used.  This routine is copied from ndown and allows integer refinement only.
 203
 204 !KAL znw is eta values on full w levels
 205 !KAL everything else is set from znw
 206 !KAL dnw is delta eta on w levels
 207 !KAL rdn is inverse delta eta on w levels
 208 !KAL fnp
 209
 210         kde_c = k_dim_c
 211         kde_n = nest%e_vert
 212 !        n_refine = nest%vert_refine_fact
 213         n_refine = (kde_n-1)/(kde_c-1)
 214
 215          kkk = 0
 216          do k = 1 , kde_c-1
 217          dznw_m = znw_c(k+1) - znw_c(k)
 218          do ii = 1,n_refine
 219          kkk = kkk + 1
 220          nest%znw(kkk) = znw_c(k) + float(ii-1)/float(n_refine)*dznw_m
 221          enddo
 222          enddo
 223          nest%znw(kde_n) = znw_c(kde_c)
 224          nest%znw(1) = znw_c(1)
 225
 226       DO k=1, kde_n-1
 227          nest%dnw(k) = nest%znw(k+1) - nest%znw(k)
 228          nest%rdnw(k) = 1./nest%dnw(k)
 229          nest%znu(k) = 0.5*(nest%znw(k+1)+nest%znw(k))
 230       END DO
 231
 232       DO k=2, kde_n-1
 233          nest%dn(k) = 0.5*(nest%dnw(k)+nest%dnw(k-1))
 234          nest%rdn(k) = 1./nest%dn(k)
 235          nest%fnp(k) = .5* nest%dnw(k  )/nest%dn(k)
 236          nest%fnm(k) = .5* nest%dnw(k-1)/nest%dn(k)
 237       END DO
 238
 239   cof1 = (2.*nest%dn(2)+nest%dn(3))/(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(2)
 240   cof2 =     nest%dn(2)        /(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(3)
 241
 242       nest%cf1  = nest%fnp(2) + cof1
 243       nest%cf2  = nest%fnm(2) - cof1 - cof2
 244       nest%cf3  = cof2
 245
 246       nest%cfn  = (.5*nest%dnw(kde_n-1)+nest%dn(kde_n-1))/nest%dn(kde_n-1)
 247       nest%cfn1 = -.5*nest%dnw(kde_n-1)/nest%dn(kde_n-1)
 248
 249       ! the variables dzs and zs are kept from the parent domain.  These are the depths and thickness of the soil layers, which are not included in vertical nesting.
 250
 251     CALL get_ijk_from_grid( nest, &
 252                             ids, ide, jds, jde, kds, kde, &
 253                             ims, ime, jms, jme, kms, kme, &
 254                             its, ite, jts, jte, kts, kte )
 255
 256     CALL compute_vcoord_1d_coeffs ( nest%ht, nest%etac, nest%znw, &
 257                                     model_config_rec%hybrid_opt, &
 258                                     r_d, g, p1000mb, &
 259                                     nest%p_top, nest%p00, nest%t00, nest%tlp, &
 260                                     ids, ide, jds, jde, kds, kde, &
 261                                     ims, ime, jms, jme, kms, kme, &
 262                                     its, ite, jts, jte, kts, kte, &
 263                                     nest%znu, &
 264                                     nest%c1f, nest%c2f, nest%c3f, nest%c4f, &
 265                                     nest%c1h, nest%c2h, nest%c3h, nest%c4h )
 266
 267       END SUBROUTINE vert_cor_vertical_nesting_integer
 268
 269 !-----------------------------------------------------------------------------------------
 270
 271 SUBROUTINE vert_cor_vertical_nesting_arbitrary(nest,znw_c,kde_c,use_baseparam_fr_nml)
 272     USE module_domain
 273     USE module_model_constants
 274     IMPLICIT NONE
 275     TYPE(domain), POINTER :: nest
 276     INTEGER, INTENT(IN   ) :: kde_c
 277     REAL, DIMENSION(kde_c), INTENT(IN   ) :: znw_c
 278     LOGICAL, INTENT(IN   ) :: use_baseparam_fr_nml
 279     INTEGER :: k, kde_n, ks, id
 280     REAL :: cof1, cof2
 281     REAL :: max_dz = 1000
 282     REAL :: p00, t00, a, tiso, a_strat, p_strat
 283     INTEGER :: ids, ide, jds, jde, kds, kde, &
 284                ims, ime, jms, jme, kms, kme, &
 285                its, ite, jts, jte, kts, kte
 286     REAL, DIMENSION(max_eta) :: eta_levels
 287
 288     IF ( use_baseparam_fr_nml ) then
 289       CALL nl_get_base_pres        ( 1 , p00        )
 290       CALL nl_get_base_temp        ( 1 , t00        )
 291       CALL nl_get_base_lapse       ( 1 , a          )
 292       CALL nl_get_iso_temp         ( 1 , tiso       )
 293       CALL nl_get_base_lapse_strat ( 1 , a_strat    )
 294       CALL nl_get_base_pres_strat  ( 1 , p_strat    )
 295       IF ((t00 .LT. 100.0) .OR. (p00 .LT. 10000.0)) THEN
 296         WRITE(wrf_err_message,*) '--- ERROR: bad base state for T00 or P00 in namelist.input file'
 297         CALL wrf_error_fatal(TRIM(wrf_err_message))
 298       END IF
 299     ELSE
 300       p00     = nest%p00
 301       t00     = nest%t00
 302       a       = nest%tlp
 303       tiso    = nest%tiso
 304       a_strat = nest%tlp_strat
 305       p_strat = nest%p_strat
 306       IF ((t00 .LT. 100.0) .OR. (p00 .LT. 10000.0)) THEN
 307         WRITE(wrf_err_message,*) '--- ERROR: did not find base state parameters in nest. Add use_baseparam_fr_nml = .true. in &dynamics and rerun'
 308         CALL wrf_error_fatal(TRIM(wrf_err_message))
 309       ENDIF
 310     ENDIF
 311
 312     kde_n = nest%e_vert
 313
 314     !DJW Added code for specifying multiple domains' eta_levels
 315     IF (nest%id .NE. 1) THEN
 316       id = 1
 317       ks = 1
 318       DO WHILE (nest%id .GT. id)
 319         id = id+1
 320         ks = ks+model_config_rec%e_vert(id-1)
 321       ENDDO
 322     ENDIF
 323     IF ((nest%this_is_an_ideal_run) .AND. (model_config_rec%eta_levels(1) .EQ. -1.0)) THEN
 324       !DJW If we're running an ideal case and do not have levels set in the
 325       !namelist then we set znw using constant spacing in eta
 326       DO k=1,kde_n
 327         CALL wrf_debug(0, "nest_init_utils: eta_levels are not specified in the namelist, setting levels with constant spacing in eta.")
 328         nest%znw(k) = 1.0-(k-1)/FLOAT((kde_n-1))
 329       ENDDO
 330     ELSEIF (.NOT.(nest%this_is_an_ideal_run) .AND. (model_config_rec%eta_levels(1) .EQ. -1.0)) THEN
 331       write(*,'(A,I2,A)') "--- WARNING: eta_levels are not specified in the namelist for grid_id=",nest%grid_id,", using WRF's default levels."
 332       CALL get_ijk_from_grid( nest, &
 333                               ids, ide, jds, jde, kds, kde, &
 334                               ims, ime, jms, jme, kms, kme, &
 335                               its, ite, jts, jte, kts, kte )
 336       write(*,'(A,F10.3)') "--- USING: nest%p_top   = ",nest%p_top
 337       write(*,'(A,F10.3)') "--- USING: g            = ",g
 338       write(*,'(A,F10.3)') "--- USING: cvpm         = ",cvpm
 339       write(*,'(A,F10.3)') "--- USING: r_d          = ",r_d
 340       write(*,'(A,F10.3)') "--- USING: cp           = ",cp
 341       write(*,'(A,F10.3)') "--- USING: p1000mb      = ",p1000mb
 342       write(*,'(A,F10.3)') "--- USING: t0           = ",t0
 343       write(*,'(A,F10.3)') "--- USING: p00          = ",p00
 344       write(*,'(A,F10.3)') "--- USING: t00          = ",t00
 345       write(*,'(A,F10.3)') "--- USING: a            = ",a
 346       write(*,'(A,F10.3)') "--- USING: tiso         = ",tiso
 347       write(*,'(A,F10.3)') "--- USING: a_strat      = ",a_strat
 348       write(*,'(A,F10.3)') "--- USING: p_strat      = ",p_strat
 349       CALL compute_eta ( nest%znw, &
 350                          eta_levels, max_eta, max_dz, &
 351                          nest%p_top, g, p00, cvpm, a, r_d, cp, &
 352                          t00, p1000mb, t0, tiso, p_strat, a_strat, &
 353                          ids, ide, jds, jde, kds, kde, &
 354                          ims, ime, jms, jme, kms, kme, &
 355                          its, ite, jts, jte, kts, kte )
 356     ELSE
 357       !DJW If we're in here then we are suppose to read in eta_levels
 358       !from the namelist. We do so and then check to make sure they make sense
 359     DO k=1,kde_n
 360       nest%znw(k) = model_config_rec%eta_levels(ks+k-1)
 361         write(*,'(A,I3,A,F6.3)') "nest%znw(",k,") = ",nest%znw(k)
 362     ENDDO
 363     !Check the value of the first and last eta level for our domain,
 364     !then check that the vector of eta levels is only decreasing
 365     IF (nest%znw(1) .NE. 1.0) THEN
 366         write(wrf_err_message,'(A,I2,A)') "--- ERROR: first eta_level for grid_id=",nest%grid_id," is not 1.0. Check namelist."
 367         CALL wrf_error_fatal( wrf_err_message )
 368     ENDIF
 369       write(*,'(A,F10.3)') "--- USING: g            = ",g
 370       write(*,'(A,F10.3)') "--- USING: cvpm         = ",cvpm
 371       write(*,'(A,F10.3)') "--- USING: r_d          = ",r_d
 372       write(*,'(A,F10.3)') "--- USING: cp           = ",cp
 373       write(*,'(A,F10.3)') "--- USING: p1000mb      = ",p1000mb
 374       write(*,'(A,F10.3)') "--- USING: t0           = ",t0
 375     IF (nest%znw(kde_n) .NE. 0.0) THEN
 376         write(wrf_err_message,'(A,I2,A)') "--- ERROR: last eta_level for grid_id=",nest%grid_id," is not 0.0. Check namelist."
 377         CALL wrf_error_fatal( wrf_err_message )
 378     ENDIF
 379     DO k=2,kde_n
 380       IF (nest%znw(k) .GT. nest%znw(k-1)) THEN
 381           write(wrf_err_message,'(A,I2,A)') "--- ERROR: eta_level for grid_id=",nest%grid_id," are not monotonically decreasing. Check namelist."
 382           CALL wrf_error_fatal( wrf_err_message )
 383       ENDIF
 384     ENDDO
 385     ENDIF
 386     !DJW End of added code for specifying eta_levels
 387
 388     DO k=1,kde_n-1
 389       nest%dnw(k) = nest%znw(k+1)-nest%znw(k)
 390       nest%rdnw(k) = 1./nest%dnw(k)
 391       nest%znu(k) = 0.5*(nest%znw(k+1)+nest%znw(k))
 392     ENDDO
 393     nest%znu(kde_n) = 0.0
 394
 395     DO k=2,kde_n-1
 396       nest%dn(k) = 0.5*(nest%dnw(k)+nest%dnw(k-1))
 397       nest%rdn(k) = 1./nest%dn(k)
 398       nest%fnp(k) = .5* nest%dnw(k  )/nest%dn(k)
 399       nest%fnm(k) = .5* nest%dnw(k-1)/nest%dn(k)
 400     ENDDO
 401
 402     cof1 = (2.*nest%dn(2)+nest%dn(3))/(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(2)
 403     cof2 =     nest%dn(2)            /(nest%dn(2)+nest%dn(3))*nest%dnw(1)/nest%dn(3)
 404
 405     nest%cf1  = nest%fnp(2) + cof1
 406     nest%cf2  = nest%fnm(2) - cof1 - cof2
 407     nest%cf3  = cof2
 408     nest%cfn  = (.5*nest%dnw(kde_n-1)+nest%dn(kde_n-1))/nest%dn(kde_n-1)
 409     nest%cfn1 = -.5*nest%dnw(kde_n-1)/nest%dn(kde_n-1)
 410
 411     CALL get_ijk_from_grid( nest, &
 412                             ids, ide, jds, jde, kds, kde, &
 413                             ims, ime, jms, jme, kms, kme, &
 414                             its, ite, jts, jte, kts, kte )
 415
 416     CALL compute_vcoord_1d_coeffs ( nest%ht, nest%etac, nest%znw, &
 417                                     model_config_rec%hybrid_opt, &
 418                                     r_d, g, p1000mb, &
 419                                     nest%p_top, nest%p00, nest%t00, nest%tlp, &
 420                                     ids, ide, jds, jde, kds, kde, &
 421                                     ims, ime, jms, jme, kms, kme, &
 422                                     its, ite, jts, jte, kts, kte, &
 423                                     nest%znu, &
 424                                     nest%c1f, nest%c2f, nest%c3f, nest%c4f, &
 425                                     nest%c1h, nest%c2h, nest%c3h, nest%c4h )
 426
 427 END SUBROUTINE vert_cor_vertical_nesting_arbitrary
 428
 429 SUBROUTINE compute_eta ( znw , &
 430                          eta_levels , max_eta , max_dz ,        &
 431                          p_top_def , g_def , p00_def ,          &
 432                          cvpm_def , a_def , r_d_def , cp_def ,  &
 433                          t00_def , p1000mb_def , t0_def ,       &
 434                          tiso_def , p_strat_def , a_strat_def , &
 435                          ids , ide , jds , jde , kds , kde ,    &
 436                          ims , ime , jms , jme , kms , kme ,    &
 437                          its , ite , jts , jte , kts , kte )
 438
 439       !  Compute eta levels, either using given values from the namelist (hardly
 440       !  a computation, yep, I know), or assuming a constant dz above the PBL,
 441       !  knowing p_top and the number of eta levels.
 442
 443       IMPLICIT NONE
 444
 445       INTEGER , INTENT(IN)        :: ids , ide , jds , jde , kds , kde , &
 446                                      ims , ime , jms , jme , kms , kme , &
 447                                      its , ite , jts , jte , kts , kte
 448       REAL , INTENT(IN)           :: max_dz
 449       REAL , INTENT(IN)           :: p_top_def , g_def , p00_def , cvpm_def , &
 450                                      a_def , r_d_def , cp_def , t00_def ,     &
 451                                      p1000mb_def , t0_def , tiso_def ,        &
 452                                      p_strat_def , a_strat_def
 453       INTEGER , INTENT(IN)        :: max_eta
 454       REAL , DIMENSION (max_eta)  :: eta_levels
 455
 456       REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw
 457
 458       !  Local vars
 459
 460       INTEGER :: k , kk
 461       REAL(KIND=8) :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
 462       REAL(KIND=8) , DIMENSION(kts:kte) :: dnw
 463       REAL(KIND=8) :: p_top , g , p00 , cvpm , &
 464                       a , r_d , cp , t00 ,     &
 465                       p1000mb , t0 , tiso ,        &
 466                       p_strat , a_strat
 467
 468       INTEGER , PARAMETER :: prac_levels = 59
 469       INTEGER :: loop , loop1
 470       REAL(KIND=8) , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
 471       REAL(KIND=8) , DIMENSION(MAX(prac_levels,kde)) :: alb , phb
 472       REAL(KIND=8) :: alb_max, t_init_max, pb_max, phb_max
 473
 474       CHARACTER(LEN=256) :: message
 475
 476       !  Compute top of the atmosphere with some silly levels.  We just want to
 477       !  integrate to get a reasonable value for ztop.  We use the planned
 478       !  PBL-esque
 479       !  levels, and then just coarse resolution above that.  We know p_top,
 480       !  and we
 481       !  have the base state vars.
 482
 483       p_top   = p_top_def
 484       g       = g_def
 485       p00     = p00_def
 486       cvpm    = cvpm_def
 487       a       = a_def
 488       r_d     = r_d_def
 489       cp      = cp_def
 490       t00     = t00_def
 491       p1000mb = p1000mb_def
 492       t0      = t0_def
 493       tiso    = tiso_def
 494       p_strat = p_strat_def
 495       a_strat = a_strat_def
 496
 497       p_surf = p00
 498
 499       znw_prac = (/ 1.0000_8 , 0.9930_8 , 0.9830_8 , 0.9700_8 , 0.9540_8 , 0.9340_8 , 0.9090_8 , 0.8800_8 , &
 500                     0.8500_8 , 0.8000_8 , 0.7500_8 , 0.7000_8 , 0.6500_8 , 0.6000_8 , 0.5500_8 , 0.5000_8 , &
 501                     0.4500_8 , 0.4000_8 , 0.3500_8 , 0.3000_8 , 0.2500_8 , 0.2000_8 , 0.1500_8 , 0.1000_8 , &
 502                     0.0800_8 , 0.0600_8 , 0.0400_8 , 0.0200_8 , &
 503                     0.0150_8 , 0.0100_8 , 0.0090_8 , 0.0080_8 , 0.0070_8 , 0.0060_8 , 0.0050_8 , 0.0040_8 , &
 504                     0.0035_8 , 0.0030_8 , &
 505                     0.0028_8 , 0.0026_8 , 0.0024_8 , 0.0022_8 , 0.0020_8 , &
 506                     0.0018_8 , 0.0016_8 , 0.0014_8 , 0.0012_8 , 0.0010_8 , &
 507                     0.0009_8 , 0.0008_8 , 0.0007_8 , 0.0006_8 , 0.0005_8 , 0.0004_8 , 0.0003_8 , &
 508                     0.0002_8 , 0.0001_8 , 0.00005_8, 0.0000_8 /)
 509
 510       DO k = 1 , prac_levels - 1
 511         znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5_8
 512         dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
 513       END DO
 514
 515       DO k = 1, prac_levels-1
 516         pb = znu_prac(k)*(p_surf - p_top) + p_top
 517         temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 518         IF ( pb .LT. p_strat ) THEN
 519           temp = tiso + A_strat*LOG(pb/p_strat)
 520         END IF
 521         t_init = temp*(p00/pb)**(r_d/cp) - t0
 522         alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 523       END DO
 524
 525       !  Base state mu is defined as base state surface pressure minus p_top
 526
 527       mub = p_surf - p_top
 528
 529       !  Integrate base geopotential, starting at terrain elevation.
 530
 531       phb(1) = 0._8
 532       DO k  = 2,prac_levels
 533         phb(k) = phb(k-1) - dnw_prac(k-1)*mub*alb(k-1)
 534       END DO
 535
 536       !  So, now we know the model top in meters.  Get the average depth above the PBL
 537       !  of each of the remaining levels.  We are going for a constant delta z thickness.
 538
 539       ztop     = phb(prac_levels) / g
 540       ztop_pbl = phb(8          ) / g
 541       dz = ( ztop - ztop_pbl ) / REAL ( kde - 8 )
 542
 543       IF ( dz .GE. max_dz ) THEN
 544         WRITE (message,FMT='("With a requested ",F7.1," Pa model top, the model lid will be about ",F7.1," m.")') p_top, ztop
 545         CALL wrf_message ( message )
 546         WRITE (message,FMT='("With ",I3," levels above the PBL, the level thickness will be about ",F6.1," m.")') kde-8, dz
 547         CALL wrf_message ( message )
 548         WRITE (message,FMT='("Thicknesses greater than ",F7.1," m are not recommended.")') max_dz
 549         CALL wrf_message ( message )
 550         CALL wrf_error_fatal ( 'Add more levels to namelist.input for e_vert' )
 551       END IF
 552
 553       !  Standard levels near the surface so no one gets in trouble.
 554
 555       DO k = 1 , 8
 556         eta_levels(k) = znw_prac(k)
 557       END DO
 558
 559       !  Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
 560       !  Skamarock et al, NCAR TN 468.  Use full levels, so
 561       !  use twice the thickness.
 562
 563       DO k = 8, kte-1-2
 564
 565         find_prac : DO kk = 1 , prac_levels
 566           IF (znw_prac(kk) .LT. eta_levels(k) ) THEN
 567             EXIT find_prac
 568           END IF
 569         end do find_prac
 570
 571         pb = 0.5*(eta_levels(k)+znw_prac(kk)) * (p_surf - p_top) + p_top
 572
 573         temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 574         IF ( pb .LT. p_strat ) THEN
 575           temp = tiso + A_strat * LOG ( pb/p_strat )
 576         END IF
 577 !       temp =             t00 + A*LOG(pb/p00)
 578         t_init = temp*(p00/pb)**(r_d/cp) - t0
 579         alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 580         eta_levels(k+1) = eta_levels(k) - dz*g / ( mub*alb(k) )
 581         pb = 0.5*(eta_levels(k)+eta_levels(k+1)) * (p_surf - p_top) + p_top
 582
 583         temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 584         IF ( pb .LT. p_strat ) THEN
 585           temp = tiso + A_strat * LOG ( pb/p_strat )
 586         END IF
 587 !       temp =             t00 + A*LOG(pb/p00)
 588         t_init = temp*(p00/pb)**(r_d/cp) - t0
 589         alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 590         eta_levels(k+1) = eta_levels(k) - dz*g / ( mub*alb(k) )
 591         pb = 0.5*(eta_levels(k)+eta_levels(k+1)) * (p_surf - p_top) + p_top
 592
 593         phb(k+1) = phb(k) - (eta_levels(k+1)-eta_levels(k)) * mub*alb(k)
 594       END DO
 595
 596       alb_max = alb(kte-1-2)
 597       t_init_max = t_init
 598       pb_max = pb
 599       phb_max = phb(kte-1)
 600
 601       DO k = 1 , kte-1-2
 602         znw(k) = eta_levels(k)
 603       END DO
 604       znw(kte-2) = 0.000
 605
 606       !  There is some iteration.  We want the top level, ztop, to be
 607       !  consistent with the delta z, and we want the half level values
 608       !  to be consistent with the eta levels.  The inner loop to 10 gets
 609       !  the eta levels very accurately, but has a residual at the top, due
 610       !  to dz changing.  We reset dz five times, and then things seem OK.
 611
 612       DO loop1 = 1 , 5
 613         DO loop = 1 , 10
 614           DO k = 8, kte-1-2-1
 615             pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
 616             temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 617             IF ( pb .LT. p_strat ) THEN
 618               temp = tiso + A_strat * LOG ( pb/p_strat )
 619             END IF
 620             t_init = temp*(p00/pb)**(r_d/cp) - t0
 621             alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 622             znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
 623           END DO
 624           pb = pb_max
 625           t_init = t_init_max
 626           alb(kte-1-2) = alb_max
 627           znw(kte-2) = znw(kte-1-2) - dz*g / ( mub*alb(kte-1-2) )
 628           IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
 629             print *,'Converged znw(kte) should be about 0.0 = ',znw(kte-2)
 630           END IF
 631           znw(kte-2) = 0.000
 632         END DO
 633
 634         !  Here is where we check the eta levels values we just computed.
 635
 636         DO k = 1, kde-1-2
 637           pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
 638           temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 639           IF ( pb .LT. p_strat ) THEN
 640             temp = tiso + A_strat * LOG ( pb/p_strat )
 641           END IF
 642           t_init = temp*(p00/pb)**(r_d/cp) - t0
 643           alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 644         END DO
 645
 646         phb(1) = 0.
 647         DO k  = 2,kde-2
 648           phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
 649         END DO
 650
 651         !  Reset the model top and the dz, and iterate.
 652
 653         ztop = phb(kde-2)/g
 654         ztop_pbl = phb(8)/g
 655         dz = ( ztop - ztop_pbl ) / REAL ( (kde-2) - 8 )
 656       END DO
 657
 658       IF ( dz .GT. max_dz ) THEN
 659         print *,'z (m)            = ',phb(1)/g
 660         do k = 2 ,kte-2
 661           print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
 662         end do
 663         print *,'dz (m) above fixed eta levels = ',dz
 664         print *,'namelist max_dz (m) = ',max_dz
 665         print *,'namelist p_top (Pa) = ',p_top
 666         CALL wrf_debug ( 0, 'You need one of three things:' )
 667         CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
 668         CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
 669         CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
 670         CALL wrf_debug ( 0, 'All are namelist options')
 671         CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
 672       END IF
 673
 674       !  Add those 2 levels back into the middle, just above the 8 levels
 675       !  that semi define a boundary layer.  After we open up the levels,
 676       !  then we just linearly interpolate in znw.  So now levels 1-8 are
 677       !  specified as the fixed boundary layer levels given in this routine.
 678       !  The top levels, 12 through kte are those computed.  The middle
 679       !  levels 9, 10, and 11 are equi-spaced in znw, and are each 1/2 the
 680       !  the znw thickness of levels 11 through 12.
 681
 682       DO k = kte-2 , 9 , -1
 683         znw(k+2) = znw(k)
 684       END DO
 685
 686       znw( 9) = 0.75 * znw( 8) + 0.25 * znw(12)
 687       znw(10) = 0.50 * znw( 8) + 0.50 * znw(12)
 688       znw(11) = 0.25 * znw( 8) + 0.75 * znw(12)
 689
 690       DO k = 8, kte-1
 691         pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
 692         temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
 693         IF ( pb .LT. p_strat ) THEN
 694           temp = tiso + A_strat * LOG ( pb/p_strat )
 695         END IF
 696         t_init = temp*(p00/pb)**(r_d/cp) - t0
 697         alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
 698         phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
 699       END DO
 700       phb(kte) = phb(kte-1) - (znw(kte)-znw(kte-1)) * mub*alb(kte-1)
 701
 702       k=1
 703       WRITE (*,FMT='("Full level index = ",I4,"     Height = ",F7.1," m")') k,phb(1)/g
 704       do k = 2 ,kte
 705         WRITE (*,FMT='("Full level index = ",I4,"     Height = ",F7.1," m      Thickness = ",F6.1," m")') k,phb(k)/g,(phb(k)-phb(k-1))/g
 706       end do
 707
 708 END SUBROUTINE compute_eta
 709
 710 !-----------------------------------------------------------------------------------------
 711
 712 SUBROUTINE blend_terrain ( ter_interpolated , ter_input , &
 713                            ids , ide , jds , jde , kds , kde , &
 714                            ims , ime , jms , jme , kms , kme , &
 715                            ips , ipe , jps , jpe , kps , kpe )
 716
 717    USE module_configure
 718    IMPLICIT NONE
 719
 720    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
 721                                                  ims , ime , jms , jme , kms , kme , &
 722                                                  ips , ipe , jps , jpe , kps , kpe
 723    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN)    :: ter_interpolated
 724    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: ter_input
 725
 726    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: ter_temp
 727    INTEGER :: i , j , k , spec_bdy_width
 728    REAL    :: r_blend_zones
 729    INTEGER blend_cell, blend_width
 730
 731    !  The fine grid elevation comes from the horizontally interpolated
 732    !  parent elevation for the first spec_bdy_width row/columns, so we need
 733    !  to get that value.  We blend the coarse and fine in the next blend_width
 734    !  rows and columns.  After that, in the interior, it is 100% fine grid.
 735
 736    CALL nl_get_spec_bdy_width ( 1, spec_bdy_width)
 737    CALL nl_get_blend_width ( 1, blend_width)
 738
 739    !  Initialize temp values to the nest ter elevation.  This fills in the values
 740    !  that will not be modified below.
 741
 742    DO j = jps , MIN(jpe, jde-1)
 743       DO k = kps , kpe
 744          DO i = ips , MIN(ipe, ide-1)
 745             ter_temp(i,k,j) = ter_input(i,k,j)
 746          END DO
 747       END DO
 748    END DO
 749
 750    !  To avoid some tricky indexing, we fill in the values inside out.  This allows
 751    !  us to overwrite incorrect assignments.  There are replicated assignments, and
 752    !  there is much unnecessary "IF test inside of a loop" stuff.  For a large
 753    !  domain, this is only a patch; for a small domain, this is not a biggy.
 754
 755    r_blend_zones = 1./(blend_width+1)
 756    DO j = jps , MIN(jpe, jde-1)
 757       DO k = kps , kpe
 758          DO i = ips , MIN(ipe, ide-1)
 759             DO blend_cell = blend_width,1,-1
 760                IF   ( ( i .EQ.       spec_bdy_width + blend_cell ) .OR.  ( j .EQ.       spec_bdy_width + blend_cell ) .OR. &
 761                       ( i .EQ. ide - spec_bdy_width - blend_cell ) .OR.  ( j .EQ. jde - spec_bdy_width - blend_cell ) ) THEN
 762                   ter_temp(i,k,j) = ( (blend_cell)*ter_input(i,k,j) + (blend_width+1-blend_cell)*ter_interpolated(i,k,j) ) &
 763                                     * r_blend_zones
 764                END IF
 765             ENDDO
 766             IF      ( ( i .LE.       spec_bdy_width     ) .OR.  ( j .LE.       spec_bdy_width     ) .OR. &
 767                       ( i .GE. ide - spec_bdy_width     ) .OR.  ( j .GE. jde - spec_bdy_width     ) ) THEN
 768                ter_temp(i,k,j) =      ter_interpolated(i,k,j)
 769             END IF
 770          END DO
 771       END DO
 772    END DO
 773
 774    !  Set nest elevation with temp values.  All values not overwritten in the above
 775    !  loops have been previously set in the initial assignment.
 776
 777    DO j = jps , MIN(jpe, jde-1)
 778       DO k = kps , kpe
 779          DO i = ips , MIN(ipe, ide-1)
 780             ter_input(i,k,j) = ter_temp(i,k,j)
 781          END DO
 782       END DO
 783    END DO
 784
 785 END SUBROUTINE blend_terrain
 786
 787 SUBROUTINE copy_3d_field ( ter_interpolated , ter_input , &
 788                            ids , ide , jds , jde , kds , kde , &
 789                            ims , ime , jms , jme , kms , kme , &
 790                            ips , ipe , jps , jpe , kps , kpe )
 791
 792    IMPLICIT NONE
 793
 794    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
 795                                                  ims , ime , jms , jme , kms , kme , &
 796                                                  ips , ipe , jps , jpe , kps , kpe
 797    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: ter_interpolated
 798    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN)  :: ter_input
 799
 800    INTEGER :: i , j , k
 801
 802    DO j = jps , MIN(jpe, jde-1)
 803       DO k = kps , kpe
 804          DO i = ips , MIN(ipe, ide-1)
 805             ter_interpolated(i,k,j) = ter_input(i,k,j)
 806          END DO
 807       END DO
 808    END DO
 809
 810 END SUBROUTINE copy_3d_field
 811
 812 SUBROUTINE adjust_tempqv ( mub, save_mub, c3, c4, znw, p_top, &
 813                            th, pp, qv,  &
 814                            use_theta_m, &
 815                            ids , ide , jds , jde , kds , kde , &
 816                            ims , ime , jms , jme , kms , kme , &
 817                            ips , ipe , jps , jpe , kps , kpe )
 818
 819    !USE module_configure
 820    !USE module_domain
 821    USE module_model_constants
 822
 823    !USE module_bc
 824    !USE module_io_domain
 825    !USE module_state_description
 826    !USE module_timing
 827    !USE module_soil_pre
 828    IMPLICIT NONE
 829
 830    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
 831                                                  ims , ime , jms , jme , kms , kme , &
 832                                                  ips , ipe , jps , jpe , kps , kpe
 833    REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: mub, save_mub
 834    REAL , DIMENSION(kms:kme) , INTENT(IN)    :: znw
 835    REAL , DIMENSION(kms:kme) , INTENT(IN)    :: c3, c4
 836    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: th, pp, qv
 837    INTEGER , INTENT(IN)                       :: use_theta_m
 838
 839    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: p_old, p_new, rh
 840    REAL :: es,dth,tc,e,dth1,thloc
 841    INTEGER :: i , j , k
 842
 843    real p_top
 844
 845
 846 ! p_old = full pressure before terrain blending; also compute initial RH
 847 ! which is going to be conserved during terrain blending
 848    DO j = jps , MIN(jpe, jde-1)
 849       DO k = kps , kpe-1
 850          DO i = ips , MIN(ipe, ide-1)
 851             p_old(i,k,j) = c4(k) +         c3(k)*save_mub(i,j) + p_top + pp(i,k,j)
 852             IF ( use_theta_m .EQ. 1 ) THEN
 853                tc = (th(i,k,j)+300.)*(p_old(i,k,j)/1.e5)**(2./7.)/(1.+R_v/R_d*qv(i,k,j)) - 273.15
 854             ELSE
 855                tc = (th(i,k,j)+300.)*(p_old(i,k,j)/1.e5)**(2./7.)                        - 273.15
 856             END IF
 857             es = 610.78*exp(17.0809*tc/(234.175+tc))
 858             e = qv(i,k,j)*p_old(i,k,j)/(0.622+qv(i,k,j))
 859             rh(i,k,j) = e/es
 860          END DO
 861       END DO
 862    END DO
 863 ! p_new = full pressure after terrain blending; also compute temperature correction and convert RH back to QV
 864    DO j = jps , MIN(jpe, jde-1)
 865       DO k = kps , kpe-1
 866          DO i = ips , MIN(ipe, ide-1)
 867             p_new(i,k,j) = c4(k) +         c3(k)*mub(i,j) + p_top + pp(i,k,j)
 868 ! 2*(g/cp-6.5e-3)*R_dry/g = -191.86e-3
 869             IF ( use_theta_m .EQ. 1 ) THEN
 870                thloc = (th(i,k,j)+300.)/(1.+R_v/R_d*qv(i,k,j))
 871             ELSE
 872                thloc = (th(i,k,j)+300.)
 873             END IF
 874             dth1 = -191.86e-3*(thloc         )/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
 875             dth  = -191.86e-3*(thloc+0.5*dth1)/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
 876             IF ( use_theta_m .EQ. 1 ) THEN
 877                th(i,k,j) = (thloc+dth)*(1.+R_v/R_d*qv(i,k,j))-300.
 878             ELSE
 879                th(i,k,j) = (thloc+dth)                       -300.
 880             END IF
 881             tc = (thloc+dth)*(p_new(i,k,j)/1.e5)**(2./7.) - 273.15
 882             es = 610.78*exp(17.0809*tc/(234.175+tc))
 883             e = rh(i,k,j)*es
 884             qv(i,k,j) = 0.622*e/(p_new(i,k,j)-e)
 885          END DO
 886       END DO
 887    END DO
 888
 889
 890 END SUBROUTINE adjust_tempqv
 891
 892 SUBROUTINE input_terrain_rsmas ( grid ,                        &
 893                            ids , ide , jds , jde , kds , kde , &
 894                            ims , ime , jms , jme , kms , kme , &
 895                            ips , ipe , jps , jpe , kps , kpe )
 896
 897    USE module_domain, ONLY : domain
 898    IMPLICIT NONE
 899    TYPE ( domain ) :: grid
 900
 901    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
 902                                                  ims , ime , jms , jme , kms , kme , &
 903                                                  ips , ipe , jps , jpe , kps , kpe
 904
 905    LOGICAL, EXTERNAL ::  wrf_dm_on_monitor
 906
 907    INTEGER :: i , j , k , myproc
 908    INTEGER, DIMENSION(256) :: ipath  ! array for integer coded ascii for passing path down to get_terrain
 909    CHARACTER*256 :: message, message2
 910    CHARACTER*256 :: rsmas_data_path
 911
 912 #if DM_PARALLEL
 913 ! Local globally sized arrays
 914    REAL , DIMENSION(ids:ide,jds:jde) :: ht_g, xlat_g, xlon_g
 915 #endif
 916
 917    CALL wrf_get_myproc ( myproc )
 918
 919 #if 0
 920 CALL domain_clock_get ( grid, current_timestr=message2 )
 921 WRITE ( message , FMT = '(A," HT before ",I3)' ) TRIM(message2), grid%id
 922 write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
 923 do j = jps,jpe
 924 do i = ips,ipe
 925 write(30+myproc,*)grid%ht(i,j)
 926 enddo
 927 enddo
 928 #endif
 929
 930    CALL nl_get_rsmas_data_path(1,rsmas_data_path)
 931    do i = 1, LEN(TRIM(rsmas_data_path))
 932       ipath(i) = ICHAR(rsmas_data_path(i:i))
 933    enddo
 934
 935 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
 936
 937    CALL wrf_patch_to_global_real ( grid%xlat , xlat_g , grid%domdesc, ' ' , 'xy' ,       &
 938                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
 939                                    ims, ime   , jms , jme   , 1 , 1 , &
 940                                    ips, ipe   , jps , jpe   , 1 , 1   )
 941    CALL wrf_patch_to_global_real ( grid%xlong , xlon_g , grid%domdesc, ' ' , 'xy' ,       &
 942                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
 943                                    ims, ime   , jms , jme   , 1 , 1 , &
 944                                    ips, ipe   , jps , jpe   , 1 , 1   )
 945
 946    IF ( wrf_dm_on_monitor() ) THEN
 947      CALL get_terrain ( grid%dx/1000., xlat_g(ids:ide,jds:jde), xlon_g(ids:ide,jds:jde), ht_g(ids:ide,jds:jde), &
 948                         ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
 949      WHERE ( ht_g(ids:ide,jds:jde) < -1000. ) ht_g(ids:ide,jds:jde) = 0.
 950    ENDIF
 951
 952    CALL wrf_global_to_patch_real ( ht_g , grid%ht , grid%domdesc, ' ' , 'xy' ,         &
 953                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
 954                                    ims, ime   , jms , jme   , 1 , 1 , &
 955                                    ips, ipe   , jps , jpe   , 1 , 1   )
 956 #else
 957
 958    CALL get_terrain ( grid%dx/1000., grid%xlat(ids:ide,jds:jde), grid%xlong(ids:ide,jds:jde), grid%ht(ids:ide,jds:jde), &
 959                        ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
 960    WHERE ( grid%ht(ids:ide,jds:jde) < -1000. ) grid%ht(ids:ide,jds:jde) = 0.
 961
 962 #endif
 963
 964 #if 0
 965 CALL domain_clock_get ( grid, current_timestr=message2 )
 966 WRITE ( message , FMT = '(A," HT after ",I3)' ) TRIM(message2), grid%id
 967 write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
 968 do j = jps,jpe
 969 do i = ips,ipe
 970 write(30+myproc,*)grid%ht(i,j)
 971 enddo
 972 enddo
 973 #endif
 974
 975 END SUBROUTINE input_terrain_rsmas
 976
 977 SUBROUTINE update_after_feedback_em ( grid  &
 978 !
 979 #include "dummy_new_args.inc"
 980 !
 981                  )
 982 !
 983 ! perform core specific updates, exchanges after
 984 ! model feedback  (called from med_feedback_domain) -John
 985 !
 986
 987 ! Driver layer modules
 988    USE module_domain, ONLY : domain, get_ijk_from_grid
 989    USE module_configure
 990    USE module_driver_constants
 991    USE module_machine
 992    USE module_tiles
 993 #ifdef DM_PARALLEL
 994    USE module_dm, ONLY : ntasks, ntasks_x, ntasks_y, itrace, local_communicator, mytask
 995    USE module_comm_dm, ONLY : HALO_EM_FEEDBACK_sub
 996 #else
 997    USE module_dm
 998 #endif
 999    USE module_bc
1000 ! Mediation layer modules
1001 ! Registry generated module
1002    USE module_state_description
1003
1004    IMPLICIT NONE
1005
1006    !  Subroutine interface block.
1007
1008    TYPE(domain) , TARGET         :: grid
1009
1010    !  Definitions of dummy arguments
1011 #include "dummy_new_decl.inc"
1012
1013    INTEGER                         :: ids , ide , jds , jde , kds , kde , &
1014                                       ims , ime , jms , jme , kms , kme , &
1015                                       ips , ipe , jps , jpe , kps , kpe
1016
1017   CALL wrf_debug( 500, "entering update_after_feedback_em" )
1018
1019 !  Obtain dimension information stored in the grid data structure.
1020   CALL get_ijk_from_grid (  grid ,                   &
1021                             ids, ide, jds, jde, kds, kde,    &
1022                             ims, ime, jms, jme, kms, kme,    &
1023                             ips, ipe, jps, jpe, kps, kpe    )
1024
1025   CALL wrf_debug( 500, "before HALO_EM_FEEDBACK.inc in update_after_feedback_em" )
1026 #ifdef DM_PARALLEL
1027 #include "HALO_EM_FEEDBACK.inc"
1028 #endif
1029   CALL wrf_debug( 500, "leaving update_after_feedback_em" )
1030
1031 END SUBROUTINE update_after_feedback_em
1032
1033 SUBROUTINE compute_vcoord_1d_coeffs ( ht, etac, znw, &
1034                                       hybrid_opt, &
1035                                       r_d, g, p1000mb, &
1036                                       p_top, p00, t00, a, &
1037                                       ids, ide, jds, jde, kds, kde, &
1038                                       ims, ime, jms, jme, kms, kme, &
1039                                       its, ite, jts, jte, kts, kte, &
1040                                       znu, &
1041                                       c1f, c2f, c3f, c4f, &
1042                                       c1h, c2h, c3h, c4h )
1043
1044    IMPLICIT NONE
1045
1046    !  Input
1047
1048    INTEGER ,                      INTENT(IN ) :: hybrid_opt
1049    REAL    ,                      INTENT(IN ) :: etac
1050    REAL    ,                      INTENT(IN ) :: r_d, g, p1000mb, p_top, p00, t00, a
1051    INTEGER ,                      INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1052                                                  ims, ime, jms, jme, kms, kme, &
1053                                                  its, ite, jts, jte, kts, kte
1054    REAL    , DIMENSION(kms:kme) ,           INTENT(IN ) :: znw
1055    REAL    , DIMENSION(ims:ime , jms:jme) , INTENT(IN ) :: ht
1056
1057    !  Output
1058
1059    REAL    , DIMENSION(kms:kme) , INTENT(OUT) :: znu
1060    REAL    , DIMENSION(kms:kme) , INTENT(OUT) :: c1f, c2f, c3f, c4f, &
1061                                                  c1h, c2h, c3h, c4h
1062
1063    !  Local vars
1064
1065    INTEGER :: i, j, k
1066    REAL    :: B1, B2, B3, B4, B5
1067    REAL    :: p_surf, press_above, press_below
1068    CHARACTER (LEN=256) :: a_message
1069
1070    DO k=kts, kte
1071       IF      ( hybrid_opt .EQ. 0 ) THEN
1072          c3f(k) = znw(k)
1073       ELSE IF ( hybrid_opt .EQ. 1 ) THEN
1074          c3f(k) = znw(k)
1075       ELSE IF ( hybrid_opt .EQ. 2 ) THEN
1076          B1 = 2. * etac**2 * ( 1. - etac )
1077          B2 = -etac * ( 4. - 3. * etac - etac**3 )
1078          B3 = 2. * ( 1. - etac**3 )
1079          B4 = - ( 1. - etac**2 )
1080          B5 = (1.-etac)**4
1081          c3f(k) = ( B1 + B2*znw(k) + B3*znw(k)**2 + B4*znw(k)**3 ) / B5
1082          IF ( znw(k) .LT. etac ) THEN
1083             c3f(k) = 0.
1084          END IF
1085          IF ( k .EQ. kds ) THEN
1086             c3f(k) = 1.
1087          ELSE IF ( k .EQ. kde ) THEN
1088             c3f(k) = 0.
1089          END IF
1090       ELSE IF ( hybrid_opt .EQ. 3 ) THEN
1091          c3f(k) = znw(k)*sin(0.5*3.14159*znw(k))**2
1092          IF      ( k .EQ. kds ) THEN
1093             c3f(k) = 1.
1094          ELSE IF ( k .EQ. kde ) THEN
1095             c3f(kde) = 0.
1096          END IF
1097       ELSE
1098          CALL wrf_message     ( 'ERROR: --- hybrid_opt' )
1099          CALL wrf_message     ( 'ERROR: --- hybrid_opt=0    ==> Standard WRF terrain-following coordinate' )
1100          CALL wrf_message     ( 'ERROR: --- hybrid_opt=1    ==> Standard WRF terrain-following coordinate, hybrid c1, c2, c3, c4' )
1101          CALL wrf_message     ( 'ERROR: --- hybrid_opt=2    ==> Hybrid, Klemp polynomial' )
1102          CALL wrf_message     ( 'ERROR: --- hybrid_opt=3    ==> Hybrid, sin^2' )
1103          CALL wrf_error_fatal ( 'ERROR: --- Invalid option' )
1104       END IF
1105    END DO
1106
1107    !  c4 is a function of c3 and eta.
1108
1109    DO k=1, kde
1110       c4f(k) = ( znw(k) - c3f(k) ) * ( p1000mb - p_top )
1111    ENDDO
1112
1113    !  Now on half levels, just add up and divide by 2 (for c3h).  Use (eta-c3)*(p00-pt) for c4 on half levels.
1114
1115    DO k=1, kde-1
1116       znu(k) = ( znw(k+1) + znw(k) ) * 0.5
1117       c3h(k) = ( c3f(k+1) + c3f(k) ) * 0.5
1118       c4h(k) = ( znu(k) - c3h(k) ) * ( p1000mb - p_top )
1119    ENDDO
1120
1121    !  c1 = d(B)/d(eta).  We define c1f as c1 on FULL levels.  For a vertical difference,
1122    !  we need to use B and eta on half levels.  The k-loop ends up referring to the
1123    !  full levels, neglecting the top and bottom.
1124
1125    DO k=kds+1, kde-1
1126       c1f(k) = ( c3h(k) - c3h(k-1) ) / ( znu(k) - znu(k-1) )
1127    ENDDO
1128
1129    !  The boundary conditions to get the coefficients:
1130    !  1) At k=kts: define d(B)/d(eta) = 1.  This gives us the same value of B and d(B)/d(eta)
1131    !     when doing the sigma-only B=eta.
1132    !  2) At k=kte: define d(B)/d(eta) = 0.  The curve B SMOOTHLY goes to zero, and at the very
1133    !     top, B continues to SMOOTHLY go to zero.  Note that for almost all cases of non B=eta,
1134    !     B is ALREADY=ZERO at the top, so this is a reasonable BC to assume.
1135
1136    c1f(kds) = 1.
1137    IF ( ( hybrid_opt .EQ. 0 ) .OR. ( hybrid_opt .EQ. 1 ) ) THEN
1138       c1f(kde) = 1.
1139    ELSE
1140       c1f(kde) = 0.
1141    END IF
1142
1143    !  c2 = ( 1. - c1(k) ) * (p00 - pt).  There is no vertical differencing, so we can do the
1144    !  full kds to kde looping.
1145
1146    DO k=kds, kde
1147       c2f(k) = ( 1. - c1f(k) ) * ( p1000mb - p_top )
1148    ENDDO
1149
1150    !  Now on half levels for c1 and c2.  The c1h will result from the full level c3 and full
1151    !  level eta differences.  The c2 value use the half level c1(k).
1152
1153    DO k=1, kde-1
1154       c1h(k) = ( c3f(k+1) - c3f(k) ) / ( znw(k+1) - znw(k) )
1155       c2h(k) = ( 1. - c1h(k) ) * ( p1000mb - p_top )
1156    ENDDO
1157
1158    !  With c3 and c4 on full levels, we can verify that we are maintaining monotonically
1159    !  decreasing reference pressure.
1160    !  P_dry_ref(k)  = c3f(k) * ( psurf - p_top ) + c4f(k) + p_top
1161
1162    DO j = jts, MIN(jde-1,jte)
1163       DO i = its, MIN(ide-1,ite)
1164          p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*ht(i,j)/a/r_d ) **0.5 )
1165          DO k=1, kde-1
1166             press_above = c3f(k+1)*(p_surf-p_top) + c4f(k+1) + p_top
1167             press_below = c3f(k  )*(p_surf-p_top) + c4f(k  ) + p_top
1168             IF ( press_below - press_above .LE. 0.0 ) THEN
1169                CALL wrf_message ( '----- ERROR: The reference pressure is not monotonically decreasing' )
1170                CALL wrf_message ( '             This tends to be caused by very high topography')
1171                WRITE (a_message,*) '(i,j) = ',i,j,', topography = ',ht(i,j),' m'
1172                CALL wrf_message ( TRIM(a_message) )
1173                WRITE (a_message,*) 'k = ',k,', reference pressure = ',press_below,' Pa'
1174                CALL wrf_message ( TRIM(a_message) )
1175                WRITE (a_message,*) 'k = ',k+1,', reference pressure = ',press_above,' Pa'
1176                CALL wrf_message ( TRIM(a_message) )
1177                WRITE (a_message,*) 'In the dynamics namelist record, reduce etac from ',etac
1178                CALL wrf_error_fatal ( TRIM(a_message) )
1179             END IF
1180          ENDDO
1181       ENDDO
1182    ENDDO
1183
1184 END SUBROUTINE compute_vcoord_1d_coeffs