!REAL:MODEL_LAYER:INITIALIZATION
! This MODULE holds the routines which are used to perform various initializations
! for the individual domains, specifically for the Eulerian, mass-based coordinate.
!-----------------------------------------------------------------------
MODULE module_initialize 6
USE module_bc
USE module_configure
USE module_domain
USE module_io_domain
USE module_model_constants
USE module_si_io_em
USE module_state_description
USE module_timing
USE module_soil_pre
#ifdef DM_PARALLEL
USE module_dm
#endif
CONTAINS
!-------------------------------------------------------------------
SUBROUTINE init_domain ( grid ) 4,21
IMPLICIT NONE
! Input space and data. No gridded meteorological data has been stored, though.
! TYPE (domain), POINTER :: grid
TYPE (domain) :: grid
! Local data.
INTEGER :: dyn_opt
INTEGER :: idum1, idum2
#ifdef DEREF_KLUDGE
INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
#endif
#ifdef DEREF_KLUDGE
sm31 = grid%sm31
em31 = grid%em31
sm32 = grid%sm32
em32 = grid%em32
sm33 = grid%sm33
em33 = grid%em33
#endif
CALL get_dyn_opt ( dyn_opt )
CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
IF ( dyn_opt .eq. 1 &
.or. dyn_opt .eq. 2 &
.or. dyn_opt .eq. 3 &
) THEN
CALL init_domain_rk( grid, &
!
#include <em_actual_args.inc>
!
)
ELSE
WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
STOP 'ERROR-dyn_opt-wrong-in-namelist'
ENDIF
END SUBROUTINE init_domain
!-------------------------------------------------------------------
SUBROUTINE init_domain_rk ( grid, & 7,143
!
# include <em_dummy_args.inc>
!
)
USE module_optional_si_input
IMPLICIT NONE
! Input space and data. No gridded meteorological data has been stored, though.
! TYPE (domain), POINTER :: grid
TYPE (domain) :: grid
# include <em_dummy_decl.inc>
TYPE (grid_config_rec_type) :: config_flags
! Local domain indices and counters.
INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
INTEGER :: &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
i, j, k
! Local data
INTEGER :: error
REAL :: p_surf, p_level
REAL :: cof1, cof2
REAL :: qvf , qvf1 , qvf2 , pd_surf
REAL :: p00 , t00 , a
REAL :: hold_znw
LOGICAL :: stretch_grid, dry_sounding, debug
! INTEGER , PARAMETER :: nl_max = 1000
! REAL , DIMENSION(nl_max) :: dn
integer::oops1,oops2
#define COPY_IN
#include <em_scalar_derefs.inc>
#ifdef DM_PARALLEL
# include <em_data_calls.inc>
#endif
SELECT CASE ( model_data_order )
CASE ( DATA_ORDER_ZXY )
kds = grid%sd31 ; kde = grid%ed31 ;
ids = grid%sd32 ; ide = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
kms = grid%sm31 ; kme = grid%em31 ;
ims = grid%sm32 ; ime = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XYZ )
ids = grid%sd31 ; ide = grid%ed31 ;
jds = grid%sd32 ; jde = grid%ed32 ;
kds = grid%sd33 ; kde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
jms = grid%sm32 ; jme = grid%em32 ;
kms = grid%sm33 ; kme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XZY )
ids = grid%sd31 ; ide = grid%ed31 ;
kds = grid%sd32 ; kde = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
kms = grid%sm32 ; kme = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
END SELECT
CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
! Check to see if the boundary conditions are set properly in the namelist file.
! This checks for sufficiency and redundancy.
CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
! Some sort of "this is the first time" initialization. Who knows.
step_number = 0
grid%itimestep=0
! Pull in the info in the namelist to compare it to the input data.
real_data_init_type = model_config_rec%real_data_init_type
! Take the data from the input file and store it in the variables that
! use the WRF naming and ordering conventions.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
IF ( snow(i,j) .GE. 10. ) then
snowc(i,j) = 1.
ELSE
snowc(i,j) = 0.0
END IF
END DO
END DO
! Set flag integers for presence of snowh and soilw fields
ifndsnowh = flag_snowh
IF (num_sw_levels_input .GE. 1) THEN
ifndsoilw = 1
ELSE
ifndsoilw = 0
END IF
! For bl_surface_physics = 1, we want to use close to a 10 cm value
! for the bottom level of the soil temps.
fix_bottom_level_for_temp : SELECT CASE ( model_config_rec%bl_surface_physics(grid%id) )
CASE (SLABSCHEME)
IF ( flag_st000010 .EQ. 1 ) THEN
DO j = jts , MIN(jde-1,jte)
DO i = its , MIN(ide-1,ite)
tmn(i,j) = st000010(i,j)
END DO
END DO
ELSE IF ( flag_soilt020 .EQ. 1 ) THEN
DO j = jts , MIN(jde-1,jte)
DO i = its , MIN(ide-1,ite)
tmn(i,j) = soilt020(i,j)
END DO
END DO
END IF
CASE (LSMSCHEME)
END SELECT fix_bottom_level_for_temp
! Generate the vegetation and soil category information from the fractional input
! data, or use the existing dominant category fields if they exist.
print *,'---> ---> forcing use of static soil and veg data, due to SI v1.3.1'
vegcat(its,jts) = 0
soilcat(its,jts) = 0
IF ( ( soilcat(its,jts) .LT. 0.5 ) .AND. ( vegcat(its,jts) .LT. 0.5 ) ) THEN
num_veg_cat = SIZE ( landusef , DIM=2 )
num_soil_top_cat = SIZE ( soilctop , DIM=2 )
num_soil_bot_cat = SIZE ( soilcbot , DIM=2 )
CALL process_percent_cat_new ( landmask , &
landusef , soilctop , soilcbot , &
isltyp , ivgtyp , &
num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
ids , ide , jds , jde , kds , kde , &
ims , ime , jms , jme , kms , kme , &
its , ite , jts , jte , kts , kte , &
model_config_rec%iswater(grid%id) )
! Make all the veg/soil parms the same so as not to confuse the developer.
DO j = jts , MIN(jde-1,jte)
DO i = its , MIN(ide-1,ite)
vegcat(i,j) = ivgtyp(i,j)
soilcat(i,j) = isltyp(i,j)
END DO
END DO
ELSE
! Do we have dominant soil and veg data from the input already?
IF ( soilcat(its,jts) .GT. 0.5 ) THEN
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
isltyp(i,j) = NINT( soilcat(i,j) )
END DO
END DO
END IF
IF ( vegcat(its,jts) .GT. 0.5 ) THEN
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
ivgtyp(i,j) = NINT( vegcat(i,j) )
END DO
END DO
END IF
END IF
! Land use assignment.
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
lu_index(i,j) = ivgtyp(i,j)
IF ( lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
landmask(i,j) = 1
xland(i,j) = 1
ELSE
landmask(i,j) = 0
xland(i,j) = 2
END IF
END DO
END DO
! Adjust the various soil temperature values depending on the difference in
! in elevation between the current model's elevation and the incoming data's
! orography.
IF ( flag_toposoil .EQ. 1 ) THEN
adjust_soil : SELECT CASE ( model_config_rec%bl_surface_physics(grid%id) )
CASE ( SLABSCHEME , LSMSCHEME )
CALL adjust_soil_temp_new ( tmn , model_config_rec%bl_surface_physics(grid%id) , &
tsk , ht , toposoil , landmask , flag_toposoil , &
st000010 , st010040 , st040100 , st100200 , st010200 , &
flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , &
flag_soilt160 , flag_soilt300 , &
ids , ide , jds , jde , kds , kde , &
ims , ime , jms , jme , kms , kme , &
its , ite , jts , jte , kts , kte )
END SELECT adjust_soil
END IF
! Fix tmn and tsk.
fix_tsk_tmn : SELECT CASE ( model_config_rec%bl_surface_physics(grid%id) )
CASE ( SLABSCHEME , LSMSCHEME )
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
tmn(i,j) = sst(i,j)
tsk(i,j) = sst(i,j)
ELSE IF ( landmask(i,j) .LT. 0.5 ) THEN
tmn(i,j) = tsk(i,j)
END IF
END DO
END DO
END SELECT fix_tsk_tmn
! Is the TSK reasonable?
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
IF ( tsk(i,j) .LT. 170 ) THEN
print *,'error in the TSK'
print *,'i,j=',i,j
print *,'landmask=',landmask(i,j)
print *,'tsk, sst, tmn=',tsk(i,j),sst(i,j),tmn(i,j)
stop
END IF
END DO
END DO
interpolate_soil_tmw : SELECT CASE ( model_config_rec%bl_surface_physics(grid%id) )
CASE ( SLABSCHEME , LSMSCHEME )
CALL process_soil_real ( tsk , tmn , &
landmask , sst , &
st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
zs , dzs , tslb , smois , sh2o , &
flag_sst , &
ids , ide , jds , jde , kds , kde , &
ims , ime , jms , jme , kms , kme , &
its , ite , jts , jte , kts , kte , &
model_config_rec%bl_surface_physics(grid%id) , &
model_config_rec%num_soil_layers , &
model_config_rec%real_data_init_type , &
num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
END SELECT interpolate_soil_tmw
! Let us make sure (again) that the landmask and the veg/soil categories match.
oops1=0
oops2=0
DO j = jts, MIN(jde-1,jte)
DO i = its, MIN(ide-1,ite)
IF ( ( ( landmask(i,j) .LT. 0.5 ) .AND. ( ivgtyp(i,j) .NE. config_flags%iswater .OR. isltyp(i,j) .NE. 14 ) ) .OR. &
( ( landmask(i,j) .GT. 0.5 ) .AND. ( ivgtyp(i,j) .EQ. config_flags%iswater .OR. isltyp(i,j) .EQ. 14 ) ) ) THEN
IF ( tslb(i,1,j) .GT. 1. ) THEN
oops1=oops1+1
ivgtyp(i,j) = 5
isltyp(i,j) = 8
landmask(i,j) = 1
xland(i,j) = 1
ELSE IF ( sst(i,j) .GT. 1. ) THEN
oops2=oops2+1
ivgtyp(i,j) = config_flags%iswater
isltyp(i,j) = 14
landmask(i,j) = 0
xland(i,j) = 2
ELSE
print *,'the landmask and soil/veg cats do not match'
print *,'i,j=',i,j
print *,'landmask=',landmask(i,j)
print *,'ivgtyp=',ivgtyp(i,j)
print *,'isltyp=',isltyp(i,j)
print *,'iswater=', config_flags%iswater
print *,'tslb=',tslb(i,:,j)
print *,'sst=',sst(i,j)
STOP 'mismatch_landmask_ivgtyp'
END IF
END IF
END DO
END DO
if (oops1.gt.0) then
print *,'points artificially set to land : ',oops1
endif
if(oops2.gt.0) then
print *,'points artificially set to water: ',oops2
endif
! From the full level data, we can get the half levels, reciprocals, and layer
! thicknesses. These are all defined at half level locations, so one less level.
! We allow the vertical coordinate to *accidently* come in upside down. We want
! the first full level to be the ground surface.
IF ( znw(1) .LT. znw(kde) ) THEN
DO k=1, kde/2
hold_znw = znw(k)
znw(k)=znw(kde+1-k)
znw(kde+1-k)=hold_znw
END DO
END IF
DO k=1, kde-1
dnw(k) = znw(k+1) - znw(k)
rdnw(k) = 1./dnw(k)
znu(k) = 0.5*(znw(k+1)+znw(k))
END DO
! Now the same sort of computations with the half eta levels, even ANOTHER
! level less than the one above.
DO k=2, kde-1
dn(k) = 0.5*(dnw(k)+dnw(k-1))
rdn(k) = 1./dn(k)
fnp(k) = .5* dnw(k )/dn(k)
fnm(k) = .5* dnw(k-1)/dn(k)
END DO
! Scads of vertical coefficients.
cof1 = (2.*dn(2)+dn(3))/(dn(2)+dn(3))*dnw(1)/dn(2)
cof2 = dn(2) /(dn(2)+dn(3))*dnw(1)/dn(3)
cf1 = fnp(2) + cof1
cf2 = fnm(2) - cof1 - cof2
cf3 = cof2
cfn = (.5*dnw(kde-1)+dn(kde-1))/dn(kde-1)
cfn1 = -.5*dnw(kde-1)/dn(kde-1)
! Inverse grid distances.
rdx = 1./dx
rdy = 1./dy
! Some of the many weird geopotential initializations that we'll see today: ph0 is total,
! and ph_2 is a perturbation from the base state geopotential. We set the base geopotential
! at the lowest level to terrain elevation * gravity.
DO j=jts,jte
DO i=its,ite
ph0(i,1,j) = ht(i,j) * g
ph_2(i,1,j) = 0.
END DO
END DO
! To define the base state, we call a USER MODIFIED routine to set the three
! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
! and A (temperature difference, from 1000 mb to 300 mb, K).
CALL const_module_initialize ( p00 , t00 , a )
! Base state potential temperature and inverse density (alpha = 1/rho) from
! the half eta levels and the base-profile surface pressure. Compute 1/rho
! from equation of state. The potential temperature is a perturbation from t0.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
! Base state pressure is a function of eta level and terrain, only, plus
! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*ht(i,j)/a/r_d ) **0.5 )
DO k = 1, kte-1
pb(i,k,j) = znu(k)*(p_surf - p_top) + p_top
t_init(i,k,j) = (t00 + A*LOG(pb(i,k,j)/p00))*(p00/pb(i,k,j))**(r_d/cp) - t0
alb(i,k,j) = (r_d/p1000mb)*(t_init(i,k,j)+t0)*(pb(i,k,j)/p1000mb)**cvpm
END DO
! Base state mu is defined as base state surface pressure minus p_top
mub(i,j) = p_surf - p_top
! Dry surface pressure is defined as the following (this mu is from the input file
! computed from the dry pressure). Here the dry pressure is just reconstituted.
pd_surf = mu0(i,j) + p_top
! Compute the perturbation mass and the full mass from the base surface pressure
! (a function of elevation) and the dry surface pressure.
mu_2(i,j) = pd_surf-p_surf
! Integrate base geopotential, starting at terrain elevation. This assures that
! the base state is in exact hydrostatic balance with respect to the model equations.
! This field is on full levels.
phb(i,1,j) = ht(i,j) * g
DO k = 2,kte
phb(i,k,j) = phb(i,k-1,j) - dnw(k-1)*mub(i,j)*alb(i,k-1,j)
END DO
END DO
END DO
! Fill in the outer rows and columns to allow us to be sloppy.
IF ( ite .EQ. ide ) THEN
i = ide
DO j = jts, MIN(jde-1,jte)
mub(i,j) = mub(i-1,j)
mu_2(i,j) = mu_2(i-1,j)
DO k = 1, kte-1
pb(i,k,j) = pb(i-1,k,j)
t_init(i,k,j) = t_init(i-1,k,j)
alb(i,k,j) = alb(i-1,k,j)
END DO
DO k = 1, kte
phb(i,k,j) = phb(i-1,k,j)
END DO
END DO
END IF
IF ( jte .EQ. jde ) THEN
j = jde
DO i = its, ite
mub(i,j) = mub(i,j-1)
mu_2(i,j) = mu_2(i,j-1)
DO k = 1, kte-1
pb(i,k,j) = pb(i,k,j-1)
t_init(i,k,j) = t_init(i,k,j-1)
alb(i,k,j) = alb(i,k,j-1)
END DO
DO k = 1, kte
phb(i,k,j) = phb(i,k,j-1)
END DO
END DO
END IF
DO j = jts, min(jde-1,jte)
DO i = its, min(ide-1,ite)
! Assign the potential temperature (perturbation from t0) and qv on all the mass
! point locations.
DO k = 1 , kde-1
t_2(i,k,j) = t_2(i,k,j) - t0
END DO
! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
! equation) down from the top to get the pressure perturbation. First get the pressure
! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
k = kte-1
qvf1 = 0.5*(moist_2(i,k,j,P_QV)+moist_2(i,k,j,P_QV))
qvf2 = 1./(1.+qvf1)
qvf1 = qvf1*qvf2
p(i,k,j) = - 0.5*(mu_2(i,j)+qvf1*mub(i,j))/rdnw(k)/qvf2
qvf = 1. + rvovrd*moist_2(i,k,j,P_QV)
alt(i,k,j) = (r_d/p1000mb)*(t_2(i,k,j)+t0)*qvf*(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm)
al(i,k,j) = alt(i,k,j) - alb(i,k,j)
! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
! inverse density fields (total and perturbation).
DO k=kte-2,1,-1
qvf1 = 0.5*(moist_2(i,k,j,P_QV)+moist_2(i,k+1,j,P_QV))
qvf2 = 1./(1.+qvf1)
qvf1 = qvf1*qvf2
p(i,k,j) = p(i,k+1,j) - (mu_2(i,j) + qvf1*mub(i,j))/qvf2/rdn(k+1)
qvf = 1. + rvovrd*moist_2(i,k,j,P_QV)
alt(i,k,j) = (r_d/p1000mb)*(t_2(i,k,j)+t0)*qvf* &
(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm)
al(i,k,j) = alt(i,k,j) - alb(i,k,j)
END DO
! This is the hydrostatic equation used in the model after the small timesteps. In
! the model, al (inverse density) is computed from the geopotential.
DO k = 2,kte
ph_2(i,k,j) = ph_2(i,k-1,j) - &
dnw(k-1) * ( (mub(i,j)+mu_2(i,j))*al(i,k-1,j) + mu_2(i,j)*alb(i,k-1,j) )
ph0(i,k,j) = ph_2(i,k,j) + phb(i,k,j)
END DO
END DO
END DO
#ifdef DM_PARALLEL
# include "HALO_EM_INIT.inc"
#endif
#define COPY_OUT
#include <em_scalar_derefs.inc>
RETURN
END SUBROUTINE init_domain_rk
!---------------------------------------------------------------------
SUBROUTINE const_module_initialize ( p00 , t00 , a ) 2
IMPLICIT NONE
REAL , PARAMETER :: sea_level_pressure_base = 100000.
REAL , PARAMETER :: sea_level_temperature_base = 290.
REAL , PARAMETER :: temp_diff_1000_to_300_mb = 50.
REAL , INTENT(OUT) :: p00 , t00 , a
p00 = sea_level_pressure_base
t00 = sea_level_temperature_base
a = temp_diff_1000_to_300_mb
END SUBROUTINE const_module_initialize
!-------------------------------------------------------------------
SUBROUTINE rebalance_driver ( grid ) 1,1
IMPLICIT NONE
TYPE (domain) :: grid
#ifdef DEREF_KLUDGE
INTEGER :: sm31 , em31 , sm32 , em32 , sm33 , em33
#endif
#ifdef DEREF_KLUDGE
sm31 = grid%sm31
em31 = grid%em31
sm32 = grid%sm32
em32 = grid%em32
sm33 = grid%sm33
em33 = grid%em33
#endif
CALL rebalance( grid, &
!
#include <em_actual_args.inc>
!
)
END SUBROUTINE rebalance_driver
!---------------------------------------------------------------------
SUBROUTINE rebalance ( grid , & 1,1
!
#include <em_dummy_args.inc>
!
)
IMPLICIT NONE
TYPE (domain) :: grid
#include <em_dummy_decl.inc>
TYPE (grid_config_rec_type) :: config_flags
REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold
REAL :: qvf , qvf1 , qvf2
REAL :: p00 , t00 , a
REAL , DIMENSION(:,:,:) , ALLOCATABLE :: t_init_int
! Local domain indices and counters.
INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
INTEGER :: &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte, &
i, j, k
#define COPY_IN
#include <em_scalar_derefs.inc>
#ifdef DM_PARALLEL
# include <em_data_calls.inc>
#endif
SELECT CASE ( model_data_order )
CASE ( DATA_ORDER_ZXY )
kds = grid%sd31 ; kde = grid%ed31 ;
ids = grid%sd32 ; ide = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
kms = grid%sm31 ; kme = grid%em31 ;
ims = grid%sm32 ; ime = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XYZ )
ids = grid%sd31 ; ide = grid%ed31 ;
jds = grid%sd32 ; jde = grid%ed32 ;
kds = grid%sd33 ; kde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
jms = grid%sm32 ; jme = grid%em32 ;
kms = grid%sm33 ; kme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
CASE ( DATA_ORDER_XZY )
ids = grid%sd31 ; ide = grid%ed31 ;
kds = grid%sd32 ; kde = grid%ed32 ;
jds = grid%sd33 ; jde = grid%ed33 ;
ims = grid%sm31 ; ime = grid%em31 ;
kms = grid%sm32 ; kme = grid%em32 ;
jms = grid%sm33 ; jme = grid%em33 ;
its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
END SELECT
ALLOCATE ( t_init_int(ims:ime,kms:kme,jms:jme) )
! Some of the many weird geopotential initializations that we'll see today: ph0 is total,
! and ph_2 is a perturbation from the base state geopotential. We set the base geopotential
! at the lowest level to terrain elevation * gravity.
DO j=jts,jte
DO i=its,ite
ph0(i,1,j) = ht_fine(i,j) * g
ph_2(i,1,j) = 0.
END DO
END DO
! To define the base state, we call a USER MODIFIED routine to set the three
! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
! and A (temperature difference, from 1000 mb to 300 mb, K).
CALL const_module_initialize ( p00 , t00 , a )
! Base state potential temperature and inverse density (alpha = 1/rho) from
! the half eta levels and the base-profile surface pressure. Compute 1/rho
! from equation of state. The potential temperature is a perturbation from t0.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
! Base state pressure is a function of eta level and terrain, only, plus
! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
! The fine grid terrain is ht_fine, the interpolated is ht.
p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*ht_fine(i,j)/a/r_d ) **0.5 )
p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*ht(i,j) /a/r_d ) **0.5 )
DO k = 1, kte-1
pb(i,k,j) = znu(k)*(p_surf - p_top) + p_top
pb_int = znu(k)*(p_surf_int - p_top) + p_top
t_init(i,k,j) = (t00 + A*LOG(pb(i,k,j)/p00))*(p00/pb(i,k,j))**(r_d/cp) - t0
t_init_int(i,k,j)= (t00 + A*LOG(pb_int /p00))*(p00/pb_int )**(r_d/cp) - t0
alb(i,k,j) = (r_d/p1000mb)*(t_init(i,k,j)+t0)*(pb(i,k,j)/p1000mb)**cvpm
END DO
! Base state mu is defined as base state surface pressure minus p_top
mub(i,j) = p_surf - p_top
! Dry surface pressure is defined as the following (this mu is from the input file
! computed from the dry pressure). Here the dry pressure is just reconstituted.
pd_surf = ( mub(i,j) + mu_2(i,j) ) + p_top
! Integrate base geopotential, starting at terrain elevation. This assures that
! the base state is in exact hydrostatic balance with respect to the model equations.
! This field is on full levels.
phb(i,1,j) = ht_fine(i,j) * g
DO k = 2,kte
phb(i,k,j) = phb(i,k-1,j) - dnw(k-1)*mub(i,j)*alb(i,k-1,j)
END DO
END DO
END DO
! Replace interpolated terrain with fine grid values.
DO j = jts, MIN(jte,jde-1)
DO i = its, MIN(ite,ide-1)
ht(i,j) = ht_fine(i,j)
END DO
END DO
! Perturbation fields.
DO j = jts, min(jde-1,jte)
DO i = its, min(ide-1,ite)
! The potential temperature is THETAnest = THETAinterp + ( TBARnest - TBARinterp)
DO k = 1 , kde-1
t_2(i,k,j) = t_2(i,k,j) + ( t_init(i,k,j) - t_init_int(i,k,j) )
END DO
! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
! equation) down from the top to get the pressure perturbation. First get the pressure
! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
k = kte-1
qvf1 = 0.5*(moist_2(i,k,j,P_QV)+moist_2(i,k,j,P_QV))
qvf2 = 1./(1.+qvf1)
qvf1 = qvf1*qvf2
p(i,k,j) = - 0.5*(mu_2(i,j)+qvf1*mub(i,j))/rdnw(k)/qvf2
qvf = 1. + rvovrd*moist_2(i,k,j,P_QV)
alt(i,k,j) = (r_d/p1000mb)*(t_2(i,k,j)+t0)*qvf*(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm)
al(i,k,j) = alt(i,k,j) - alb(i,k,j)
! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
! inverse density fields (total and perturbation).
DO k=kte-2,1,-1
qvf1 = 0.5*(moist_2(i,k,j,P_QV)+moist_2(i,k+1,j,P_QV))
qvf2 = 1./(1.+qvf1)
qvf1 = qvf1*qvf2
p(i,k,j) = p(i,k+1,j) - (mu_2(i,j) + qvf1*mub(i,j))/qvf2/rdn(k+1)
qvf = 1. + rvovrd*moist_2(i,k,j,P_QV)
alt(i,k,j) = (r_d/p1000mb)*(t_2(i,k,j)+t0)*qvf* &
(((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm)
al(i,k,j) = alt(i,k,j) - alb(i,k,j)
END DO
! This is the hydrostatic equation used in the model after the small timesteps. In
! the model, al (inverse density) is computed from the geopotential.
DO k = 2,kte
ph_2(i,k,j) = ph_2(i,k-1,j) - &
dnw(k-1) * ( (mub(i,j)+mu_2(i,j))*al(i,k-1,j) + mu_2(i,j)*alb(i,k-1,j) )
ph0(i,k,j) = ph_2(i,k,j) + phb(i,k,j)
END DO
END DO
END DO
DEALLOCATE ( t_init_int )
END SUBROUTINE rebalance
!---------------------------------------------------------------------
SUBROUTINE init_module_initialize,8
END SUBROUTINE init_module_initialize
!---------------------------------------------------------------------
END MODULE module_initialize