/wrfv2_fire/dyn_em/module_initialize_real.F
FORTRAN Legacy | 6042 lines | 3975 code | 1083 blank | 984 comment | 36 complexity | 0c1a13b5e2da8ab3fb4da9ef56917932 MD5 | raw file
Possible License(s): AGPL-1.0
- !REAL:MODEL_LAYER:INITIALIZATION
- #ifndef VERT_UNIT
- ! 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_real
- USE module_bc
- USE module_configure
- USE module_domain
- USE module_io_domain
- USE module_model_constants
- USE module_state_description
- USE module_timing
- USE module_soil_pre
- USE module_date_time
- USE module_llxy
- #ifdef DM_PARALLEL
- USE module_dm
- USE module_comm_dm, ONLY : &
- HALO_EM_INIT_1_sub &
- ,HALO_EM_INIT_2_sub &
- ,HALO_EM_INIT_3_sub &
- ,HALO_EM_INIT_4_sub &
- ,HALO_EM_INIT_5_sub &
- ,HALO_EM_VINTERP_UV_1_sub
- #endif
- REAL , SAVE :: p_top_save
- INTEGER :: internal_time_loop
- CONTAINS
- !-------------------------------------------------------------------
- SUBROUTINE init_domain ( grid )
- IMPLICIT NONE
- ! Input space and data. No gridded meteorological data has been stored, though.
- ! TYPE (domain), POINTER :: grid
- TYPE (domain) :: grid
- ! Local data.
- INTEGER :: idum1, idum2
- CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
- CALL init_domain_rk( grid &
- !
- #include "actual_new_args.inc"
- !
- )
- END SUBROUTINE init_domain
- !-------------------------------------------------------------------
- SUBROUTINE init_domain_rk ( grid &
- !
- #include "dummy_new_args.inc"
- !
- )
- USE module_optional_input
- IMPLICIT NONE
- ! Input space and data. No gridded meteorological data has been stored, though.
- ! TYPE (domain), POINTER :: grid
- TYPE (domain) :: grid
- #include "dummy_new_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 :: loop , num_seaice_changes
- INTEGER :: ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- its, ite, jts, jte, kts, kte, &
- ips, ipe, jps, jpe, kps, kpe, &
- i, j, k
- INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, &
- ipsx, ipex, jpsx, jpex, kpsx, kpex, &
- imsy, imey, jmsy, jmey, kmsy, kmey, &
- ipsy, ipey, jpsy, jpey, kpsy, kpey
- INTEGER :: ns
- ! Local data
- INTEGER :: error
- INTEGER :: im, num_3d_m, num_3d_s
- REAL :: p_surf, p_level
- REAL :: cof1, cof2
- REAL :: qvf , qvf1 , qvf2 , pd_surf
- REAL :: p00 , t00 , a , tiso
- REAL :: hold_znw , ptemp
- REAL :: vap_pres_mb , sat_vap_pres_mb
- LOGICAL :: were_bad
- LOGICAL :: stretch_grid, dry_sounding, debug
- INTEGER IICOUNT
- REAL :: p_top_requested , temp
- INTEGER :: num_metgrid_levels
- REAL , DIMENSION(max_eta) :: eta_levels
- REAL :: max_dz
- ! INTEGER , PARAMETER :: nl_max = 1000
- ! REAL , DIMENSION(nl_max) :: grid%dn
- integer::oops1,oops2
- REAL :: zap_close_levels
- INTEGER :: force_sfc_in_vinterp
- INTEGER :: interp_type , lagrange_order , extrap_type , t_extrap_type
- LOGICAL :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
- LOGICAL :: we_have_tavgsfc , we_have_tsk
- INTEGER :: lev500 , loop_count
- REAL :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu
- REAL :: pfu, pfd, phm
- LOGICAL , PARAMETER :: want_full_levels = .TRUE.
- LOGICAL , PARAMETER :: want_half_levels = .FALSE.
- CHARACTER (LEN=80) :: a_message
- REAL :: max_mf
-
- ! Excluded middle.
- LOGICAL :: any_valid_points
- INTEGER :: i_valid , j_valid
- !-- Carsel and Parrish [1988]
- REAL , DIMENSION(100) :: lqmi
- REAL :: t_start , t_end
- ! Dimension information stored in grid data structure.
- CALL cpu_time(t_start)
- CALL get_ijk_from_grid ( grid , &
- ids, ide, jds, jde, kds, kde, &
- ims, ime, jms, jme, kms, kme, &
- ips, ipe, jps, jpe, kps, kpe, &
- imsx, imex, jmsx, jmex, kmsx, kmex, &
- ipsx, ipex, jpsx, jpex, kpsx, kpex, &
- imsy, imey, jmsy, jmey, kmsy, kmey, &
- ipsy, ipey, jpsy, jpey, kpsy, kpey )
- its = ips ; ite = ipe ; jts = jps ; jte = jpe ; kts = kps ; kte = kpe
- CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
- ! Send out a quick message about the time steps based on the map scale factors.
- IF ( ( internal_time_loop .EQ. 1 ) .AND. ( grid%id .EQ. 1 ) .AND. &
- ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) ) THEN
- max_mf = grid%msft(its,jts)
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- max_mf = MAX ( max_mf , grid%msft(i,j) )
- END DO
- END DO
- #if ( defined(DM_PARALLEL) && ! defined(STUBMPI) )
- max_mf = wrf_dm_max_real ( max_mf )
- #endif
- WRITE ( a_message , FMT='(A,F5.2,A)' ) 'Max map factor in domain 1 = ',max_mf, &
- '. Scale the dt in the model accordingly.'
- CALL wrf_message ( a_message )
- END IF
- ! 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.
- grid%step_number = 0
- grid%itimestep=0
- ! Pull in the info in the namelist to compare it to the input data.
- grid%real_data_init_type = model_config_rec%real_data_init_type
- ! 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 , tiso )
- ! Save these constants to write out in model output file
- grid%t00 = t00
- grid%p00 = p00
- grid%tlp = a
- grid%tiso = tiso
- ! Are there any hold-ups to us bypassing the middle of the domain? These
- ! holdups would be situations where we need data in the middle of the domain.
- ! FOr example, if this si the first time period, we need the full domain
- ! processed for ICs. Also, if there is some sort of gridded FDDA turned on, or
- ! if the SST update is activated, then we can't just blow off the middle of the
- ! domain all willy-nilly. Other cases of these hold-ups? Sure - what if the
- ! user wants to smooth the CG topo, we need several rows and columns available.
- ! What if the lat/lon proj is used, then we need to run a spectral filter on
- ! the topo. Both are killers when trying to ignore data in the middle of the
- ! domain.
- ! If hold_ups = .F., then there are no hold-ups to excluding the middle
- ! domain processing. If hold_ups = .T., then there are hold-ups, and we
- ! must process the middle of the domain.
- hold_ups = ( internal_time_loop .EQ. 1 ) .OR. &
- ( config_flags%grid_fdda .NE. 0 ) .OR. &
- ( config_flags%sst_update .EQ. 1 ) .OR. &
- ( config_flags%all_ic_times ) .OR. &
- ( config_flags%smooth_cg_topo ) .OR. &
- ( config_flags%map_proj .EQ. PROJ_CASSINI )
- ! There are a few checks that we need to do when the input data comes in with the middle
- ! excluded by WPS.
- IF ( flag_excluded_middle .NE. 0 ) THEN
- ! If this time period of data from WPS has the middle excluded, it had better be OK for
- ! us to have a hole.
- IF ( hold_ups ) THEN
- WRITE ( a_message,* ) 'None of the following are allowed to be TRUE : '
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( internal_time_loop .EQ. 1 ) ', ( internal_time_loop .EQ. 1 )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( config_flags%grid_fdda .NE. 0 ) ', ( config_flags%grid_fdda .NE. 0 )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( config_flags%sst_update .EQ. 1 ) ', ( config_flags%sst_update .EQ. 1 )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( config_flags%all_ic_times ) ', ( config_flags%all_ic_times )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( config_flags%smooth_cg_topo ) ', ( config_flags%smooth_cg_topo )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) ' ( config_flags%map_proj .EQ. PROJ_CASSINI ) ', ( config_flags%map_proj .EQ. PROJ_CASSINI )
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) 'Problems, we cannot have excluded middle data from WPS'
- CALL wrf_error_fatal ( a_message )
- END IF
- ! Make sure that the excluded middle data from metgrid is "wide enough". We only have to check
- ! when the excluded middle was actually used in WPS.
- IF ( config_flags%spec_bdy_width .GT. flag_excluded_middle ) THEN
- WRITE ( a_message,* ) 'The WRF &bdy_control namelist.input spec_bdy_width = ', config_flags%spec_bdy_width
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) 'The WPS &metgrid namelist.wps process_only_bdy width = ',flag_excluded_middle
- CALL wrf_message ( a_message )
- WRITE ( a_message,* ) 'WPS process_only_bdy must be >= WRF spec_bdy_width'
- CALL wrf_error_fatal ( a_message )
- END IF
- END IF
- em_width = config_flags%spec_bdy_width
- ! We need to find if there are any valid non-excluded-middle points in this
- ! tile. If so, then we need to hang on to a valid i,j location.
- any_valid_points = .false.
- find_valid : DO j = jts,jte
- DO i = its,ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- any_valid_points = .true.
- i_valid = i
- j_valid = j
- EXIT find_valid
- END DO
- END DO find_valid
- ! Replace traditional seaice field with optional seaice (AFWA source)
- IF ( flag_icefrac .EQ. 1 ) THEN
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%xice(i,j) = grid%icefrac_gc(i,j)
- END DO
- END DO
- END IF
- ! Fix the snow (water equivalent depth, kg/m^2) and the snowh (physical snow
- ! depth, m) fields.
- IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%snow(i,j) = 0.
- grid%snowh(i,j) = 0.
- END DO
- END DO
- ELSE IF ( ( flag_snow .EQ. 0 ) .AND. ( flag_snowh .EQ. 1 ) ) THEN
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- ! ( m -> kg/m^2 ) & ( reduce to liquid, 5:1 ratio )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%snow(i,j) = grid%snowh(i,j) * 1000. / 5.
- END DO
- END DO
- ELSE IF ( ( flag_snow .EQ. 1 ) .AND. ( flag_snowh .EQ. 0 ) ) THEN
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- ! ( kg/m^2 -> m) & ( liquid to snow depth, 5:1 ratio )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%snowh(i,j) = grid%snow(i,j) / 1000. * 5.
- END DO
- END DO
- END IF
- ! For backward compatibility, we might need to assign the map factors from
- ! what they were, to what they are.
- IF ( ( config_flags%polar ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
- DO j=max(jds+1,jts),min(jde-1,jte)
- DO i=its,min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
- END DO
- END DO
- IF(jts == jds) THEN
- DO i=its,ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx(i,jts) = 0.
- grid%msfvx_inv(i,jts) = 0.
- END DO
- END IF
- IF(jte == jde) THEN
- DO i=its,ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx(i,jte) = 0.
- grid%msfvx_inv(i,jte) = 0.
- END DO
- END IF
- ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
- DO j=jts,jte
- DO i=its,min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
- END DO
- END DO
- ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
- DO j=jts,jte
- DO i=its,ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx(i,j) = grid%msfv(i,j)
- grid%msfvy(i,j) = grid%msfv(i,j)
- grid%msfux(i,j) = grid%msfu(i,j)
- grid%msfuy(i,j) = grid%msfu(i,j)
- grid%msftx(i,j) = grid%msft(i,j)
- grid%msfty(i,j) = grid%msft(i,j)
- ENDDO
- ENDDO
- DO j=jts,min(jde,jte)
- DO i=its,min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
- END DO
- END DO
- ELSE IF ( ( .NOT. config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .EQ. 1 ) ) THEN
- IF ( grid%msfvx(its,jts) .EQ. 0 ) THEN
- CALL wrf_error_fatal ( 'Maybe you do not have the new map factors, try re-running geogrid' )
- END IF
- DO j=jts,min(jde,jte)
- DO i=its,min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%msfvx_inv(i,j) = 1./grid%msfvx(i,j)
- END DO
- END DO
- ELSE IF ( ( config_flags%map_proj .EQ. PROJ_CASSINI ) .AND. ( flag_mf_xy .NE. 1 ) ) THEN
- CALL wrf_error_fatal ( 'Neither SI data nor older metgrid data can initialize a global domain' )
- ENDIF
- ! Check to see what available surface temperatures we have.
- IF ( flag_tavgsfc .EQ. 1 ) THEN
- we_have_tavgsfc = .TRUE.
- ELSE
- we_have_tavgsfc = .FALSE.
- END IF
- IF ( flag_tsk .EQ. 1 ) THEN
- we_have_tsk = .TRUE.
- ELSE
- we_have_tsk = .FALSE.
- END IF
-
- IF ( config_flags%use_tavg_for_tsk ) THEN
- IF ( we_have_tsk .OR. we_have_tavgsfc ) THEN
- ! we are OK
- ELSE
- CALL wrf_error_fatal ( 'We either need TSK or TAVGSFC, verify these fields are coming from WPS' )
- END IF
-
- ! Since we require a skin temperature in the model, we can use the average 2-m temperature if provided.
-
- IF ( we_have_tavgsfc ) THEN
- DO j=jts,min(jde,jte)
- DO i=its,min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%tsk(i,j) = grid%tavgsfc(i,j)
- END DO
- END DO
- END IF
- END IF
- ! Is there any vertical interpolation to do? The "old" data comes in on the correct
- ! vertical locations already.
- IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
- ! If this is data from the PINTERP program, it is emulating METGRID output.
- ! One of the caveats of this data is the way that the vertical structure is
- ! handled. We take the k=1 level and toss it (it is disposable), and we
- ! swap in the surface data. This is done for all of the 3d fields about
- ! which we show some interest: u, v, t, rh, ght, and p. For u, v, and rh,
- ! we assume no interesting vertical structure, and just assign the 1000 mb
- ! data. We directly use the 2-m temp for surface temp. We use the surface
- ! pressure field and the topography elevation for the lowest level of
- ! pressure and height, respectively.
- IF ( flag_pinterp .EQ. 1 ) THEN
- WRITE ( a_message , * ) 'Data from P_INTERP program, filling k=1 level with artificial surface fields.'
- CALL wrf_message ( a_message )
- DO j=jts,jte
- DO i=its,ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%u_gc(i,1,j) = grid%u_gc(i,2,j)
- grid%v_gc(i,1,j) = grid%v_gc(i,2,j)
- grid%rh_gc(i,1,j) = grid%rh_gc(i,2,j)
- grid%t_gc(i,1,j) = grid%t2(i,j)
- grid%ght_gc(i,1,j) = grid%ht(i,j)
- grid%p_gc(i,1,j) = grid%psfc(i,j)
- END DO
- END DO
- flag_psfc = 0
- END IF
- ! Variables that are named differently between SI and WPS.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%tsk(i,j) = grid%tsk_gc(i,j)
- grid%tmn(i,j) = grid%tmn_gc(i,j)
- grid%xlat(i,j) = grid%xlat_gc(i,j)
- grid%xlong(i,j) = grid%xlong_gc(i,j)
- grid%ht(i,j) = grid%ht_gc(i,j)
- END DO
- END DO
- ! A user could request that the most coarse grid has the
- ! topography along the outer boundary smoothed. This smoothing
- ! is similar to the coarse/nest interface. The outer rows and
- ! cols come from the existing large scale topo, and then the
- ! next several rows/cols are a linear ramp of the large scale
- ! model and the hi-res topo from WPS. We only do this for the
- ! coarse grid since we are going to make the interface consistent
- ! in the model betwixt the CG and FG domains.
- IF ( ( config_flags%smooth_cg_topo ) .AND. &
- ( grid%id .EQ. 1 ) .AND. &
- ( flag_soilhgt .EQ. 1) ) THEN
- CALL blend_terrain ( grid%toposoil , grid%ht , &
- ids , ide , jds , jde , 1 , 1 , &
- ims , ime , jms , jme , 1 , 1 , &
- ips , ipe , jps , jpe , 1 , 1 )
- END IF
- ! Filter the input topography if this is a polar projection.
- IF ( ( config_flags%polar ) .AND. ( grid%fft_filter_lat .GT. 90 ) ) THEN
- CALL wrf_error_fatal ( 'If the polar boundary condition is used, then fft_filter_lat must be set in namelist.input' )
- END IF
- IF ( config_flags%map_proj .EQ. PROJ_CASSINI ) THEN
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- ! We stick the topo and map fac in an unused 3d array. The map scale
- ! factor and computational latitude are passed along for the ride
- ! (part of the transpose process - we only do 3d arrays) to determine
- ! "how many" values are used to compute the mean. We want a number
- ! that is consistent with the original grid resolution.
- DO j = jts, MIN(jte,jde-1)
- DO k = kts, kte
- DO i = its, MIN(ite,ide-1)
- grid%t_init(i,k,j) = 1.
- END DO
- END DO
- DO i = its, MIN(ite,ide-1)
- grid%t_init(i,1,j) = grid%ht(i,j)
- grid%t_init(i,2,j) = grid%msftx(i,j)
- grid%t_init(i,3,j) = grid%clat(i,j)
- END DO
- END DO
- # include "XPOSE_POLAR_FILTER_TOPO_z2x.inc"
- ! Retrieve the 2d arrays for topo, map factors, and the
- ! computational latitude.
- DO j = jpsx, MIN(jpex,jde-1)
- DO i = ipsx, MIN(ipex,ide-1)
- grid%ht_xxx(i,j) = grid%t_xxx(i,1,j)
- grid%mf_xxx(i,j) = grid%t_xxx(i,2,j)
- grid%clat_xxx(i,j) = grid%t_xxx(i,3,j)
- END DO
- END DO
- ! Get a mean topo field that is consistent with the grid
- ! distance on each computational latitude loop.
- CALL filter_topo ( grid%ht_xxx , grid%clat_xxx , grid%mf_xxx , &
- grid%fft_filter_lat , &
- ids, ide, jds, jde, 1 , 1 , &
- imsx, imex, jmsx, jmex, 1, 1, &
- ipsx, ipex, jpsx, jpex, 1, 1 )
- ! Stick the filtered topo back into the dummy 3d array to
- ! transpose it back to "all z on a patch".
- DO j = jpsx, MIN(jpex,jde-1)
- DO i = ipsx, MIN(ipex,ide-1)
- grid%t_xxx(i,1,j) = grid%ht_xxx(i,j)
- END DO
- END DO
- # include "XPOSE_POLAR_FILTER_TOPO_x2z.inc"
- ! Get the un-transposed topo data.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- grid%ht(i,j) = grid%t_init(i,1,j)
- END DO
- END DO
- #else
- CALL filter_topo ( grid%ht , grid%clat , grid%msftx , &
- grid%fft_filter_lat , &
- ids, ide, jds, jde, 1,1, &
- ims, ime, jms, jme, 1,1, &
- its, ite, jts, jte, 1,1 )
- #endif
- END IF
- ! If we have any input low-res surface pressure, we store it.
- IF ( flag_psfc .EQ. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%psfc_gc(i,j) = grid%psfc(i,j)
- grid%p_gc(i,1,j) = grid%psfc(i,j)
- END DO
- END DO
- END IF
- ! If we have the low-resolution surface elevation, stick that in the
- ! "input" locations of the 3d height. We still have the "hi-res" topo
- ! stuck in the grid%ht array. The grid%landmask if test is required as some sources
- ! have ZERO elevation over water (thank you very much).
- IF ( flag_soilhgt .EQ. 1) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- ! IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%ght_gc(i,1,j) = grid%toposoil(i,j)
- grid%ht_gc(i,j)= grid%toposoil(i,j)
- ! END IF
- END DO
- END DO
- END IF
- ! The number of vertical levels in the input data. There is no staggering for
- ! different variables.
- num_metgrid_levels = grid%num_metgrid_levels
- ! For UM data, swap incoming extra (theta-based) pressure with the standardly
- ! named (rho-based) pressure.
- IF ( flag_ptheta .EQ. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO k = 1 , num_metgrid_levels
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ptemp = grid%p_gc(i,k,j)
- grid%p_gc(i,k,j) = grid%prho_gc(i,k,j)
- grid%prho_gc(i,k,j) = ptemp
- END DO
- END DO
- END DO
- ! For UM data, the "surface" and the "first hybrid" level for the theta-level data fields are the same.
- ! Average the surface (k=1) and the second hybrid level (k=num_metgrid_levels-1) to get the first hybrid
- ! layer. We only do this for the theta-level data: pressure, temperature, specific humidity, and
- ! geopotential height (i.e. we do not modify u, v, or the rho-based pressure).
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid% p_gc(i,num_metgrid_levels,j) = ( grid% p_gc(i,1,j) + grid% p_gc(i,num_metgrid_levels-1,j) ) * 0.5
- grid% t_gc(i,num_metgrid_levels,j) = ( grid% t_gc(i,1,j) + grid% t_gc(i,num_metgrid_levels-1,j) ) * 0.5
- grid% sh_gc(i,num_metgrid_levels,j) = ( grid% sh_gc(i,1,j) + grid% sh_gc(i,num_metgrid_levels-1,j) ) * 0.5
- grid%ght_gc(i,num_metgrid_levels,j) = ( grid%ght_gc(i,1,j) + grid%ght_gc(i,num_metgrid_levels-1,j) ) * 0.5
- END DO
- END DO
- END IF
- IF ( any_valid_points ) THEN
- ! Check for and semi-fix missing surface fields.
- IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
- k = 2
- ELSE
- k = num_metgrid_levels
- END IF
- IF ( grid%t_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%t_gc(i,1,j) = grid%t_gc(i,k,j)
- END DO
- END DO
- config_flags%use_surface = .FALSE.
- grid%use_surface = .FALSE.
- WRITE ( a_message , * ) 'Missing surface temp, replaced with closest level, use_surface set to false.'
- CALL wrf_message ( a_message )
- END IF
- IF ( grid%rh_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%rh_gc(i,1,j) = grid%rh_gc(i,k,j)
- END DO
- END DO
- config_flags%use_surface = .FALSE.
- grid%use_surface = .FALSE.
- WRITE ( a_message , * ) 'Missing surface RH, replaced with closest level, use_surface set to false.'
- CALL wrf_message ( a_message )
- END IF
- IF ( grid%u_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, ite
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%u_gc(i,1,j) = grid%u_gc(i,k,j)
- END DO
- END DO
- config_flags%use_surface = .FALSE.
- grid%use_surface = .FALSE.
- WRITE ( a_message , * ) 'Missing surface u wind, replaced with closest level, use_surface set to false.'
- CALL wrf_message ( a_message )
- END IF
- IF ( grid%v_gc(i_valid,1,j_valid) .EQ. -1.E30 ) THEN
- DO j = jts, jte
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%v_gc(i,1,j) = grid%v_gc(i,k,j)
- END DO
- END DO
- config_flags%use_surface = .FALSE.
- grid%use_surface = .FALSE.
- WRITE ( a_message , * ) 'Missing surface v wind, replaced with closest level, use_surface set to false.'
- CALL wrf_message ( a_message )
- END IF
- ! Compute the mixing ratio from the input relative humidity.
- IF ( ( flag_qv .NE. 1 ) .AND. ( flag_sh .NE. 1 ) ) THEN
- IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
- k = 2
- ELSE
- k = num_metgrid_levels
- END IF
- IF ( config_flags%rh2qv_method .eq. 1 ) THEN
- CALL rh_to_mxrat1(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc , &
- config_flags%rh2qv_wrt_liquid , &
- config_flags%qv_max_p_safe , &
- config_flags%qv_max_flag , config_flags%qv_max_value , &
- config_flags%qv_min_p_safe , &
- config_flags%qv_min_flag , config_flags%qv_min_value , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ELSE IF ( config_flags%rh2qv_method .eq. 2 ) THEN
- CALL rh_to_mxrat2(grid%rh_gc, grid%t_gc, grid%p_gc, grid%qv_gc , &
- config_flags%rh2qv_wrt_liquid , &
- config_flags%qv_max_p_safe , &
- config_flags%qv_max_flag , config_flags%qv_max_value , &
- config_flags%qv_min_p_safe , &
- config_flags%qv_min_flag , config_flags%qv_min_value , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- END IF
- ELSE IF ( flag_sh .EQ. 1 ) THEN
- IF ( grid%p_gc(i_valid,num_metgrid_levels,j_valid) .LT. grid%p_gc(i_valid,2,j_valid) ) THEN
- k = 2
- ELSE
- k = num_metgrid_levels
- END IF
- IF ( grid%sh_gc(i_valid,kts,j_valid) .LT. 1.e-6 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%sh_gc(i,1,j) = grid%sh_gc(i,k,j)
- END DO
- END DO
- END IF
- DO j = jts, MIN(jte,jde-1)
- DO k = 1 , num_metgrid_levels
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%qv_gc(i,k,j) = grid%sh_gc(i,k,j) /( 1. - grid%sh_gc(i,k,j) )
- sat_vap_pres_mb = 0.6112*10.*EXP(17.67*(grid%t_gc(i,k,j)-273.15)/(grid%t_gc(i,k,j)-29.65))
- vap_pres_mb = grid%qv_gc(i,k,j) * grid%p_gc(i,k,j)/100. / (grid%qv_gc(i,k,j) + 0.622 )
- IF ( sat_vap_pres_mb .GT. 0 ) THEN
- grid%rh_gc(i,k,j) = ( vap_pres_mb / sat_vap_pres_mb ) * 100.
- ELSE
- grid%rh_gc(i,k,j) = 0.
- END IF
- END DO
- END DO
- END DO
- END IF
- ! Some data sets do not provide a 3d geopotential height field.
- IF ( grid%ght_gc(i_valid,grid%num_metgrid_levels/2,j_valid) .LT. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO k = kts+1 , grid%num_metgrid_levels
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%ght_gc(i,k,j) = grid%ght_gc(i,k-1,j) - &
- R_d / g * 0.5 * ( grid%t_gc(i,k ,j) * ( 1 + 0.608 * grid%qv_gc(i,k ,j) ) + &
- grid%t_gc(i,k-1,j) * ( 1 + 0.608 * grid%qv_gc(i,k-1,j) ) ) * &
- LOG ( grid%p_gc(i,k,j) / grid%p_gc(i,k-1,j) )
- END DO
- END DO
- END DO
- END IF
- ! If the pressure levels in the middle of the atmosphere are upside down, then
- ! this is hybrid data. Computing the new surface pressure should use sfcprs2.
- IF ( grid%p_gc(i_valid,num_metgrid_levels/2,j_valid) .LT. grid%p_gc(i_valid,num_metgrid_levels/2+1,j_valid) ) THEN
- config_flags%sfcp_to_sfcp = .TRUE.
- END IF
- END IF
- ! Assign surface fields with original input values. If this is hybrid data,
- ! the values are not exactly representative. However - this is only for
- ! plotting purposes and such at the 0h of the forecast, so we are not all that
- ! worried.
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%u10(i,j)=grid%u_gc(i,1,j)
- END DO
- END DO
- DO j = jts, min(jde,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%v10(i,j)=grid%v_gc(i,1,j)
- END DO
- END DO
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%t2(i,j)=grid%t_gc(i,1,j)
- END DO
- END DO
- IF ( flag_qv .EQ. 1 ) THEN
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%q2(i,j)=grid%qv_gc(i,1,j)
- END DO
- END DO
- END IF
- ! The requested ptop for real data cases.
- p_top_requested = grid%p_top_requested
- ! Compute the top pressure, grid%p_top. For isobaric data, this is just the
- ! top level. For the generalized vertical coordinate data, we find the
- ! max pressure on the top level. We have to be careful of two things:
- ! 1) the value has to be communicated, 2) the value can not increase
- ! at subsequent times from the initial value.
- IF ( internal_time_loop .EQ. 1 ) THEN
- CALL find_p_top ( grid%p_gc , grid%p_top , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- grid%p_top = wrf_dm_max_real ( grid%p_top )
- #endif
- ! Compare the requested grid%p_top with the value available from the input data.
- IF ( p_top_requested .LT. grid%p_top ) THEN
- print *,'p_top_requested = ',p_top_requested
- print *,'allowable grid%p_top in data = ',grid%p_top
- CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' )
- END IF
- ! The grid%p_top valus is the max of what is available from the data and the
- ! requested value. We have already compared <, so grid%p_top is directly set to
- ! the value in the namelist.
- grid%p_top = p_top_requested
- ! For subsequent times, we have to remember what the grid%p_top for the first
- ! time was. Why? If we have a generalized vert coordinate, the grid%p_top value
- ! could fluctuate.
- p_top_save = grid%p_top
- ELSE
- CALL find_p_top ( grid%p_gc , grid%p_top , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
- grid%p_top = wrf_dm_max_real ( grid%p_top )
- #endif
- IF ( grid%p_top .GT. p_top_save ) THEN
- print *,'grid%p_top from last time period = ',p_top_save
- print *,'grid%p_top from this time period = ',grid%p_top
- CALL wrf_error_fatal ( 'grid%p_top > previous value' )
- END IF
- grid%p_top = p_top_save
- ENDIF
- ! Get the monthly values interpolated to the current date for the traditional monthly
- ! fields of green-ness fraction and background albedo.
- CALL monthly_interp_to_date ( grid%greenfrac , current_date , grid%vegfra , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- CALL monthly_interp_to_date ( grid%albedo12m , current_date , grid%albbck , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Get the min/max of each i,j for the monthly green-ness fraction.
- CALL monthly_min_max ( grid%greenfrac , grid%shdmin , grid%shdmax , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! The model expects the green-ness values in percent, not fraction.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
- grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
- grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
- END DO
- END DO
- ! The model expects the albedo fields as a fraction, not a percent. Set the
- ! water values to 8%.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%albbck(i,j) = grid%albbck(i,j) / 100.
- grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
- IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
- grid%albbck(i,j) = 0.08
- grid%snoalb(i,j) = 0.08
- END IF
- END DO
- END DO
- ! Two ways to get the surface pressure. 1) If we have the low-res input surface
- ! pressure and the low-res topography, then we can do a simple hydrostatic
- ! relation. 2) Otherwise we compute the surface pressure from the sea-level
- ! pressure.
- ! Note that on output, grid%psfc is now hi-res. The low-res surface pressure and
- ! elevation are grid%psfc_gc and grid%ht_gc (same as grid%ght_gc(k=1)).
- IF ( ( flag_psfc .EQ. 1 ) .AND. &
- ( flag_soilhgt .EQ. 1 ) .AND. &
- ( flag_slp .EQ. 1 ) .AND. &
- ( .NOT. config_flags%sfcp_to_sfcp ) ) THEN
- WRITE(a_message,FMT='(A)') 'Using sfcprs3 to compute psfc'
- CALL wrf_message ( a_message )
- CALL sfcprs3(grid%ght_gc, grid%p_gc, grid%ht, &
- grid%pslv_gc, grid%psfc, &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ELSE IF ( ( flag_psfc .EQ. 1 ) .AND. &
- ( flag_soilhgt .EQ. 1 ) .AND. &
- ( config_flags%sfcp_to_sfcp ) ) THEN
- WRITE(a_message,FMT='(A)') 'Using sfcprs2 to compute psfc'
- CALL wrf_message ( a_message )
- CALL sfcprs2(grid%t_gc, grid%qv_gc, grid%ght_gc, grid%psfc_gc, grid%ht, &
- grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ELSE IF ( flag_slp .EQ. 1 ) THEN
- WRITE(a_message,FMT='(A)') 'Using sfcprs to compute psfc'
- CALL wrf_message ( a_message )
- CALL sfcprs (grid%t_gc, grid%qv_gc, grid%ght_gc, grid%pslv_gc, grid%ht, &
- grid%tavgsfc, grid%p_gc, grid%psfc, we_have_tavgsfc, &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ELSE
- WRITE(a_message,FMT='(3(A,I2),A,L1)') 'ERROR in psfc: flag_psfc = ',flag_psfc, &
- ', flag_soilhgt = ',flag_soilhgt , &
- ', flag_slp = ',flag_slp , &
- ', sfcp_to_sfcp = ',config_flags%sfcp_to_sfcp
- CALL wrf_message ( a_message )
- CALL wrf_error_fatal ( 'not enough info for a p sfc computation' )
- END IF
- ! If we have no input surface pressure, we'd better stick something in there.
- IF ( flag_psfc .NE. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%psfc_gc(i,j) = grid%psfc(i,j)
- grid%p_gc(i,1,j) = grid%psfc(i,j)
- END DO
- END DO
- END IF
- ! Integrate the mixing ratio to get the vapor pressure.
- CALL integ_moist ( grid%qv_gc , grid%p_gc , grid%pd_gc , grid%t_gc , grid%ght_gc , grid%intq_gc , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ! If this is UM data, the same moisture removed from the "theta" level pressure data can
- ! be removed from the "rho" level pressures. This is an approximation. We'll revisit to
- ! see if this is a bad idea.
- IF ( flag_ptheta .EQ. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO k = num_metgrid_levels-1 , 1 , -1
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ptemp = ((grid%p_gc(i,k,j) - grid%pd_gc(i,k,j)) + (grid%p_gc(i,k+1,j) - grid%pd_gc(i,k+1,j)))/2
- grid%pdrho_gc(i,k,j) = grid%prho_gc(i,k,j) - ptemp
- END DO
- END DO
- END DO
- END IF
- ! Compute the difference between the dry, total surface pressure (input) and the
- ! dry top pressure (constant).
- CALL p_dts ( grid%mu0 , grid%intq_gc , grid%psfc , grid%p_top , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- ! Compute the dry, hydrostatic surface pressure.
- CALL p_dhs ( grid%pdhs , grid%ht , p00 , t00 , a , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute the eta levels if not defined already.
- IF ( grid%znw(1) .NE. 1.0 ) THEN
- eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
- max_dz = model_config_rec%max_dz
- CALL compute_eta ( grid%znw , &
- eta_levels , max_eta , max_dz , &
- grid%p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- IF ( config_flags%interp_theta ) THEN
- ! The input field is temperature, we want potential temp.
- CALL t_to_theta ( grid%t_gc , grid%p_gc , p00 , &
- ids , ide , jds , jde , 1 , num_metgrid_levels , &
- ims , ime , jms , jme , 1 , num_metgrid_levels , &
- its , ite , jts , jte , 1 , num_metgrid_levels )
- END IF
- IF ( flag_slp .EQ. 1 ) THEN
- ! On the eta surfaces, compute the dry pressure = mu eta, stored in
- ! grid%pb, since it is a pressure, and we don't need another kms:kme 3d
- ! array floating around. The grid%pb array is re-computed as the base pressure
- ! later after the vertical interpolations are complete.
- CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_full_levels , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! All of the vertical interpolations are done in dry-pressure space. The
- ! input data has had the moisture removed (grid%pd_gc). The target levels (grid%pb)
- ! had the vapor pressure removed from the surface pressure, then they were
- ! scaled by the eta levels.
- interp_type = 2
- lagrange_order = grid%lagrange_order
- lowest_lev_from_sfc = .FALSE.
- use_levels_below_ground = .TRUE.
- use_surface = .TRUE.
- zap_close_levels = grid%zap_close_levels
- force_sfc_in_vinterp = 0
- t_extrap_type = grid%t_extrap_type
- extrap_type = 1
- ! For the height field, the lowest level pressure is the slp (approximately "dry"). The
- ! lowest level of the input height field (to be associated with slp) then is an array
- ! of zeros.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- grid%psfc_gc(i,j) = grid%pd_gc(i,1,j)
- grid%pd_gc(i,1,j) = grid%pslv_gc(i,j) - ( grid%p_gc(i,1,j) - grid%pd_gc(i,1,j) )
- grid%ht_gc(i,j) = grid%ght_gc(i,1,j)
- grid%ght_gc(i,1,j) = 0.
- END DO
- END DO
- CALL vert_interp ( grid%ght_gc , grid%pd_gc , grid%ph0 , grid%pb , &
- num_metgrid_levels , 'Z' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Put things back to normal.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%pd_gc(i,1,j) = grid%psfc_gc(i,j)
- grid%ght_gc(i,1,j) = grid%ht_gc(i,j)
- END DO
- END DO
- END IF
- ! Now the rest of the variables on half-levels to inteprolate.
- CALL p_dry ( grid%mu0 , grid%znw , grid%p_top , grid%pb , want_half_levels , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- interp_type = grid%interp_type
- lagrange_order = grid%lagrange_order
- lowest_lev_from_sfc = grid%lowest_lev_from_sfc
- use_levels_below_ground = grid%use_levels_below_ground
- use_surface = grid%use_surface
- zap_close_levels = grid%zap_close_levels
- force_sfc_in_vinterp = grid%force_sfc_in_vinterp
- t_extrap_type = grid%t_extrap_type
- extrap_type = grid%extrap_type
- ! Interpolate RH, diagnose Qv later when have temp and pressure. Temporarily
- ! store this in the u_1 space, for later diagnosis into Qv and stored into moist.
- CALL vert_interp ( grid%rh_gc , grid%pd_gc , grid%u_1 , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Depending on the setting of interp_theta = T/F, t_gc is is either theta Xor
- ! temperature, and that means that the t_2 field is also the associated field.
- ! It is better to interpolate temperature and potential temperature in LOG(p),
- ! regardless of requested default.
- interp_type = 2
- CALL vert_interp ( grid%t_gc , grid%pd_gc , grid%t_2 , grid%pb , &
- num_metgrid_levels , 'T' , &
- interp_type , lagrange_order , t_extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- interp_type = grid%interp_type
-
- ! It is better to interpolate pressure in p regardless default options
- interp_type = 1
- CALL vert_interp ( grid%p_gc , grid%pd_gc , grid%p , grid%pb , &
- num_metgrid_levels , 'T' , &
- interp_type , lagrange_order , t_extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- interp_type = grid%interp_type
- ! Do not have full pressure on eta levels, get a first guess at Qv by using
- ! dry pressure. The use of u_1 (rh) and v_1 (temperature) is temporary.
- ! We fix the approximation to Qv after the total pressure is available on
- ! eta surfaces.
- grid%v_1 = grid%t_2
- IF ( config_flags%interp_theta ) THEN
- CALL theta_to_t ( grid%v_1 , grid%p , p00 , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- IF ( config_flags%rh2qv_method .eq. 1 ) THEN
- CALL rh_to_mxrat1(grid%u_1, grid%v_1, grid%p , moist(:,:,:,P_QV) , &
- config_flags%rh2qv_wrt_liquid , &
- config_flags%qv_max_p_safe , &
- config_flags%qv_max_flag , config_flags%qv_max_value , &
- config_flags%qv_min_p_safe , &
- config_flags%qv_min_flag , config_flags%qv_min_value , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte-1 )
- ELSE IF ( config_flags%rh2qv_method .eq. 2 ) THEN
- CALL rh_to_mxrat2(grid%u_1, grid%v_1, grid%p , moist(:,:,:,P_QV) , &
- config_flags%rh2qv_wrt_liquid , &
- config_flags%qv_max_p_safe , &
- config_flags%qv_max_flag , config_flags%qv_max_value , &
- config_flags%qv_min_p_safe , &
- config_flags%qv_min_flag , config_flags%qv_min_value , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte-1 )
- END IF
- IF ( .NOT. config_flags%interp_theta ) THEN
- CALL t_to_theta ( grid%t_2 , grid%p , p00 , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- num_3d_m = num_moist
- num_3d_s = num_scalar
- IF ( flag_qr .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QR ) THEN
- CALL vert_interp ( grid%qr_gc , grid%pd_gc , moist(:,:,:,P_QR) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qc .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QC ) THEN
- CALL vert_interp ( grid%qc_gc , grid%pd_gc , moist(:,:,:,P_QC) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qi .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QI ) THEN
- CALL vert_interp ( grid%qi_gc , grid%pd_gc , moist(:,:,:,P_QI) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qs .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QS ) THEN
- CALL vert_interp ( grid%qs_gc , grid%pd_gc , moist(:,:,:,P_QS) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qg .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QG ) THEN
- CALL vert_interp ( grid%qg_gc , grid%pd_gc , moist(:,:,:,P_QG) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qh .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_m
- IF ( im .EQ. P_QH ) THEN
- CALL vert_interp ( grid%qh_gc , grid%pd_gc , moist(:,:,:,P_QH) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- IF ( flag_qni .EQ. 1 ) THEN
- DO im = PARAM_FIRST_SCALAR, num_3d_s
- IF ( im .EQ. P_QNI ) THEN
- CALL vert_interp ( grid%qni_gc , grid%pd_gc , scalar(:,:,:,P_QNI) , grid%pb , &
- num_metgrid_levels , 'Q' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF
- END DO
- END IF
- ! If this is UM data, put the dry rho-based pressure back into the dry pressure array.
- ! Since the dry pressure is no longer needed, no biggy.
- IF ( flag_ptheta .EQ. 1 ) THEN
- DO j = jts, MIN(jte,jde-1)
- DO k = 1 , num_metgrid_levels
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%pd_gc(i,k,j) = grid%prho_gc(i,k,j)
- END DO
- END DO
- END DO
- END IF
- #ifdef DM_PARALLEL
- ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
- ! For the U and V vertical interpolation, we need the pressure defined
- ! at both the locations for the horizontal momentum, which we get by
- ! averaging two pressure values (i and i-1 for U, j and j-1 for V). The
- ! pressure field on input (grid%pd_gc) and the pressure of the new coordinate
- ! (grid%pb) are both communicated with an 8 stencil.
- # include "HALO_EM_VINTERP_UV_1.inc"
- #endif
- CALL vert_interp ( grid%u_gc , grid%pd_gc , grid%u_2 , grid%pb , &
- num_metgrid_levels , 'U' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- CALL vert_interp ( grid%v_gc , grid%pd_gc , grid%v_2 , grid%pb , &
- num_metgrid_levels , 'V' , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
- ! Set the temperature of the inland lakes to tavgsfc if the temperature is available
- ! and islake is > num_veg_cat
- num_veg_cat = SIZE ( grid%landusef , DIM=2 )
- CALL nl_get_iswater ( grid%id , grid%iswater )
- CALL nl_get_islake ( grid%id , grid%islake )
- IF ( grid%islake < 0 ) THEN
- CALL wrf_debug ( 0 , 'Old data, no inland lake information')
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( ( grid%landusef(i,grid%iswater,j) >= 0.5 ) .OR. ( grid%lu_index(i,j) == grid%iswater ) ) .AND. &
- ( ( grid%sst(i,j) .LT. 150 ) .OR. ( grid%sst(i,j) .GT. 400 ) ) ) THEN
- IF ( we_have_tavgsfc ) THEN
- grid%sst(i,j) = grid%tavgsfc(i,j)
- END IF
- IF ( ( grid%sst(i,j) .LT. 150 ) .OR. ( grid%sst(i,j) .GT. 400 ) ) THEN
- grid%sst(i,j) = grid%tsk(i,j)
- END IF
- IF ( ( grid%sst(i,j) .LT. 150 ) .OR. ( grid%sst(i,j) .GT. 400 ) ) THEN
- grid%sst(i,j) = grid%t2(i,j)
- END IF
- END IF
- END DO
- END DO
- ELSE
- IF ( we_have_tavgsfc ) THEN
- CALL wrf_debug ( 0 , 'Using inland lakes with average surface temperature')
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( grid%landusef(i,grid%islake,j) >= 0.5 ) .OR. ( grid%lu_index(i,j) == grid%islake ) ) THEN
- grid%sst(i,j) = grid%tavgsfc(i,j)
- grid%tsk(i,j) = grid%tavgsfc(i,j)
- END IF
- IF ( ( grid%sst(i,j) .LT. 150 ) .OR. ( grid%sst(i,j) .GT. 400 ) ) THEN
- grid%sst(i,j) = grid%t2(i,j)
- END IF
- END DO
- END DO
- ELSE ! We don't have tavgsfc
- CALL wrf_debug ( 0 , 'No average surface temperature for use with inland lakes')
- END IF
- DO j=jts,MIN(jde-1,jte)
- DO i=its,MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%landusef(i,grid%iswater,j) = grid%landusef(i,grid%iswater,j) + &
- grid%landusef(i,grid%islake,j)
- grid%landusef(i,grid%islake,j) = 0.
- END DO
- END DO
- END IF
- ! Save the grid%tsk field for later use in the sea ice surface temperature
- ! for the Noah LSM scheme.
- DO j = jts, MIN(jte,jde-1)
- DO i = its, MIN(ite,ide-1)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%tsk_save(i,j) = grid%tsk(i,j)
- END DO
- END DO
- ! Protect against bad grid%tsk values over water by supplying grid%sst (if it is
- ! available, and if the grid%sst is reasonable).
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
- ( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
- grid%tsk(i,j) = grid%sst(i,j)
- ENDIF
- END DO
- END DO
- ! 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 ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( grid%snow(i,j) .GE. 10. ) then
- grid%snowc(i,j) = 1.
- ELSE
- grid%snowc(i,j) = 0.0
- END IF
- END DO
- END DO
- ! Set flag integers for presence of snowh and soilw fields
- grid%ifndsnowh = flag_snowh
- IF (num_sw_levels_input .GE. 1) THEN
- grid%ifndsoilw = 1
- ELSE
- grid%ifndsoilw = 0
- END IF
- ! We require input data for the various LSM schemes.
- enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
- CASE ( LSMSCHEME, NOAHMPSCHEME )
- IF ( num_st_levels_input .LT. 2 ) THEN
- CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
- END IF
- CASE (RUCLSMSCHEME)
- IF ( num_st_levels_input .LT. 2 ) THEN
- CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
- END IF
- CASE (PXLSMSCHEME)
- IF ( num_st_levels_input .LT. 2 ) THEN
- CALL wrf_error_fatal ( 'Not enough soil temperature data for P-X LSM scheme.')
- END IF
- !---------- fds (06/2010) ---------------------------------
- CASE (SSIBSCHEME)
- IF ( num_st_levels_input .LT. 2 ) THEN
- CALL wrf_error_fatal ( 'Not enough soil temperature data for SSIB LSM scheme.')
- END IF
- !--------------------------------------------------------
- END SELECT enough_data
- interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
- CASE ( SLABSCHEME , LSMSCHEME, NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME, SSIBSCHEME )
- CALL process_soil_real ( grid%tsk , grid%tmn , grid%tavgsfc, &
- grid%landmask , grid%sst , grid%ht, grid%toposoil, &
- st_input , sm_input , sw_input , &
- st_levels_input , sm_levels_input , sw_levels_input , &
- grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
- flag_sst , flag_tavgsfc, flag_soilhgt,&
- flag_soil_layers, flag_soil_levels, &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte , &
- model_config_rec%sf_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
- ! surface_input_source=1 => use data from static file (fractional category as input)
- ! surface_input_source=2 => use data from grib file (dominant category as input)
- ! surface_input_source=3 => use dominant data from static file (dominant category as input)
- IF ( any_valid_points ) THEN
- IF ( config_flags%surface_input_source .EQ. 1 ) THEN
- ! Generate the vegetation and soil category information from the fractional input
- ! data, or use the existing dominant category fields if they exist.
- grid%vegcat (its,jts) = 0
- grid%soilcat(its,jts) = 0
- num_veg_cat = SIZE ( grid%landusef , DIM=2 )
- num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
- num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
- CALL process_percent_cat_new ( grid%landmask , &
- grid%landusef , grid%soilctop , grid%soilcbot , &
- grid%isltyp , grid%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)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%vegcat(i,j) = grid%ivgtyp(i,j)
- grid%soilcat(i,j) = grid%isltyp(i,j)
- END DO
- END DO
- ELSE IF ( config_flags%surface_input_source .EQ. 2 ) THEN
- ! Do we have dominant soil and veg data from the input already?
- IF ( grid%soilcat(i_valid,j_valid) .GT. 0.5 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
- END DO
- END DO
- END IF
- IF ( grid%vegcat(i_valid,j_valid) .GT. 0.5 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
- END DO
- END DO
- END IF
- ELSE IF ( config_flags%surface_input_source .EQ. 3 ) THEN
- ! Do we have dominant soil and veg data from the static input already?
- IF ( grid%sct_dom_gc(i_valid,j_valid) .GT. 0.5 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%isltyp(i,j) = NINT( grid%sct_dom_gc(i,j) )
- grid%soilcat(i,j) = grid%isltyp(i,j)
- END DO
- END DO
- ELSE
- WRITE ( a_message , * ) 'You have set surface_input_source = 3,'// &
- ' but your geogrid data does not have valid dominant soil data.'
- CALL wrf_error_fatal ( a_message )
- END IF
- IF ( grid%lu_index(i_valid,j_valid) .GT. 0.5 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%ivgtyp(i,j) = NINT( grid%lu_index(i,j) )
- grid%vegcat(i,j) = grid%ivgtyp(i,j)
- END DO
- END DO
- ELSE
- WRITE ( a_message , * ) 'You have set surface_input_source = 3,'//&
- ' but your geogrid data does not have valid dominant land use data.'
- CALL wrf_error_fatal ( a_message )
- END IF
- END IF
- END IF
- ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
- ! is for the 5-layer scheme.
- num_veg_cat = SIZE ( grid%landusef , DIM=2 )
- num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
- num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
- CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
- CALL nl_get_isice ( grid%id , grid%isice )
- CALL nl_get_iswater ( grid%id , grid%iswater )
- CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
- grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
- grid%soilcbot , grid%tmn , &
- grid%seaice_threshold , &
- config_flags%fractional_seaice, &
- num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
- grid%iswater , grid%isice , &
- model_config_rec%sf_surface_physics(grid%id) , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Land use assignment.
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%lu_index(i,j) = grid%ivgtyp(i,j)
- IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
- grid%landmask(i,j) = 1
- grid%xland(i,j) = 1
- ELSE
- grid%landmask(i,j) = 0
- grid%xland(i,j) = 2
- END IF
- END DO
- END DO
- ! Fix grid%tmn and grid%tsk.
- fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
- CASE ( SLABSCHEME , LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME, SSIBSCHEME )
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
- ( grid%sst(i,j) .GT. 170. ) .AND. ( grid%sst(i,j) .LT. 400. ) ) THEN
- grid%tmn(i,j) = grid%sst(i,j)
- grid%tsk(i,j) = grid%sst(i,j)
- ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
- grid%tmn(i,j) = grid%tsk(i,j)
- END IF
- END DO
- END DO
- END SELECT fix_tsk_tmn
- ! Is the grid%tsk reasonable?
- IF ( internal_time_loop .NE. 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
- grid%tsk(i,j) = grid%t_2(i,1,j)
- END IF
- END DO
- END DO
- ELSE
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
- print *,'error in the grid%tsk'
- print *,'i,j=',i,j
- print *,'grid%landmask=',grid%landmask(i,j)
- print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
- if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
- grid%tsk(i,j)=grid%tmn(i,j)
- else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
- grid%tsk(i,j)=grid%sst(i,j)
- else
- CALL wrf_error_fatal ( 'grid%tsk unreasonable' )
- end if
- END IF
- END DO
- END DO
- END IF
- ! Is the grid%tmn reasonable?
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
- .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
- IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .and. &
- ( model_config_rec%sf_surface_physics(grid%id) .NE. NOAHMPSCHEME ) ) THEN
- print *,'error in the grid%tmn'
- print *,'i,j=',i,j
- print *,'grid%landmask=',grid%landmask(i,j)
- print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
- END IF
- if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
- grid%tmn(i,j)=grid%tsk(i,j)
- else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
- grid%tmn(i,j)=grid%sst(i,j)
- else
- CALL wrf_error_fatal ( 'grid%tmn unreasonable' )
- endif
- END IF
- END DO
- END DO
-
- ! Minimum soil values, residual, from RUC LSM scheme. For input from Noah or EC, and using
- ! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
- ! input RUC data and using the Noah LSM scheme, this value must be added to the soil
- ! moisture input.
- lqmi(1:num_soil_top_cat) = &
- (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
- 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
- 0.004, 0.065 /)
- ! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
- ! At the initial time we care about values of soil moisture and temperature, other times are
- ! ignored by the model, so we ignore them, too.
- IF ( domain_ClockIsStartTime(grid) ) THEN
- account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
- CASE ( LSMSCHEME , NOAHMPSCHEME )
- iicount = 0
- IF ( flag_soil_layers == 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
- ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
- print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
- iicount = iicount + 1
- grid%smois(i,:,j) = 0.005
- END IF
- END DO
- END DO
- IF ( iicount .GT. 0 ) THEN
- print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
- END IF
- ELSE IF ( flag_soil_levels == 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) , 0.005 )
- ! grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) + lqmi(grid%isltyp(i,j)) , 0.005 )
- END DO
- END DO
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 170 ) .and. &
- ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
- print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
- iicount = iicount + 1
- grid%smois(i,:,j) = 0.005
- END IF
- END DO
- END DO
- IF ( iicount .GT. 0 ) THEN
- print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
- END IF
- END IF
- CASE ( RUCLSMSCHEME )
- iicount = 0
- IF ( flag_soil_layers == 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) , 0.005 )
- ! grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0.005 )
- END DO
- END DO
- ELSE IF ( flag_soil_levels == 1 ) THEN
- ! no op
- END IF
- CASE ( PXLSMSCHEME )
- iicount = 0
- IF ( flag_soil_layers == 1 ) THEN
- ! no op
- ELSE IF ( flag_soil_levels == 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) , 0.005 )
- ! grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) + lqmi(grid%isltyp(i,j)) , 0.005 )
- END DO
- END DO
- END IF
- END SELECT account_for_zero_soil_moisture
- END IF
- ! Is the grid%tslb reasonable?
- IF ( internal_time_loop .NE. 1 ) THEN
- DO j = jts, MIN(jde-1,jte)
- DO ns = 1 , model_config_rec%num_soil_layers
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
- grid%tslb(i,ns,j) = grid%t_2(i,1,j)
- grid%smois(i,ns,j) = 0.3
- END IF
- END DO
- END DO
- END DO
- ELSE
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
- ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
- IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
- ( model_config_rec%sf_surface_physics(grid%id) .NE. NOAHMPSCHEME ) .AND. &
- ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ).AND. &
- ( model_config_rec%sf_surface_physics(grid%id) .NE. SSIBSCHEME ).AND. & !fds
- ( model_config_rec%sf_surface_physics(grid%id) .NE. PXLSMSCHEME ) ) THEN
- print *,'error in the grid%tslb'
- print *,'i,j=',i,j
- print *,'grid%landmask=',grid%landmask(i,j)
- print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
- print *,'grid%tslb = ',grid%tslb(i,:,j)
- print *,'old grid%smois = ',grid%smois(i,:,j)
- grid%smois(i,1,j) = 0.3
- grid%smois(i,2,j) = 0.3
- grid%smois(i,3,j) = 0.3
- grid%smois(i,4,j) = 0.3
- END IF
- IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
- (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
- fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
- CASE ( SLABSCHEME )
- DO ns = 1 , model_config_rec%num_soil_layers
- grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
- grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
- END DO
- CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME, SSIBSCHEME )
- ! CALL wrf_error_fatal ( 'Assigned constant soil moisture to 0.3, stopping')
- DO ns = 1 , model_config_rec%num_soil_layers
- grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
- grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
- END DO
- END SELECT fake_soil_temp
- else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
- CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
- DO ns = 1 , model_config_rec%num_soil_layers
- grid%tslb(i,ns,j)=grid%tsk(i,j)
- END DO
- else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
- CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
- DO ns = 1 , model_config_rec%num_soil_layers
- grid%tslb(i,ns,j)=grid%sst(i,j)
- END DO
- else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
- CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
- DO ns = 1 , model_config_rec%num_soil_layers
- grid%tslb(i,ns,j)=grid%tmn(i,j)
- END DO
- else
- CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
- endif
- END IF
- END DO
- END DO
- END IF
- ! Adjustments for the seaice field AFTER the grid%tslb computations. This is
- ! is for the Noah LSM scheme.
- num_veg_cat = SIZE ( grid%landusef , DIM=2 )
- num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
- num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
- CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
- CALL nl_get_isice ( grid%id , grid%isice )
- CALL nl_get_iswater ( grid%id , grid%iswater )
- CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
- grid%ivgtyp , grid%vegcat , grid%lu_index , &
- grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
- grid%soilctop , &
- grid%soilcbot , grid%tmn , grid%vegfra , &
- grid%tslb , grid%smois , grid%sh2o , &
- grid%seaice_threshold , &
- grid%sst,flag_sst, &
- config_flags%fractional_seaice, &
- num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
- model_config_rec%num_soil_layers , &
- grid%iswater , grid%isice , &
- model_config_rec%sf_surface_physics(grid%id) , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Let us make sure (again) that the grid%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 ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
- ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
- ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
- ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
- IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
- oops1=oops1+1
- grid%ivgtyp(i,j) = 5
- grid%isltyp(i,j) = 8
- grid%landmask(i,j) = 1
- grid%xland(i,j) = 1
- ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
- oops2=oops2+1
- grid%ivgtyp(i,j) = config_flags%iswater
- grid%isltyp(i,j) = 14
- grid%landmask(i,j) = 0
- grid%xland(i,j) = 2
- ELSE
- print *,'the grid%landmask and soil/veg cats do not match'
- print *,'i,j=',i,j
- print *,'grid%landmask=',grid%landmask(i,j)
- print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
- print *,'grid%isltyp=',grid%isltyp(i,j)
- print *,'iswater=', config_flags%iswater
- print *,'grid%tslb=',grid%tslb(i,:,j)
- print *,'grid%sst=',grid%sst(i,j)
- CALL wrf_error_fatal ( '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
- ! fill grid%sst array with grid%tsk if missing in real input (needed for time-varying grid%sst in wrf)
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( flag_sst .NE. 1 ) THEN
- grid%sst(i,j) = grid%tsk(i,j)
- ENDIF
- END DO
- END DO
- !tgs set snoalb to land value if the water point is covered with ice
- DO j = jts, MIN(jde-1,jte)
- DO i = its, MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( grid%ivgtyp(i,j) .EQ. config_flags%isice) THEN
- grid%snoalb(i,j) = 0.75
- ENDIF
- END DO
- END DO
- ! 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.
- ! Check whether grid%znw (full level) data are truly full levels. If not, we need to adjust them
- ! to be full levels.
- ! in this test, we check if grid%znw(1) is neither 0 nor 1 (within a tolerance of 10**-5)
- were_bad = .false.
- IF ( ( (grid%znw(1).LT.(1-1.E-5) ) .OR. ( grid%znw(1).GT.(1+1.E-5) ) ).AND. &
- ( (grid%znw(1).LT.(0-1.E-5) ) .OR. ( grid%znw(1).GT.(0+1.E-5) ) ) ) THEN
- were_bad = .true.
- print *,'Your grid%znw input values are probably half-levels. '
- print *,grid%znw
- print *,'WRF expects grid%znw values to be full levels. '
- print *,'Adjusting now to full levels...'
- ! We want to ignore the first value if it's negative
- IF (grid%znw(1).LT.0) THEN
- grid%znw(1)=0
- END IF
- DO k=2,kde
- grid%znw(k)=2*grid%znw(k)-grid%znw(k-1)
- END DO
- END IF
- ! Let's check our changes
- IF ( ( ( grid%znw(1) .LT. (1-1.E-5) ) .OR. ( grid%znw(1) .GT. (1+1.E-5) ) ).AND. &
- ( ( grid%znw(1) .LT. (0-1.E-5) ) .OR. ( grid%znw(1) .GT. (0+1.E-5) ) ) ) THEN
- print *,'The input grid%znw height values were half-levels or erroneous. '
- print *,'Attempts to treat the values as half-levels and change them '
- print *,'to valid full levels failed.'
- CALL wrf_error_fatal("bad grid%znw values from input files")
- ELSE IF ( were_bad ) THEN
- print *,'...adjusted. grid%znw array now contains full eta level values. '
- ENDIF
- IF ( grid%znw(1) .LT. grid%znw(kde) ) THEN
- DO k=1, kde/2
- hold_znw = grid%znw(k)
- grid%znw(k)=grid%znw(kde+1-k)
- grid%znw(kde+1-k)=hold_znw
- END DO
- END IF
- DO k=1, kde-1
- grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
- grid%rdnw(k) = 1./grid%dnw(k)
- grid%znu(k) = 0.5*(grid%znw(k+1)+grid%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
- grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
- grid%rdn(k) = 1./grid%dn(k)
- grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k)
- grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
- END DO
- ! Scads of vertical coefficients.
- cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
- cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
- grid%cf1 = grid%fnp(2) + cof1
- grid%cf2 = grid%fnm(2) - cof1 - cof2
- grid%cf3 = cof2
- grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
- grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
- ! Inverse grid distances.
- grid%rdx = 1./config_flags%dx
- grid%rdy = 1./config_flags%dy
- ! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
- ! and grid%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
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%ph0(i,1,j) = grid%ht(i,j) * g
- grid%ph_2(i,1,j) = 0.
- END DO
- END DO
- ! 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)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
-
- ! 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*grid%ht(i,j)/a/r_d ) **0.5 )
- DO k = 1, kte-1
- grid%php(i,k,j) = grid%znw(k)*(p_surf - grid%p_top) + grid%p_top ! temporary, full lev base pressure
- grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
- temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
- ! temp = t00 + A*LOG(grid%pb(i,k,j)/p00)
- grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
- grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
- END DO
- ! Base state mu is defined as base state surface pressure minus grid%p_top
- grid%mub(i,j) = p_surf - grid%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 = grid%mu0(i,j) + grid%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.
- grid%phb(i,1,j) = grid%ht(i,j) * g
- IF (grid%hypsometric_opt == 1) THEN
- DO k = 2,kte
- grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
- END DO
- ELSE IF (grid%hypsometric_opt == 2) THEN
- DO k = 2,kte
- pfu = grid%mub(i,j)*grid%znw(k) + grid%p_top
- pfd = grid%mub(i,j)*grid%znw(k-1) + grid%p_top
- phm = grid%mub(i,j)*grid%znu(k-1) + grid%p_top
- grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
- END DO
- ELSE
- CALL wrf_error_fatal( 'initialize_real: hypsometric_opt should be 1 or 2' )
- END IF
- 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)
- grid%mub(i,j) = grid%mub(i-1,j)
- grid%mu_2(i,j) = grid%mu_2(i-1,j)
- DO k = 1, kte-1
- grid%pb(i,k,j) = grid%pb(i-1,k,j)
- grid%t_init(i,k,j) = grid%t_init(i-1,k,j)
- grid%alb(i,k,j) = grid%alb(i-1,k,j)
- END DO
- DO k = 1, kte
- grid%phb(i,k,j) = grid%phb(i-1,k,j)
- END DO
- END DO
- END IF
- IF ( jte .EQ. jde ) THEN
- j = jde
- DO i = its, ite
- grid%mub(i,j) = grid%mub(i,j-1)
- grid%mu_2(i,j) = grid%mu_2(i,j-1)
- DO k = 1, kte-1
- grid%pb(i,k,j) = grid%pb(i,k,j-1)
- grid%t_init(i,k,j) = grid%t_init(i,k,j-1)
- grid%alb(i,k,j) = grid%alb(i,k,j-1)
- END DO
- DO k = 1, kte
- grid%phb(i,k,j) = grid%phb(i,k,j-1)
- END DO
- END DO
- END IF
- ! Compute the perturbation dry pressure (grid%mub + grid%mu_2 + ptop = dry grid%psfc).
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%mu_2(i,j) = grid%mu0(i,j) - grid%mub(i,j)
- 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)
- grid%mu_2(i,j) = grid%mu_2(i-1,j)
- END DO
- END IF
- IF ( jte .EQ. jde ) THEN
- j = jde
- DO i = its, ite
- grid%mu_2(i,j) = grid%mu_2(i,j-1)
- END DO
- END IF
- lev500 = 0
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! Assign the potential temperature (perturbation from t0) and qv on all the mass
- ! point locations.
- DO k = 1 , kde-1
- grid%t_2(i,k,j) = grid%t_2(i,k,j) - t0
- END DO
- dpmu = 10001.
- loop_count = 0
- DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
- ( loop_count .LT. 5 ) )
- loop_count = loop_count + 1
- ! 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(i,k,j,P_QV)+moist(i,k,j,P_QV))
- qvf2 = 1./(1.+qvf1)
- qvf1 = qvf1*qvf2
- grid%p(i,k,j) = - 0.5*(grid%mu_2(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
- qvf = 1. + rvovrd*moist(i,k,j,P_QV)
- grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf&
- *(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
- grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
- grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(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(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
- qvf2 = 1./(1.+qvf1)
- qvf1 = qvf1*qvf2
- grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_2(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
- qvf = 1. + rvovrd*moist(i,k,j,P_QV)
- grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
- (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
- grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
- grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
- END DO
- #if 1
- ! This is the hydrostatic equation used in the model after the small timesteps. In
- ! the model, grid%al (inverse density) is computed from the geopotential.
- IF (grid%hypsometric_opt == 1) THEN
- DO k = 2,kte
- grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
- grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
- + grid%mu_2(i,j)*grid%alb(i,k-1,j) )
- grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
- END DO
- ELSE IF (grid%hypsometric_opt == 2) THEN
- ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
- ! Note that al*p approximates Rd*T and dLOG(p) does z.
- ! Here T varies mostly linear with z, the first-order integration produces better result.
- grid%ph_2(i,1,j) = grid%phb(i,1,j)
- DO k = 2,kte
- pfu = grid%mu0(i,j)*grid%znw(k) + grid%p_top
- pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
- phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
- grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
- END DO
- DO k = 1,kte
- grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
- END DO
- END IF
- #else
- ! Get the perturbation geopotential from the 3d height array from WPS.
- DO k = 2,kte
- grid%ph_2(i,k,j) = grid%ph0(i,k,j)*g - grid%phb(i,k,j)
- END DO
- #endif
- ! Adjust the column pressure so that the computed 500 mb height is close to the
- ! input value (of course, not when we are doing hybrid input).
- IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. i_valid ) .AND. ( j .EQ. j_valid ) ) THEN
- DO k = 1 , num_metgrid_levels
- IF ( ABS ( grid%p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
- lev500 = k
- EXIT
- END IF
- END DO
- END IF
- ! We only do the adjustment of height if we have the input data on pressure
- ! surfaces, and folks have asked to do this option.
- IF ( ( flag_metgrid .EQ. 1 ) .AND. &
- ( flag_ptheta .EQ. 0 ) .AND. &
- ( config_flags%adjust_heights ) .AND. &
- ( lev500 .NE. 0 ) ) THEN
- DO k = 2 , kte-1
- ! Get the pressures on the full eta levels (grid%php is defined above as
- ! the full-lev base pressure, an easy array to use for 3d space).
- pl = grid%php(i,k ,j) + &
- ( grid%p(i,k-1 ,j) * ( grid%znw(k ) - grid%znu(k ) ) + &
- grid%p(i,k ,j) * ( grid%znu(k-1 ) - grid%znw(k ) ) ) / &
- ( grid%znu(k-1 ) - grid%znu(k ) )
- pu = grid%php(i,k+1,j) + &
- ( grid%p(i,k-1+1,j) * ( grid%znw(k +1) - grid%znu(k+1) ) + &
- grid%p(i,k +1,j) * ( grid%znu(k-1+1) - grid%znw(k+1) ) ) / &
- ( grid%znu(k-1+1) - grid%znu(k+1) )
- ! If these pressure levels trap 500 mb, use them to interpolate
- ! to the 500 mb level of the computed height.
- IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
- zl = ( grid%ph_2(i,k ,j) + grid%phb(i,k ,j) ) / g
- zu = ( grid%ph_2(i,k+1,j) + grid%phb(i,k+1,j) ) / g
- z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
- zu * ( LOG(pl ) - LOG(50000.) ) ) / &
- ( LOG(pl) - LOG(pu) )
- ! z500 = ( zl * ( (50000.) - (pu ) ) + &
- ! zu * ( (pl ) - (50000.) ) ) / &
- ! ( (pl) - (pu) )
- ! Compute the difference of the 500 mb heights (computed minus input), and
- ! then the change in grid%mu_2. The grid%php is still full-levels, base pressure.
- dz500 = z500 - grid%ght_gc(i,lev500,j)
- tvsfc = ((grid%t_2(i,1,j)+t0)*((grid%p(i,1,j)+grid%php(i,1,j))/p1000mb)**(r_d/cp)) * &
- (1.+0.6*moist(i,1,j,P_QV))
- dpmu = ( grid%php(i,1,j) + grid%p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
- dpmu = dpmu - ( grid%php(i,1,j) + grid%p(i,1,j) )
- grid%mu_2(i,j) = grid%mu_2(i,j) - dpmu
- EXIT
- END IF
- END DO
- ELSE
- dpmu = 0.
- END IF
- END DO
- END DO
- END DO
-
- ! If this is data from the SI, then we probably do not have the original
- ! surface data laying around. Note that these are all the lowest levels
- ! of the respective 3d arrays. For surface pressure, we assume that the
- ! vertical gradient of grid%p prime is zilch. This is not all that important.
- ! These are filled in so that the various plotting routines have something
- ! to play with at the initial time for the model.
- IF ( flag_metgrid .NE. 1 ) THEN
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%u10(i,j)=grid%u_2(i,1,j)
- END DO
- END DO
- DO j = jts, min(jde,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- grid%v10(i,j)=grid%v_2(i,1,j)
- END DO
- END DO
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
- grid%psfc(i,j)=p_surf + grid%p(i,1,j)
- grid%q2(i,j)=moist(i,1,j,P_QV)
- grid%th2(i,j)=grid%t_2(i,1,j)+300.
- grid%t2(i,j)=grid%th2(i,j)*(((grid%p(i,1,j)+grid%pb(i,1,j))/p00)**(r_d/cp))
- END DO
- END DO
- ! If this data is from WPS, then we have previously assigned the surface
- ! data for u, v, and t. If we have an input qv, welp, we assigned that one,
- ! too. Now we pick up the left overs, and if RH came in - we assign the
- ! mixing ratio.
- ELSE IF ( flag_metgrid .EQ. 1 ) THEN
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
- ! grid%psfc(i,j)=p_surf + grid%p(i,1,j)
- grid%th2(i,j)=grid%t2(i,j)*(p00/(grid%p(i,1,j)+grid%pb(i,1,j)))**(r_d/cp)
- END DO
- END DO
- IF ( flag_qv .NE. 1 ) THEN
- DO j = jts, min(jde-1,jte)
- DO i = its, min(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! grid%q2(i,j)=moist(i,1,j,P_QV)
- grid%q2(i,j)=grid%qv_gc(i,1,j)
- END DO
- END DO
- END IF
- END IF
- CALL cpu_time(t_end)
- ! Set flag to denote that we are saving original values of HT, MUB, and
- ! PHB for 2-way nesting and cycling.
- grid%save_topo_from_real=1
- ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
- #ifdef DM_PARALLEL
- # include "HALO_EM_INIT_1.inc"
- # include "HALO_EM_INIT_2.inc"
- # include "HALO_EM_INIT_3.inc"
- # include "HALO_EM_INIT_4.inc"
- # include "HALO_EM_INIT_5.inc"
- #endif
- RETURN
- END SUBROUTINE init_domain_rk
- !---------------------------------------------------------------------
- SUBROUTINE const_module_initialize ( p00 , t00 , a , tiso )
- USE module_configure
- IMPLICIT NONE
- ! For the real-data-cases only.
- REAL , INTENT(OUT) :: p00 , t00 , a , tiso
- CALL nl_get_base_pres ( 1 , p00 )
- CALL nl_get_base_temp ( 1 , t00 )
- CALL nl_get_base_lapse ( 1 , a )
- CALL nl_get_iso_temp ( 1 , tiso )
- END SUBROUTINE const_module_initialize
- !-------------------------------------------------------------------
- SUBROUTINE rebalance_driver ( grid )
- IMPLICIT NONE
- TYPE (domain) :: grid
- CALL rebalance( grid &
- !
- #include "actual_new_args.inc"
- !
- )
- END SUBROUTINE rebalance_driver
- !---------------------------------------------------------------------
- SUBROUTINE rebalance ( grid &
- !
- #include "dummy_new_args.inc"
- !
- )
- IMPLICIT NONE
- TYPE (domain) :: grid
- #include "dummy_new_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 , tiso
- 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, &
- ips, ipe, jps, jpe, kps, kpe, &
- i, j, k
- REAL :: temp, temp_int
- REAL :: pfu, pfd, phm
- REAL :: w1, w2, z0, z1, z2
- 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) )
- ! Fill config_flags the options for a particular domain
- CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
- ! Some of the many weird geopotential initializations that we'll see today: grid%ph0 is total,
- ! and grid%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
- grid%ph0(i,1,j) = grid%ht_fine(i,j) * g
- grid%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), and constant stratosphere
- ! temp (tiso, K) either from input file or from namelist (for backward compatibiliy).
- IF ( config_flags%use_baseparam_fr_nml ) then
- ! get these from namelist
- CALL wrf_message('ndown: using namelist constants')
- CALL const_module_initialize ( p00 , t00 , a , tiso )
- ELSE
- ! get these constants from model data
- CALL wrf_message('ndown: using constants from file')
- t00 = grid%t00
- p00 = grid%p00
- a = grid%tlp
- tiso = grid%tiso
- IF (t00 .LT. 100. .or. p00 .LT. 10000.) THEN
- WRITE(wrf_err_message,*)&
- 'ndown_em: did not find base state parameters in wrfout. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
- CALL wrf_error_fatal(TRIM(wrf_err_message))
- ENDIF
- ENDIF
- ! 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 grid%ht.
- p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht_fine(i,j)/a/r_d ) **0.5 )
- p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j) /a/r_d ) **0.5 )
- DO k = 1, kte-1
- grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
- pb_int = grid%znu(k)*(p_surf_int - grid%p_top) + grid%p_top
- temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
- ! temp = t00 + A*LOG(pb/p00)
- grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
- ! grid%t_init(i,k,j) = (t00 + A*LOG(grid%pb(i,k,j)/p00))*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
- temp_int = MAX ( tiso, t00 + A*LOG(pb_int /p00) )
- t_init_int(i,k,j)= temp_int*(p00/pb_int )**(r_d/cp) - t0
- ! t_init_int(i,k,j)= (t00 + A*LOG(pb_int /p00))*(p00/pb_int )**(r_d/cp) - t0
- grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
- END DO
- ! Base state mu is defined as base state surface pressure minus grid%p_top
- grid%mub(i,j) = p_surf - grid%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 = ( grid%mub(i,j) + grid%mu_2(i,j) ) + grid%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.
- grid%phb(i,1,j) = grid%ht_fine(i,j) * g
- IF (grid%hypsometric_opt == 1) THEN
- DO k = 2,kte
- grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
- END DO
- ELSE IF (grid%hypsometric_opt == 2) THEN
- DO k = 2,kte
- pfu = grid%mub(i,j)*grid%znw(k) + grid%p_top
- pfd = grid%mub(i,j)*grid%znw(k-1) + grid%p_top
- phm = grid%mub(i,j)*grid%znu(k-1) + grid%p_top
- grid%phb(i,k,j) = grid%phb(i,k-1,j) + grid%alb(i,k-1,j)*phm*LOG(pfd/pfu)
- END DO
- ELSE
- CALL wrf_error_fatal( 'initialize_real: hypsometric_opt should be 1 or 2' )
- END IF
- 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)
- grid%ht(i,j) = grid%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
- grid%t_2(i,k,j) = grid%t_2(i,k,j) + ( grid%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(i,k,j,P_QV)+moist(i,k,j,P_QV))
- qvf2 = 1./(1.+qvf1)
- qvf1 = qvf1*qvf2
- grid%p(i,k,j) = - 0.5*(grid%mu_2(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
- qvf = 1. + rvovrd*moist(i,k,j,P_QV)
- grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
- (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
- grid%al(i,k,j) = grid%alt(i,k,j) - grid%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(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
- qvf2 = 1./(1.+qvf1)
- qvf1 = qvf1*qvf2
- grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_2(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
- qvf = 1. + rvovrd*moist(i,k,j,P_QV)
- grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_2(i,k,j)+t0)*qvf* &
- (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
- grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
- END DO
- ! This is the hydrostatic equation used in the model after the small timesteps. In
- ! the model, grid%al (inverse density) is computed from the geopotential.
- IF (grid%hypsometric_opt == 1) THEN
- DO k = 2,kte
- grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) - &
- grid%dnw(k-1) * ( (grid%mub(i,j)+grid%mu_2(i,j))*grid%al(i,k-1,j) &
- + grid%mu_2(i,j)*grid%alb(i,k-1,j) )
- grid%ph0(i,k,j) = grid%ph_2(i,k,j) + grid%phb(i,k,j)
- END DO
- ELSE IF (grid%hypsometric_opt == 2) THEN
- ! Alternative hydrostatic eq.: dZ = -al*p*dLOG(p), where p is dry pressure.
- ! Note that al*p approximates Rd*T and dLOG(p) does z.
- ! Here T varies mostly linear with z, the first-order integration produces better result.
- grid%ph_2(i,1,j) = grid%phb(i,1,j)
- DO k = 2,kte
- pfu = grid%mu0(i,j)*grid%znw(k) + grid%p_top
- pfd = grid%mu0(i,j)*grid%znw(k-1) + grid%p_top
- phm = grid%mu0(i,j)*grid%znu(k-1) + grid%p_top
- grid%ph_2(i,k,j) = grid%ph_2(i,k-1,j) + grid%alt(i,k-1,j)*phm*LOG(pfd/pfu)
- END DO
- DO k = 1,kte
- grid%ph_2(i,k,j) = grid%ph_2(i,k,j) - grid%phb(i,k,j)
- END DO
- END IF
- ! update psfc in fine grid
- z0 = grid%ph0(i,1,j)/g
- z1 = 0.5*(grid%ph0(i,1,j)+grid%ph0(i,2,j))/g
- z2 = 0.5*(grid%ph0(i,2,j)+grid%ph0(i,3,j))/g
- w1 = (z0 - z2)/(z1 - z2)
- w2 = 1. - w1
- grid%psfc(i,j) = w1*(grid%p(i,1,j)+grid%pb(i,1,j))+w2*(grid%p(i,2,j)+grid%pb(i,2,j))
- END DO
- END DO
- DEALLOCATE ( t_init_int )
- ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
- #ifdef DM_PARALLEL
- # include "HALO_EM_INIT_1.inc"
- # include "HALO_EM_INIT_2.inc"
- # include "HALO_EM_INIT_3.inc"
- # include "HALO_EM_INIT_4.inc"
- # include "HALO_EM_INIT_5.inc"
- #endif
- END SUBROUTINE rebalance
- !---------------------------------------------------------------------
- RECURSIVE SUBROUTINE find_my_parent ( grid_ptr_in , grid_ptr_out , id_i_am , id_wanted , found_the_id )
- ! RAR - Modified to correct problem in which the parent of a child domain could
- ! not be found in the namelist. This condition typically occurs while using the
- ! "allow_grid" namelist option when an inactive domain comes before an active
- ! domain in the list, i.e., the domain number of the active domain is greater than
- ! that of an inactive domain at the same level.
- !
- USE module_domain
- TYPE(domain) , POINTER :: grid_ptr_in , grid_ptr_out
- TYPE(domain) , POINTER :: grid_ptr_sibling
- INTEGER :: id_wanted , id_i_am
- INTEGER :: nest ! RAR
- LOGICAL :: found_the_id
- found_the_id = .FALSE.
- grid_ptr_sibling => grid_ptr_in
- nest = 0 ! RAR
- DO WHILE ( ASSOCIATED ( grid_ptr_sibling ) )
- IF ( grid_ptr_sibling%grid_id .EQ. id_wanted ) THEN
- found_the_id = .TRUE.
- grid_ptr_out => grid_ptr_sibling
- RETURN
- ! RAR ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 ) THEN
- ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 .AND. nest .LT. grid_ptr_sibling%num_nests ) THEN
- nest = nest + 1 ! RAR
- grid_ptr_sibling => grid_ptr_sibling%nests(nest)%ptr ! RAR
- CALL find_my_parent ( grid_ptr_sibling , grid_ptr_out , id_i_am , id_wanted , found_the_id )
- IF (.NOT. found_the_id) grid_ptr_sibling => grid_ptr_sibling%parents(1)%ptr ! RAR
- ELSE
- grid_ptr_sibling => grid_ptr_sibling%sibling
- END IF
- END DO
- END SUBROUTINE find_my_parent
- !---------------------------------------------------------------------
- RECURSIVE SUBROUTINE find_my_parent2 ( grid_ptr_in , grid_ptr_out , id_wanted , found_the_id )
- USE module_domain
- TYPE(domain) , POINTER :: grid_ptr_in
- TYPE(domain) , POINTER :: grid_ptr_out
- INTEGER , INTENT(IN ) :: id_wanted
- LOGICAL , INTENT(OUT) :: found_the_id
- ! Local
- TYPE(domain) , POINTER :: grid_ptr_holder
- INTEGER :: kid
- ! Initializations
- found_the_id = .FALSE.
- grid_ptr_holder => grid_ptr_in
- ! Have we found the correct location? If so, we can just pop back up with
- ! the pointer to the right location (i.e. the parent), thank you very much.
- IF ( id_wanted .EQ. grid_ptr_in%grid_id ) THEN
- found_the_id = .TRUE.
- grid_ptr_out => grid_ptr_in
- ! We gotta keep looking.
- ELSE
- ! We drill down and process each nest from this domain. We don't have to
- ! worry about siblings, as we are running over all of the kids for this parent,
- ! so it amounts to the same set of domains being tested.
- loop_over_all_kids : DO kid = 1 , grid_ptr_in%num_nests
- IF ( ASSOCIATED ( grid_ptr_in%nests(kid)%ptr ) ) THEN
- CALL find_my_parent2 ( grid_ptr_in%nests(kid)%ptr , grid_ptr_out , id_wanted , found_the_id )
- IF ( found_the_id ) THEN
- EXIT loop_over_all_kids
- END IF
- END IF
- END DO loop_over_all_kids
- END IF
- END SUBROUTINE find_my_parent2
- #endif
- !---------------------------------------------------------------------
- #ifdef VERT_UNIT
- !This is a main program for a small unit test for the vertical interpolation.
- program vint
- implicit none
- integer , parameter :: ij = 3
- integer , parameter :: keta = 30
- integer , parameter :: kgen =20
- integer :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- integer :: generic
- real , dimension(1:ij,kgen,1:ij) :: fo , po
- real , dimension(1:ij,1:keta,1:ij) :: fn_calc , fn_interp , pn
- integer, parameter :: interp_type = 1 ! 2
- ! integer, parameter :: lagrange_order = 2 ! 1
- integer :: lagrange_order
- logical, parameter :: lowest_lev_from_sfc = .FALSE. ! .TRUE.
- logical, parameter :: use_levels_below_ground = .FALSE. ! .TRUE.
- logical, parameter :: use_surface = .FALSE. ! .TRUE.
- real , parameter :: zap_close_levels = 500. ! 100.
- integer, parameter :: force_sfc_in_vinterp = 0 ! 6
- integer :: k
- ids = 1 ; ide = ij ; jds = 1 ; jde = ij ; kds = 1 ; kde = keta
- ims = 1 ; ime = ij ; jms = 1 ; jme = ij ; kms = 1 ; kme = keta
- its = 1 ; ite = ij ; jts = 1 ; jte = ij ; kts = 1 ; kte = keta
- generic = kgen
- print *,' '
- print *,'------------------------------------'
- print *,'UNIT TEST FOR VERTICAL INTERPOLATION'
- print *,'------------------------------------'
- print *,' '
- do lagrange_order = 1 , 2
- print *,' '
- print *,'------------------------------------'
- print *,'Lagrange Order = ',lagrange_order
- print *,'------------------------------------'
- print *,' '
- call fillitup ( fo , po , fn_calc , pn , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte , &
- generic , lagrange_order )
- print *,' '
- print *,'Level Pressure Field'
- print *,' (Pa) (generic)'
- print *,'------------------------------------'
- print *,' '
- do k = 1 , generic
- write (*,fmt='(i2,2x,f12.3,1x,g15.8)' ) &
- k,po(2,k,2),fo(2,k,2)
- end do
- print *,' '
- call vert_interp ( fo , po , fn_interp , pn , &
- generic , 'T' , &
- interp_type , lagrange_order , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- print *,'Multi-Order Interpolator'
- print *,'------------------------------------'
- print *,' '
- print *,'Level Pressure Field Field Field'
- print *,' (Pa) Calc Interp Diff'
- print *,'------------------------------------'
- print *,' '
- do k = kts , kte-1
- write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
- k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
- end do
- call vert_interp_old ( fo , po , fn_interp , pn , &
- generic , 'T' , &
- interp_type , lagrange_order , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- print *,'Linear Interpolator'
- print *,'------------------------------------'
- print *,' '
- print *,'Level Pressure Field Field Field'
- print *,' (Pa) Calc Interp Diff'
- print *,'------------------------------------'
- print *,' '
- do k = kts , kte-1
- write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
- k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
- end do
- end do
- end program vint
- subroutine wrf_error_fatal (string)
- character (len=*) :: string
- print *,string
- stop
- end subroutine wrf_error_fatal
- subroutine fillitup ( fo , po , fn , pn , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte , &
- generic , lagrange_order )
- implicit none
- integer , intent(in) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- integer , intent(in) :: generic , lagrange_order
- real , dimension(ims:ime,generic,jms:jme) , intent(out) :: fo , po
- real , dimension(ims:ime,kms:kme,jms:jme) , intent(out) :: fn , pn
- integer :: i , j , k
- real , parameter :: piov2 = 3.14159265358 / 2.
- k = 1
- do j = jts , jte
- do i = its , ite
- po(i,k,j) = 102000.
- end do
- end do
- do k = 2 , generic
- do j = jts , jte
- do i = its , ite
- po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) )
- end do
- end do
- end do
- if ( lagrange_order .eq. 1 ) then
- do k = 1 , generic
- do j = jts , jte
- do i = its , ite
- fo(i,k,j) = po(i,k,j)
- ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
- end do
- end do
- end do
- else if ( lagrange_order .eq. 2 ) then
- do k = 1 , generic
- do j = jts , jte
- do i = its , ite
- fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)*((102000.-po(i,k,j))/102000.))*102000.
- ! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
- end do
- end do
- end do
- end if
- !!!!!!!!!!!!
- do k = kts , kte
- do j = jts , jte
- do i = its , ite
- pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) )
- end do
- end do
- end do
- do k = kts , kte-1
- do j = jts , jte
- do i = its , ite
- pn(i,k,j) = ( pn(i,k,j) + pn(i,k+1,j) ) /2.
- end do
- end do
- end do
- if ( lagrange_order .eq. 1 ) then
- do k = kts , kte-1
- do j = jts , jte
- do i = its , ite
- fn(i,k,j) = pn(i,k,j)
- ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
- end do
- end do
- end do
- else if ( lagrange_order .eq. 2 ) then
- do k = kts , kte-1
- do j = jts , jte
- do i = its , ite
- fn(i,k,j) = (((pn(i,k,j)-5000.)/102000.)*((102000.-pn(i,k,j))/102000.))*102000.
- ! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
- end do
- end do
- end do
- end if
- end subroutine fillitup
- #endif
- !---------------------------------------------------------------------
- SUBROUTINE vert_interp ( fo , po , fnew , pnu , &
- generic , var_type , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Vertically interpolate the new field. The original field on the original
- ! pressure levels is provided, and the new pressure surfaces to interpolate to.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
- LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
- REAL , INTENT(IN) :: zap_close_levels
- INTEGER , INTENT(IN) :: force_sfc_in_vinterp
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- INTEGER , INTENT(IN) :: generic
- CHARACTER (LEN=1) :: var_type
- REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: fo , po
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
- REAL , DIMENSION(ims:ime,generic,jms:jme) :: forig , porig
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
- ! Local vars
- INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 , knext
- INTEGER :: istart , iend , jstart , jend , kstart , kend
- INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
- INTEGER , DIMENSION(ims:ime ) :: ks
- INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
- INTEGER :: count , zap , zap_below , zap_above , kst , kcount
- INTEGER :: kinterp_start , kinterp_end , sfc_level
- LOGICAL :: any_below_ground
- REAL :: p1 , p2 , pn, hold
- REAL , DIMENSION(1:generic) :: ordered_porig , ordered_forig
- REAL , DIMENSION(kts:kte) :: ordered_pnew , ordered_fnew
-
- ! Excluded middle.
- LOGICAL :: any_valid_points
- INTEGER :: i_valid , j_valid
- LOGICAL :: flip_data_required
- ! Horiontal loop bounds for different variable types.
- IF ( var_type .EQ. 'U' ) THEN
- istart = its
- iend = ite
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = MAX(ids+1,its) , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
- END DO
- END DO
- IF ( ids .EQ. its ) THEN
- DO k = 1,generic
- porig(its,k,j) = po(its,k,j)
- END DO
- END IF
- IF ( ide .EQ. ite ) THEN
- DO k = 1,generic
- porig(ite,k,j) = po(ite-1,k,j)
- END DO
- END IF
- DO k = kstart,kend
- DO i = MAX(ids+1,its) , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
- END DO
- END DO
- IF ( ids .EQ. its ) THEN
- DO k = kstart,kend
- pnew(its,k,j) = pnu(its,k,j)
- END DO
- END IF
- IF ( ide .EQ. ite ) THEN
- DO k = kstart,kend
- pnew(ite,k,j) = pnu(ite-1,k,j)
- END DO
- END IF
- END DO
- ELSE IF ( var_type .EQ. 'V' ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = jte
- kstart = kts
- kend = kte-1
- DO i = istart,iend
- DO k = 1,generic
- DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
- END DO
- END DO
- IF ( jds .EQ. jts ) THEN
- DO k = 1,generic
- porig(i,k,jts) = po(i,k,jts)
- END DO
- END IF
- IF ( jde .EQ. jte ) THEN
- DO k = 1,generic
- porig(i,k,jte) = po(i,k,jte-1)
- END DO
- END IF
- DO k = kstart,kend
- DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
- END DO
- END DO
- IF ( jds .EQ. jts ) THEN
- DO k = kstart,kend
- pnew(i,k,jts) = pnu(i,k,jts)
- END DO
- END IF
- IF ( jde .EQ. jte ) THEN
- DO k = kstart,kend
- pnew(i,k,jte) = pnu(i,k,jte-1)
- END DO
- END IF
- END DO
- ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- ELSE
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- END IF
- ! We need to find if there are any valid non-excluded-middle points in this
- ! tile. If so, then we need to hang on to a valid i,j location.
- any_valid_points = .false.
- find_valid : DO j = jstart , jend
- DO i = istart , iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- any_valid_points = .true.
- i_valid = i
- j_valid = j
- EXIT find_valid
- END DO
- END DO find_valid
- IF ( .NOT. any_valid_points ) THEN
- RETURN
- END IF
- IF ( porig(i_valid,2,j_valid) .LT. porig(i_valid,generic,j_valid) ) THEN
- flip_data_required = .true.
- ELSE
- flip_data_required = .false.
- END IF
- DO j = jstart , jend
- ! The lowest level is the surface. Levels 2 through "generic" are supposed to
- ! be "bottom-up". Flip if they are not. This is based on the input pressure
- ! array.
- IF ( flip_data_required ) THEN
- DO kn = 2 , ( generic + 1 ) / 2
- DO i = istart , iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- hold = porig(i,kn,j)
- porig(i,kn,j) = porig(i,generic+2-kn,j)
- porig(i,generic+2-kn,j) = hold
- forig(i,kn,j) = fo (i,generic+2-kn,j)
- forig(i,generic+2-kn,j) = fo (i,kn,j)
- END DO
- END DO
- DO i = istart , iend
- forig(i,1,j) = fo (i,1,j)
- END DO
- IF ( MOD(generic,2) .EQ. 0 ) THEN
- k=generic/2 + 1
- DO i = istart , iend
- forig(i,k,j) = fo (i,k,j)
- END DO
- END IF
- ELSE
- DO kn = 1 , generic
- DO i = istart , iend
- forig(i,kn,j) = fo (i,kn,j)
- END DO
- END DO
- END IF
- ! Skip all of the levels below ground in the original data based upon the surface pressure.
- ! The ko_above_sfc is the index in the pressure array that is above the surface. If there
- ! are no levels underground, this is index = 2. The remaining levels are eligible for use
- ! in the vertical interpolation.
- DO i = istart , iend
- ko_above_sfc(i) = -1
- END DO
- DO ko = kstart+1 , generic
- DO i = istart , iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ko_above_sfc(i) .EQ. -1 ) THEN
- IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
- ko_above_sfc(i) = ko
- END IF
- END IF
- END DO
- END DO
- ! Piece together columns of the original input data. Pass the vertical columns to
- ! the iterpolator.
- DO i = istart , iend
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! If the surface value is in the middle of the array, three steps: 1) do the
- ! values below the ground (this is just to catch the occasional value that is
- ! inconsistently below the surface based on input data), 2) do the surface level, then
- ! 3) add in the levels that are above the surface. For the levels next to the surface,
- ! we check to remove any levels that are "too close". When building the column of input
- ! pressures, we also attend to the request for forcing the surface analysis to be used
- ! in a few lower eta-levels.
- ! Fill in the column from up to the level just below the surface with the input
- ! presssure and the input field (orig or old, which ever). For an isobaric input
- ! file, this data is isobaric.
- ! How many levels have we skipped in the input column.
- zap = 0
- zap_below = 0
- zap_above = 0
- IF ( ko_above_sfc(i) .GT. 2 ) THEN
- count = 1
- DO ko = 2 , ko_above_sfc(i)-1
- ordered_porig(count) = porig(i,ko,j)
- ordered_forig(count) = forig(i,ko,j)
- count = count + 1
- END DO
- ! Make sure the pressure just below the surface is not "too close", this
- ! will cause havoc with the higher order interpolators. In case of a "too close"
- ! instance, we toss out the offending level (NOT the surface one) by simply
- ! decrementing the accumulating loop counter.
- IF ( ordered_porig(count-1) - porig(i,1,j) .LT. zap_close_levels ) THEN
- count = count -1
- zap = 1
- zap_below = 1
- END IF
- ! Add in the surface values.
- ordered_porig(count) = porig(i,1,j)
- ordered_forig(count) = forig(i,1,j)
- count = count + 1
- ! A usual way to do the vertical interpolation is to pay more attention to the
- ! surface data. Why? Well it has about 20x the density as the upper air, so we
- ! hope the analysis is better there. We more strongly use this data by artificially
- ! tossing out levels above the surface that are beneath a certain number of prescribed
- ! eta levels at this (i,j). The "zap" value is how many levels of input we are
- ! removing, which is used to tell the interpolator how many valid values are in
- ! the column. The "count" value is the increment to the index of levels, and is
- ! only used for assignments.
- IF ( force_sfc_in_vinterp .GT. 0 ) THEN
- ! Get the pressure at the eta level. We want to remove all input pressure levels
- ! between the level above the surface to the pressure at this eta surface. That
- ! forces the surface value to be used through the selected eta level. Keep track
- ! of two things: the level to use above the eta levels, and how many levels we are
- ! skipping.
- knext = ko_above_sfc(i)
- find_level : DO ko = ko_above_sfc(i) , generic
- IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
- knext = ko
- exit find_level
- ELSE
- zap = zap + 1
- zap_above = zap_above + 1
- END IF
- END DO find_level
- ! No request for special interpolation, so we just assign the next level to use
- ! above the surface as, ta da, the first level above the surface. I know, wow.
- ELSE
- knext = ko_above_sfc(i)
- END IF
- ! One more time, make sure the pressure just above the surface is not "too close", this
- ! will cause havoc with the higher order interpolators. In case of a "too close"
- ! instance, we toss out the offending level above the surface (NOT the surface one) by simply
- ! incrementing the loop counter. Here, count-1 is the surface level and knext is either
- ! the next level up OR it is the level above the prescribed number of eta surfaces.
- IF ( ordered_porig(count-1) - porig(i,knext,j) .LT. zap_close_levels ) THEN
- kst = knext+1
- zap = zap + 1
- zap_above = zap_above + 1
- ELSE
- kst = knext
- END IF
- DO ko = kst , generic
- ordered_porig(count) = porig(i,ko,j)
- ordered_forig(count) = forig(i,ko,j)
- count = count + 1
- END DO
- ! This is easy, the surface is the lowest level, just stick them in, in this order. OK,
- ! there are a couple of subtleties. We have to check for that special interpolation that
- ! skips some input levels so that the surface is used for the lowest few eta levels. Also,
- ! we must make sure that we still do not have levels that are "too close" together.
- ELSE
- ! Initialize no input levels have yet been removed from consideration.
- zap = 0
- ! The surface is the lowest level, so it gets set right away to location 1.
- ordered_porig(1) = porig(i,1,j)
- ordered_forig(1) = forig(i,1,j)
- ! We start filling in the array at loc 2, as in just above the level we just stored.
- count = 2
- ! Are we forcing the interpolator to skip valid input levels so that the
- ! surface data is used through more levels? Essentially as above.
- IF ( force_sfc_in_vinterp .GT. 0 ) THEN
- knext = 2
- find_level2: DO ko = 2 , generic
- IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
- knext = ko
- exit find_level2
- ELSE
- zap = zap + 1
- zap_above = zap_above + 1
- END IF
- END DO find_level2
- ELSE
- knext = 2
- END IF
- ! Fill in the data above the surface. The "knext" index is either the one
- ! just above the surface OR it is the index associated with the level that
- ! is just above the pressure at this (i,j) of the top eta level that is to
- ! be directly impacted with the surface level in interpolation.
- DO ko = knext , generic
- IF ( ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) .AND. &
- ( ko .LT. generic ) ) THEN
- zap = zap + 1
- zap_above = zap_above + 1
- CYCLE
- END IF
- ordered_porig(count) = porig(i,ko,j)
- ordered_forig(count) = forig(i,ko,j)
- count = count + 1
- END DO
- END IF
- ! Now get the column of the "new" pressure data. So, this one is easy.
- DO kn = kstart , kend
- ordered_pnew(kn) = pnew(i,kn,j)
- END DO
- ! How many levels (count) are we shipping to the Lagrange interpolator.
- IF ( ( use_levels_below_ground ) .AND. ( use_surface ) ) THEN
- ! Use all levels, including the input surface, and including the pressure
- ! levels below ground. We know to stop when we have reached the top of
- ! the input pressure data.
- count = 0
- find_how_many_1 : DO ko = 1 , generic
- IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
- count = count + 1
- EXIT find_how_many_1
- ELSE
- count = count + 1
- END IF
- END DO find_how_many_1
- kinterp_start = 1
- kinterp_end = kinterp_start + count - 1
- ELSE IF ( ( use_levels_below_ground ) .AND. ( .NOT. use_surface ) ) THEN
- ! Use all levels (excluding the input surface) and including the pressure
- ! levels below ground. We know to stop when we have reached the top of
- ! the input pressure data.
- count = 0
- find_sfc_2 : DO ko = 1 , generic
- IF ( porig(i,1,j) .EQ. ordered_porig(ko) ) THEN
- sfc_level = ko
- EXIT find_sfc_2
- END IF
- END DO find_sfc_2
- DO ko = sfc_level , generic-1
- ordered_porig(ko) = ordered_porig(ko+1)
- ordered_forig(ko) = ordered_forig(ko+1)
- END DO
- ordered_porig(generic) = 1.E-5
- ordered_forig(generic) = 1.E10
- count = 0
- find_how_many_2 : DO ko = 1 , generic
- IF ( porig(i,generic,j) .EQ. ordered_porig(ko) ) THEN
- count = count + 1
- EXIT find_how_many_2
- ELSE
- count = count + 1
- END IF
- END DO find_how_many_2
- kinterp_start = 1
- kinterp_end = kinterp_start + count - 1
- ELSE IF ( ( .NOT. use_levels_below_ground ) .AND. ( use_surface ) ) THEN
- ! Use all levels above the input surface pressure.
- kcount = ko_above_sfc(i)-1-zap_below
- count = 0
- DO ko = 1 , generic
- IF ( porig(i,ko,j) .EQ. ordered_porig(kcount) ) THEN
- ! write (6,fmt='(f11.3,f11.3,g11.5)') porig(i,ko,j),ordered_porig(kcount),ordered_forig(kcount)
- kcount = kcount + 1
- count = count + 1
- ELSE
- ! write (6,fmt='(f11.3 )') porig(i,ko,j)
- END IF
- END DO
- kinterp_start = ko_above_sfc(i)-1-zap_below
- kinterp_end = kinterp_start + count - 1
- END IF
- ! The polynomials are either in pressure or LOG(pressure).
- IF ( interp_type .EQ. 1 ) THEN
- CALL lagrange_setup ( var_type , interp_type , &
- ordered_porig(kinterp_start:kinterp_end) , &
- ordered_forig(kinterp_start:kinterp_end) , &
- count , lagrange_order , extrap_type , &
- ordered_pnew(kstart:kend) , ordered_fnew , kend-kstart+1 ,i,j)
- ELSE
- CALL lagrange_setup ( var_type , interp_type , &
- LOG(ordered_porig(kinterp_start:kinterp_end)) , &
- ordered_forig(kinterp_start:kinterp_end) , &
- count , lagrange_order , extrap_type , &
- LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j)
- END IF
- ! Save the computed data.
- DO kn = kstart , kend
- fnew(i,kn,j) = ordered_fnew(kn)
- END DO
- ! There may have been a request to have the surface data from the input field
- ! to be assigned as to the lowest eta level. This assumes thin layers (usually
- ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
- IF ( lowest_lev_from_sfc ) THEN
- fnew(i,1,j) = forig(i,ko_above_sfc(i)-1,j)
- END IF
- END DO
- END DO
- END SUBROUTINE vert_interp
- !---------------------------------------------------------------------
- SUBROUTINE vert_interp_old ( forig , po , fnew , pnu , &
- generic , var_type , &
- interp_type , lagrange_order , extrap_type , &
- lowest_lev_from_sfc , use_levels_below_ground , use_surface , &
- zap_close_levels , force_sfc_in_vinterp , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Vertically interpolate the new field. The original field on the original
- ! pressure levels is provided, and the new pressure surfaces to interpolate to.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: interp_type , lagrange_order , extrap_type
- LOGICAL , INTENT(IN) :: lowest_lev_from_sfc , use_levels_below_ground , use_surface
- REAL , INTENT(IN) :: zap_close_levels
- INTEGER , INTENT(IN) :: force_sfc_in_vinterp
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- INTEGER , INTENT(IN) :: generic
- CHARACTER (LEN=1) :: var_type
- REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: forig , po
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
- REAL , DIMENSION(ims:ime,generic,jms:jme) :: porig
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
- ! Local vars
- INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2
- INTEGER :: istart , iend , jstart , jend , kstart , kend
- INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
- INTEGER , DIMENSION(ims:ime ) :: ks
- INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
- LOGICAL :: any_below_ground
- REAL :: p1 , p2 , pn
- integer vert_extrap
- vert_extrap = 0
- ! Horiontal loop bounds for different variable types.
- IF ( var_type .EQ. 'U' ) THEN
- istart = its
- iend = ite
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = MAX(ids+1,its) , MIN(ide-1,ite)
- porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
- END DO
- END DO
- IF ( ids .EQ. its ) THEN
- DO k = 1,generic
- porig(its,k,j) = po(its,k,j)
- END DO
- END IF
- IF ( ide .EQ. ite ) THEN
- DO k = 1,generic
- porig(ite,k,j) = po(ite-1,k,j)
- END DO
- END IF
- DO k = kstart,kend
- DO i = MAX(ids+1,its) , MIN(ide-1,ite)
- pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
- END DO
- END DO
- IF ( ids .EQ. its ) THEN
- DO k = kstart,kend
- pnew(its,k,j) = pnu(its,k,j)
- END DO
- END IF
- IF ( ide .EQ. ite ) THEN
- DO k = kstart,kend
- pnew(ite,k,j) = pnu(ite-1,k,j)
- END DO
- END IF
- END DO
- ELSE IF ( var_type .EQ. 'V' ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = jte
- kstart = kts
- kend = kte-1
- DO i = istart,iend
- DO k = 1,generic
- DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
- porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
- END DO
- END DO
- IF ( jds .EQ. jts ) THEN
- DO k = 1,generic
- porig(i,k,jts) = po(i,k,jts)
- END DO
- END IF
- IF ( jde .EQ. jte ) THEN
- DO k = 1,generic
- porig(i,k,jte) = po(i,k,jte-1)
- END DO
- END IF
- DO k = kstart,kend
- DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
- pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
- END DO
- END DO
- IF ( jds .EQ. jts ) THEN
- DO k = kstart,kend
- pnew(i,k,jts) = pnu(i,k,jts)
- END DO
- END IF
- IF ( jde .EQ. jte ) THEN
- DO k = kstart,kend
- pnew(i,k,jte) = pnu(i,k,jte-1)
- END DO
- END IF
- END DO
- ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- ELSE
- istart = its
- iend = MIN(ide-1,ite)
- jstart = jts
- jend = MIN(jde-1,jte)
- kstart = kts
- kend = kte-1
- DO j = jstart,jend
- DO k = 1,generic
- DO i = istart,iend
- porig(i,k,j) = po(i,k,j)
- END DO
- END DO
- DO k = kstart,kend
- DO i = istart,iend
- pnew(i,k,j) = pnu(i,k,j)
- END DO
- END DO
- END DO
- END IF
- DO j = jstart , jend
- ! Skip all of the levels below ground in the original data based upon the surface pressure.
- ! The ko_above_sfc is the index in the pressure array that is above the surface. If there
- ! are no levels underground, this is index = 2. The remaining levels are eligible for use
- ! in the vertical interpolation.
- DO i = istart , iend
- ko_above_sfc(i) = -1
- END DO
- DO ko = kstart+1 , kend
- DO i = istart , iend
- IF ( ko_above_sfc(i) .EQ. -1 ) THEN
- IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
- ko_above_sfc(i) = ko
- END IF
- END IF
- END DO
- END DO
- ! Initialize interpolation location. These are the levels in the original pressure
- ! data that are physically below and above the targeted new pressure level.
- DO kn = kts , kte
- DO i = its , ite
- k_above(i,kn) = -1
- k_below(i,kn) = -2
- END DO
- END DO
- ! Starting location is no lower than previous found location. This is for O(n logn)
- ! and not O(n^2), where n is the number of vertical levels to search.
- DO i = its , ite
- ks(i) = 1
- END DO
- ! Find trapping layer for interpolation. The kn index runs through all of the "new"
- ! levels of data.
- DO kn = kstart , kend
- DO i = istart , iend
- ! For each "new" level (kn), we search to find the trapping levels in the "orig"
- ! data. Most of the time, the "new" levels are the eta surfaces, and the "orig"
- ! levels are the input pressure levels.
- found_trap_above : DO ko = ks(i) , generic-1
- ! Because we can have levels in the interpolation that are not valid,
- ! let's toss out any candidate orig pressure values that are below ground
- ! based on the surface pressure. If the level =1, then this IS the surface
- ! level, so we HAVE to keep that one, but maybe not the ones above. If the
- ! level (ks) is NOT=1, then we have to just CYCLE our loop to find a legit
- ! below-pressure value. If we are not below ground, then we choose two
- ! neighboring levels to test whether they surround the new pressure level.
- ! The input trapping levels that we are trying is the surface and the first valid
- ! level above the surface.
- IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .EQ. 1 ) ) THEN
- ko_1 = ko
- ko_2 = ko_above_sfc(i)
- ! The "below" level is underground, cycle until we get to a valid pressure
- ! above ground.
- ELSE IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .NE. 1 ) ) THEN
- CYCLE found_trap_above
- ! The "below" level is above the surface, so we are in the clear to test these
- ! two levels out.
- ELSE
- ko_1 = ko
- ko_2 = ko+1
- END IF
- ! The test of the candidate levels: "below" has to have a larger pressure, and
- ! "above" has to have a smaller pressure.
- ! OK, we found the correct two surrounding levels. The locations are saved for use in the
- ! interpolation.
- IF ( ( porig(i,ko_1,j) .GE. pnew(i,kn,j) ) .AND. &
- ( porig(i,ko_2,j) .LT. pnew(i,kn,j) ) ) THEN
- k_above(i,kn) = ko_2
- k_below(i,kn) = ko_1
- ks(i) = ko_1
- EXIT found_trap_above
- ! What do we do is we need to extrapolate the data underground? This happens when the
- ! lowest pressure that we have is physically "above" the new target pressure. Our
- ! actions depend on the type of variable we are interpolating.
- ELSE IF ( porig(i,1,j) .LT. pnew(i,kn,j) ) THEN
- ! For horizontal winds and moisture, we keep a constant value under ground.
- IF ( ( var_type .EQ. 'U' ) .OR. &
- ( var_type .EQ. 'V' ) .OR. &
- ( var_type .EQ. 'Q' ) ) THEN
- k_above(i,kn) = 1
- ks(i) = 1
- ! For temperature and height, we extrapolate the data. Hopefully, we are not
- ! extrapolating too far. For pressure level input, the eta levels are always
- ! contained within the surface to p_top levels, so no extrapolation is ever
- ! required.
- ELSE IF ( ( var_type .EQ. 'Z' ) .OR. &
- ( var_type .EQ. 'T' ) ) THEN
- k_above(i,kn) = ko_above_sfc(i)
- k_below(i,kn) = 1
- ks(i) = 1
- ! Just a catch all right now.
- ELSE
- k_above(i,kn) = 1
- ks(i) = 1
- END IF
- EXIT found_trap_above
- ! The other extrapolation that might be required is when we are going above the
- ! top level of the input data. Usually this means we chose a P_PTOP value that
- ! was inappropriate, and we should stop and let someone fix this mess.
- ELSE IF ( porig(i,generic,j) .GT. pnew(i,kn,j) ) THEN
- print *,'data is too high, try a lower p_top'
- print *,'pnew=',pnew(i,kn,j)
- print *,'porig=',porig(i,:,j)
- CALL wrf_error_fatal ('requested p_top is higher than input data, lower p_top')
- END IF
- END DO found_trap_above
- END DO
- END DO
- ! Linear vertical interpolation.
- DO kn = kstart , kend
- DO i = istart , iend
- IF ( k_above(i,kn) .EQ. 1 ) THEN
- fnew(i,kn,j) = forig(i,1,j)
- ELSE
- k2 = MAX ( k_above(i,kn) , 2)
- k1 = MAX ( k_below(i,kn) , 1)
- IF ( k1 .EQ. k2 ) THEN
- CALL wrf_error_fatal ( 'identical values in the interp, bad for divisions' )
- END IF
- IF ( interp_type .EQ. 1 ) THEN
- p1 = porig(i,k1,j)
- p2 = porig(i,k2,j)
- pn = pnew(i,kn,j)
- ELSE IF ( interp_type .EQ. 2 ) THEN
- p1 = ALOG(porig(i,k1,j))
- p2 = ALOG(porig(i,k2,j))
- pn = ALOG(pnew(i,kn,j))
- END IF
- IF ( ( p1-pn) * (p2-pn) > 0. ) THEN
- ! CALL wrf_error_fatal ( 'both trapping pressures are on the same side of the new pressure' )
- ! CALL wrf_debug ( 0 , 'both trapping pressures are on the same side of the new pressure' )
- vert_extrap = vert_extrap + 1
- END IF
- fnew(i,kn,j) = ( forig(i,k1,j) * ( p2 - pn ) + &
- forig(i,k2,j) * ( pn - p1 ) ) / &
- ( p2 - p1 )
- END IF
- END DO
- END DO
- search_below_ground : DO kn = kstart , kend
- any_below_ground = .FALSE.
- DO i = istart , iend
- IF ( k_above(i,kn) .EQ. 1 ) THEN
- fnew(i,kn,j) = forig(i,1,j)
- any_below_ground = .TRUE.
- END IF
- END DO
- IF ( .NOT. any_below_ground ) THEN
- EXIT search_below_ground
- END IF
- END DO search_below_ground
- ! There may have been a request to have the surface data from the input field
- ! to be assigned as to the lowest eta level. This assumes thin layers (usually
- ! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
- DO i = istart , iend
- IF ( lowest_lev_from_sfc ) THEN
- fnew(i,1,j) = forig(i,ko_above_sfc(i),j)
- END IF
- END DO
- END DO
- print *,'VERT EXTRAP = ', vert_extrap
- END SUBROUTINE vert_interp_old
- !---------------------------------------------------------------------
- SUBROUTINE lagrange_setup ( var_type , interp_type , all_x , all_y , all_dim , n , extrap_type , &
- target_x , target_y , target_dim ,i,j)
- ! We call a Lagrange polynomial interpolator. The parallel concerns are put off as this
- ! is initially set up for vertical use. The purpose is an input column of pressure (all_x),
- ! and the associated pressure level data (all_y). These are assumed to be sorted (ascending
- ! or descending, no matter). The locations to be interpolated to are the pressures in
- ! target_x, probably the new vertical coordinate values. The field that is output is the
- ! target_y, which is defined at the target_x location. Mostly we expect to be 2nd order
- ! overlapping polynomials, with only a single 2nd order method near the top and bottom.
- ! When n=1, this is linear; when n=2, this is a second order interpolator.
- IMPLICIT NONE
- CHARACTER (LEN=1) :: var_type
- INTEGER , INTENT(IN) :: interp_type , all_dim , n , extrap_type , target_dim
- REAL, DIMENSION(all_dim) , INTENT(IN) :: all_x , all_y
- REAL , DIMENSION(target_dim) , INTENT(IN) :: target_x
- REAL , DIMENSION(target_dim) , INTENT(OUT) :: target_y
- ! Brought in for debug purposes, all of the computations are in a single column.
- INTEGER , INTENT(IN) :: i,j
- ! Local vars
- REAL , DIMENSION(n+1) :: x , y
- REAL :: a , b
- REAL :: target_y_1 , target_y_2
- LOGICAL :: found_loc
- INTEGER :: loop , loc_center_left , loc_center_right , ist , iend , target_loop
- INTEGER :: vboundb , vboundt
- ! Local vars for the problem of extrapolating theta below ground.
- REAL :: temp_1 , temp_2 , temp_3 , temp_y
- REAL :: depth_of_extrap_in_p , avg_of_extrap_p , temp_extrap_starting_point , dhdp , dh , dt
- REAL , PARAMETER :: RovCp = rcp
- REAL , PARAMETER :: CRC_const1 = 11880.516 ! m
- REAL , PARAMETER :: CRC_const2 = 0.1902632 !
- REAL , PARAMETER :: CRC_const3 = 0.0065 ! K/km
- REAL, DIMENSION(all_dim) :: all_x_full
- REAL , DIMENSION(target_dim) :: target_x_full
- IF ( all_dim .LT. n+1 ) THEN
- print *,'all_dim = ',all_dim
- print *,'order = ',n
- print *,'i,j = ',i,j
- print *,'p array = ',all_x
- print *,'f array = ',all_y
- print *,'p target= ',target_x
- CALL wrf_error_fatal ( 'troubles, the interpolating order is too large for this few input values' )
- END IF
- IF ( n .LT. 1 ) THEN
- CALL wrf_error_fatal ( 'pal, linear is about as low as we go' )
- END IF
- ! We can pinch in the area of the higher order interpolation with vbound. If
- ! vbound = 0, no pinching. If vbound = m, then we make the lower "m" and upper
- ! "m" eta levels use a linear interpolation.
- vboundb = 4
- vboundt = 0
- ! Loop over the list of target x and y values.
- DO target_loop = 1 , target_dim
- ! Find the two trapping x values, and keep the indices.
- found_loc = .FALSE.
- find_trap : DO loop = 1 , all_dim -1
- a = target_x(target_loop) - all_x(loop)
- b = target_x(target_loop) - all_x(loop+1)
- IF ( a*b .LE. 0.0 ) THEN
- loc_center_left = loop
- loc_center_right = loop+1
- found_loc = .TRUE.
- EXIT find_trap
- END IF
- END DO find_trap
- IF ( ( .NOT. found_loc ) .AND. ( target_x(target_loop) .GT. all_x(1) ) ) THEN
- ! Get full pressure back so that our extrpolations make sense.
- IF ( interp_type .EQ. 1 ) THEN
- all_x_full = all_x
- target_x_full = target_x
- ELSE
- all_x_full = EXP ( all_x )
- target_x_full = EXP ( target_x )
- END IF
- ! Isothermal extrapolation.
- IF ( ( extrap_type .EQ. 1 ) .AND. ( var_type .EQ. 'T' ) ) THEN
- temp_1 = all_y(1) * ( all_x_full(1) / 100000. ) ** RovCp
- target_y(target_loop) = temp_1 * ( 100000. / target_x_full(target_loop) ) ** RovCp
- ! Standard atmosphere -6.5 K/km lapse rate for the extrapolation.
- ELSE IF ( ( extrap_type .EQ. 2 ) .AND. ( var_type .EQ. 'T' ) ) THEN
- depth_of_extrap_in_p = target_x_full(target_loop) - all_x_full(1)
- avg_of_extrap_p = ( target_x_full(target_loop) + all_x_full(1) ) * 0.5
- temp_extrap_starting_point = all_y(1) * ( all_x_full(1) / 100000. ) ** RovCp
- dhdp = CRC_const1 * CRC_const2 * ( avg_of_extrap_p / 100. ) ** ( CRC_const2 - 1. )
- dh = dhdp * ( depth_of_extrap_in_p / 100. )
- dt = dh * CRC_const3
- target_y(target_loop) = ( temp_extrap_starting_point + dt ) * ( 100000. / target_x_full(target_loop) ) ** RovCp
- ! Adiabatic extrapolation for theta.
- ELSE IF ( ( extrap_type .EQ. 3 ) .AND. ( var_type .EQ. 'T' ) ) THEN
- target_y(target_loop) = all_y(1)
- ! Wild extrapolation for non-temperature vars.
- ELSE IF ( extrap_type .EQ. 1 ) THEN
- target_y(target_loop) = ( all_y(2) * ( target_x(target_loop) - all_x(3) ) + &
- all_y(3) * ( all_x(2) - target_x(target_loop) ) ) / &
- ( all_x(2) - all_x(3) )
- ! Use a constant value below ground.
- ELSE IF ( extrap_type .EQ. 2 ) THEN
- target_y(target_loop) = all_y(1)
- ELSE IF ( extrap_type .EQ. 3 ) THEN
- CALL wrf_error_fatal ( 'You are not allowed to use extrap_option #3 for any var except for theta.' )
- END IF
- CYCLE
- ELSE IF ( .NOT. found_loc ) THEN
- print *,'i,j = ',i,j
- print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
- DO loop = 1 , all_dim
- print *,'column of pressure and value = ',all_x(loop),all_y(loop)
- END DO
- CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
- END IF
- ! Even or odd order? We can put the value in the middle if this is
- ! an odd order interpolator. For the even guys, we'll do it twice
- ! and shift the range one index, then get an average.
- IF ( MOD(n,2) .NE. 0 ) THEN
- IF ( ( loc_center_left -(((n+1)/2)-1) .GE. 1 ) .AND. &
- ( loc_center_right+(((n+1)/2)-1) .LE. all_dim ) ) THEN
- ist = loc_center_left -(((n+1)/2)-1)
- iend = ist + n
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
- ELSE
- IF ( .NOT. found_loc ) THEN
- CALL wrf_error_fatal ( 'I doubt this will happen, I will only do 2nd order for now' )
- END IF
- END IF
- ELSE IF ( ( MOD(n,2) .EQ. 0 ) .AND. &
- ( ( target_loop .GE. 1 + vboundb ) .AND. ( target_loop .LE. target_dim - vboundt ) ) ) THEN
- IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
- ( loc_center_right+(((n )/2) ) .LE. all_dim ) .AND. &
- ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
- ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
- ist = loc_center_left -(((n )/2)-1)
- iend = ist + n
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_1 )
- ist = loc_center_left -(((n )/2) )
- iend = ist + n
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_2 )
- target_y(target_loop) = ( target_y_1 + target_y_2 ) * 0.5
- ELSE IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
- ( loc_center_right+(((n )/2) ) .LE. all_dim ) ) THEN
- ist = loc_center_left -(((n )/2)-1)
- iend = ist + n
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
- ELSE IF ( ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
- ( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
- ist = loc_center_left -(((n )/2) )
- iend = ist + n
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
- ELSE
- CALL wrf_error_fatal ( 'unauthorized area, you should not be here' )
- END IF
- ELSE IF ( MOD(n,2) .EQ. 0 ) THEN
- ist = loc_center_left
- iend = loc_center_right
- CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , 1 , target_x(target_loop) , target_y(target_loop) )
- END IF
- END DO
- END SUBROUTINE lagrange_setup
- !---------------------------------------------------------------------
- SUBROUTINE lagrange_interp ( x , y , n , target_x , target_y )
- ! Interpolation using Lagrange polynomials.
- ! P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
- ! where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
- ! ---------------------------------------------
- ! (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: n
- REAL , DIMENSION(0:n) , INTENT(IN) :: x , y
- REAL , INTENT(IN) :: target_x
- REAL , INTENT(OUT) :: target_y
- ! Local vars
- INTEGER :: i , k
- REAL :: numer , denom , Px
- REAL , DIMENSION(0:n) :: Ln
- Px = 0.
- DO i = 0 , n
- numer = 1.
- denom = 1.
- DO k = 0 , n
- IF ( k .EQ. i ) CYCLE
- numer = numer * ( target_x - x(k) )
- denom = denom * ( x(i) - x(k) )
- END DO
- IF ( denom .NE. 0. ) THEN
- Ln(i) = y(i) * numer / denom
- Px = Px + Ln(i)
- ENDIF
- END DO
- target_y = Px
- END SUBROUTINE lagrange_interp
- #ifndef VERT_UNIT
- !---------------------------------------------------------------------
- SUBROUTINE p_dry ( mu0 , eta , pdht , pdry , full_levs , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute reference pressure and the reference mu.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- LOGICAL :: full_levs
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: mu0
- REAL , DIMENSION( kms:kme ) , INTENT(IN) :: eta
- REAL :: pdht
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pdry
- ! Local vars
- INTEGER :: i , j , k
- REAL , DIMENSION( kms:kme ) :: eta_h
- IF ( full_levs ) THEN
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pdry(i,k,j) = eta(k) * mu0(i,j) + pdht
- END DO
- END DO
- END DO
- ELSE
- DO k = kts , kte-1
- eta_h(k) = ( eta(k) + eta(k+1) ) * 0.5
- END DO
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pdry(i,k,j) = eta_h(k) * mu0(i,j) + pdht
- END DO
- END DO
- END DO
- END IF
- END SUBROUTINE p_dry
- !---------------------------------------------------------------------
- SUBROUTINE p_dts ( pdts , intq , psfc , p_top , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute difference between the dry, total surface pressure and the top pressure.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , INTENT(IN) :: p_top
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: psfc
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: intq
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: pdts
- ! Local vars
- INTEGER :: i , j , k
- DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pdts(i,j) = psfc(i,j) - intq(i,j) - p_top
- END DO
- END DO
- END SUBROUTINE p_dts
- !---------------------------------------------------------------------
- SUBROUTINE p_dhs ( pdhs , ht , p0 , t0 , a , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute dry, hydrostatic surface pressure.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: ht
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: pdhs
- REAL , INTENT(IN) :: p0 , t0 , a
- ! Local vars
- INTEGER :: i , j , k
- REAL , PARAMETER :: Rd = r_d
- DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pdhs(i,j) = p0 * EXP ( -t0/a + SQRT ( (t0/a)**2 - 2. * g * ht(i,j)/(a * Rd) ) )
- END DO
- END DO
- END SUBROUTINE p_dhs
- !---------------------------------------------------------------------
- SUBROUTINE find_p_top ( p , p_top , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Find the largest pressure in the top level. This is our p_top. We are
- ! assuming that the top level is the location where the pressure is a minimum
- ! for each column. In cases where the top surface is not isobaric, a
- ! communicated value must be shared in the calling routine. Also in cases
- ! where the top surface is not isobaric, care must be taken that the new
- ! maximum pressure is not greater than the previous value. This test is
- ! also handled in the calling routine.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL :: p_top
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
- ! Local vars
- INTEGER :: i , j , k, min_lev
- i = its
- j = jts
- p_top = p(i,2,j)
- min_lev = 2
- DO k = 2 , kte
- IF ( p_top .GT. p(i,k,j) ) THEN
- p_top = p(i,k,j)
- min_lev = k
- END IF
- END DO
- k = min_lev
- p_top = p(its,k,jts)
- DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- p_top = MAX ( p_top , p(i,k,j) )
- END DO
- END DO
- END SUBROUTINE find_p_top
- !---------------------------------------------------------------------
- SUBROUTINE t_to_theta ( t , p , p00 , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute potential temperature from temperature and pressure.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , INTENT(IN) :: p00
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
- ! Local vars
- INTEGER :: i , j , k
- REAL , PARAMETER :: Rd = r_d
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- t(i,k,j) = t(i,k,j) * ( p00 / p(i,k,j) ) ** (Rd / Cp)
- END DO
- END DO
- END DO
- END SUBROUTINE t_to_theta
- !---------------------------------------------------------------------
- SUBROUTINE theta_to_t ( t , p , p00 , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute temperature from potential temp and pressure.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , INTENT(IN) :: p00
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
- ! Local vars
- INTEGER :: i , j , k
- REAL , PARAMETER :: Rd = r_d
- CHARACTER (LEN=80) :: mess
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- if ( p(i,k,j) .NE. 0. ) then
- t(i,k,j) = t(i,k,j) / ( ( p00 / p(i,k,j) ) ** (Rd / Cp) )
- else
- WRITE(mess,*) 'Troubles in theta_to_t'
- CALL wrf_debug(0,mess)
- WRITE(mess,*) "i,j,k = ", i,j,k
- CALL wrf_debug(0,mess)
- WRITE(mess,*) "p(i,k,j) = ", p(i,k,j)
- CALL wrf_debug(0,mess)
- WRITE(mess,*) "t(i,k,j) = ", t(i,k,j)
- CALL wrf_debug(0,mess)
- endif
- END DO
- END DO
- END DO
- END SUBROUTINE theta_to_t
- !---------------------------------------------------------------------
- SUBROUTINE integ_moist ( q_in , p_in , pd_out , t_in , ght_in , intq , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Integrate the moisture field vertically. Mostly used to get the total
- ! vapor pressure, which can be subtracted from the total pressure to get
- ! the dry pressure.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: q_in , p_in , t_in , ght_in
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pd_out
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: intq
- ! Local vars
- INTEGER :: i , j , k
- INTEGER , DIMENSION(ims:ime) :: level_above_sfc
- REAL , DIMENSION(ims:ime,jms:jme) :: psfc , tsfc , qsfc, zsfc
- REAL , DIMENSION(ims:ime,kms:kme) :: q , p , t , ght, pd
- REAL :: rhobar , qbar , dz
- REAL :: p1 , p2 , t1 , t2 , q1 , q2 , z1, z2
- LOGICAL :: upside_down
- LOGICAL :: already_assigned_upside_down
- REAL , PARAMETER :: Rd = r_d
- ! Is the data upside down?
- already_assigned_upside_down = .FALSE.
- find_valid : DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( p_in(i,kts+1,j) .LT. p_in(i,kte,j) ) THEN
- upside_down = .TRUE.
- already_assigned_upside_down = .TRUE.
- ELSE
- upside_down = .FALSE.
- already_assigned_upside_down = .TRUE.
- END IF
- EXIT find_valid
- END DO
- END DO find_valid
- IF ( .NOT. already_assigned_upside_down ) THEN
- upside_down = .FALSE.
- END IF
- ! Get a surface value, always the first level of a 3d field.
- DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- psfc(i,j) = p_in(i,kts,j)
- tsfc(i,j) = t_in(i,kts,j)
- qsfc(i,j) = q_in(i,kts,j)
- zsfc(i,j) = ght_in(i,kts,j)
- END DO
- END DO
- DO j = jts , MIN ( jde-1 , jte )
- ! Initialize the integrated quantity of moisture to zero.
- DO i = its , MIN (ide-1 , ite )
- intq(i,j) = 0.
- END DO
- IF ( upside_down ) THEN
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- p(i,kts) = p_in(i,kts,j)
- t(i,kts) = t_in(i,kts,j)
- q(i,kts) = q_in(i,kts,j)
- ght(i,kts) = ght_in(i,kts,j)
- DO k = kts+1,kte
- p(i,k) = p_in(i,kte+2-k,j)
- t(i,k) = t_in(i,kte+2-k,j)
- q(i,k) = q_in(i,kte+2-k,j)
- ght(i,k) = ght_in(i,kte+2-k,j)
- END DO
- END DO
- ELSE
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- DO k = kts,kte
- p(i,k) = p_in(i,k ,j)
- t(i,k) = t_in(i,k ,j)
- q(i,k) = q_in(i,k ,j)
- ght(i,k) = ght_in(i,k ,j)
- END DO
- END DO
- END IF
- ! Find the first level above the ground. If all of the levels are above ground, such as
- ! a terrain following lower coordinate, then the first level above ground is index #2.
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- level_above_sfc(i) = -1
- IF ( p(i,kts+1) .LT. psfc(i,j) ) THEN
- level_above_sfc(i) = kts+1
- ELSE
- find_k : DO k = kts+1,kte-1
- IF ( ( p(i,k )-psfc(i,j) .GE. 0. ) .AND. &
- ( p(i,k+1)-psfc(i,j) .LT. 0. ) ) THEN
- level_above_sfc(i) = k+1
- EXIT find_k
- END IF
- END DO find_k
- IF ( level_above_sfc(i) .EQ. -1 ) THEN
- print *,'i,j = ',i,j
- print *,'p = ',p(i,:)
- print *,'p sfc = ',psfc(i,j)
- CALL wrf_error_fatal ( 'Could not find level above ground')
- END IF
- END IF
- END DO
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! Account for the moisture above the ground.
- pd(i,kte) = p(i,kte)
- DO k = kte-1,level_above_sfc(i),-1
- rhobar = ( p(i,k ) / ( Rd * t(i,k ) ) + &
- p(i,k+1) / ( Rd * t(i,k+1) ) ) * 0.5
- qbar = ( q(i,k ) + q(i,k+1) ) * 0.5
- dz = ght(i,k+1) - ght(i,k)
- intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
- pd(i,k) = p(i,k) - intq(i,j)
- END DO
- ! Account for the moisture between the surface and the first level up.
- IF ( ( p(i,level_above_sfc(i)-1)-psfc(i,j) .GE. 0. ) .AND. &
- ( p(i,level_above_sfc(i) )-psfc(i,j) .LT. 0. ) .AND. &
- ( level_above_sfc(i) .GT. kts ) ) THEN
- p1 = psfc(i,j)
- p2 = p(i,level_above_sfc(i))
- t1 = tsfc(i,j)
- t2 = t(i,level_above_sfc(i))
- q1 = qsfc(i,j)
- q2 = q(i,level_above_sfc(i))
- z1 = zsfc(i,j)
- z2 = ght(i,level_above_sfc(i))
- rhobar = ( p1 / ( Rd * t1 ) + &
- p2 / ( Rd * t2 ) ) * 0.5
- qbar = ( q1 + q2 ) * 0.5
- dz = z2 - z1
- IF ( dz .GT. 0.1 ) THEN
- intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
- END IF
- ! Fix the underground values.
- DO k = level_above_sfc(i)-1,kts+1,-1
- pd(i,k) = p(i,k) - intq(i,j)
- END DO
- END IF
- pd(i,kts) = psfc(i,j) - intq(i,j)
- END DO
- IF ( upside_down ) THEN
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- pd_out(i,kts,j) = pd(i,kts)
- DO k = kts+1,kte
- pd_out(i,kte+2-k,j) = pd(i,k)
- END DO
- END DO
- ELSE
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- DO k = kts,kte
- pd_out(i,k,j) = pd(i,k)
- END DO
- END DO
- END IF
- END DO
- END SUBROUTINE integ_moist
- !---------------------------------------------------------------------
- SUBROUTINE rh_to_mxrat2(rh, t, p, q , wrt_liquid , &
- qv_max_p_safe , &
- qv_max_flag , qv_max_value , &
- qv_min_p_safe , &
- qv_min_flag , qv_min_value , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! This subroutine computes mixing ratio (q, kg/kg) from basic variables
- ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 0-100%).
- ! Phase transition, liquid water to ice, occurs over (0,-23) temperature range (Celcius).
- ! Formulation used here is based on:
- ! WMO, General meteorological standards and recommended practices,
- ! Appendix A, WMO Technical Regulations, WMO-No. 49, corrigendum,
- ! August 2000. --TKW 03/30/2011
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- LOGICAL , INTENT(IN) :: wrt_liquid
- REAL , INTENT(IN) :: qv_max_p_safe , qv_max_flag , qv_max_value
- REAL , INTENT(IN) :: qv_min_p_safe , qv_min_flag , qv_min_value
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
- ! Local vars
- REAL, PARAMETER :: T0K = 273.16
- REAL, PARAMETER :: Tice = T0K - 23.0
- REAL, PARAMETER :: cfe = 1.0/(23.0*23.0)
- REAL, PARAMETER :: eps = 0.622
- ! Coefficients for esat over liquid water
- REAL, PARAMETER :: cw1 = 10.79574
- REAL, PARAMETER :: cw2 = -5.02800
- REAL, PARAMETER :: cw3 = 1.50475E-4
- REAL, PARAMETER :: cw4 = 0.42873E-3
- REAL, PARAMETER :: cw5 = 0.78614
- ! Coefficients for esat over ice
- REAL, PARAMETER :: ci1 = -9.09685
- REAL, PARAMETER :: ci2 = -3.56654
- REAL, PARAMETER :: ci3 = 0.87682
- REAL, PARAMETER :: ci4 = 0.78614
- REAL, PARAMETER :: Tn = 273.16
- ! 1 ppm is a reasonable estimate for minimum QV even for stratospheric altitudes
- REAL, PARAMETER :: QV_MIN = 1.e-6
- ! Maximum allowed QV is computed under the extreme condition:
- ! Saturated at 40 degree in Celcius and 1000 hPa
- REAL, PARAMETER :: QV_MAX = 0.045
- ! Need to constrain WVP in the stratosphere where pressure
- ! is low but tempearure is hot (warm)
- ! Maximum ratio of e/p, = q/(0.622+q)
- REAL, PARAMETER :: EP_MAX = QV_MAX/(eps+QV_MAX)
- INTEGER :: i , j , k
- REAL :: ew , q1 , t1
- REAL :: ta, tb, pw3, pw4, pwr
- REAL :: es, esw, esi, wvp, pmb, wvpmax
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 0. ) , 100. )
- END DO
- END DO
- END DO
- IF ( wrt_liquid ) THEN
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- Ta=Tn/T(i,k,j)
- Tb=T(i,k,j)/Tn
- pw3 = -8.2969*(Tb-1.0)
- pw4 = 4.76955*(1.0-Ta)
- pwr = cw1*(1.0-Ta) + cw2*LOG10(Tb) + cw3*(1.0-10.0**pw3) + cw4*(10.0**pw4-1.0) + cw5
- es = 10.0**pwr ! Saturation WVP
- wvp = 0.01*rh(i,k,j)*es ! Actual WVP
- pmb = p(i,k,j)/100.
- wvpmax = EP_MAX*pmb ! Prevents unrealistic QV in the stratosphere
- wvp = MIN(wvp,wvpmax)
- q(i,k,j) = eps*wvp/(pmb-wvp)
- END DO
- END DO
- END DO
- ELSE
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- Ta=Tn/T(i,k,j)
- Tb=T(i,k,j)/Tn
- IF (t(i,k,j) >= T0K) THEN ! Over liquid water
- pw3 = -8.2969*(Tb-1.0)
- pw4 = 4.76955*(1.0-Ta)
- pwr = cw1*(1.0-Ta) + cw2*LOG10(Tb) + cw3*(1.0-10.0**pw3) + cw4*(10.0**pw4-1.0) + cw5
- es = 10.0**pwr
- wvp = 0.01*rh(i,k,j)*es
- ELSE IF (t(i,k,j) <= Tice) THEN ! Over ice
- pwr = ci1*(Ta-1.0) + ci2*LOG10(Ta) + ci3*(1.0-Tb) + ci4
- es = 10.0**pwr
- wvp = 0.01*rh(i,k,j)*es
- ELSE ! Mixed
- pw3 = -8.2969*(Tb-1.0)
- pw4 = 4.76955*(1.0-Ta)
- pwr = cw1*(1.0-Ta) + cw2*LOG10(Tb) + cw3*(1.0-10.0**pw3) + cw4*(10.0**pw4-1.0) + cw5
- esw = 10.0**pwr ! Over liquid water
- pwr = ci1*(Ta-1.0) + ci2*LOG10(Ta) + ci3*(1.0-Tb) + ci4
- esi = 10.0**pwr ! Over ice
- es = esi + (esw-esi)*cfe*(T(i,k,j)-Tice)*(T(i,k,j)-Tice)
- wvp = 0.01*rh(i,k,j)*es
- END IF
- pmb = p(i,k,j)/100.
- wvpmax = EP_MAX*pmb ! Prevents unrealistic QV in the stratosphere
- wvp = MIN(wvp,wvpmax)
- q(i,k,j) = eps*wvp/(pmb-wvp)
- END DO
- END DO
- END DO
- END IF
- ! For pressures above a defined level, reasonable Qv values should be
- ! a certain value or smaller. If they are larger than this, the input data
- ! probably had "missing" RH, and we filled in some values. This is an
- ! attempt to catch those. Also, set the minimum value for the entire
- ! domain that is above the selected pressure level.
-
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( p(i,k,j) .LT. qv_max_p_safe ) THEN
- IF ( q(i,k,j) .GT. qv_max_flag ) THEN
- q(i,k,j) = qv_max_value
- END IF
- END IF
- IF ( p(i,k,j) .LT. qv_min_p_safe ) THEN
- IF ( q(i,k,j) .LT. qv_min_flag ) THEN
- q(i,k,j) = qv_min_value
- END IF
- END IF
- END DO
- END DO
- END DO
- END SUBROUTINE rh_to_mxrat2
- !---------------------------------------------------------------------
- SUBROUTINE rh_to_mxrat1(rh, t, p, q , wrt_liquid , &
- qv_max_p_safe , &
- qv_max_flag , qv_max_value , &
- qv_min_p_safe , &
- qv_min_flag , qv_min_value , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- LOGICAL , INTENT(IN) :: wrt_liquid
- REAL , INTENT(IN) :: qv_max_p_safe , qv_max_flag , qv_max_value
- REAL , INTENT(IN) :: qv_min_p_safe , qv_min_flag , qv_min_value
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
- REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
- ! Local vars
- INTEGER :: i , j , k
- REAL :: ew , q1 , t1
- REAL, PARAMETER :: T_REF = 0.0
- REAL, PARAMETER :: MW_AIR = 28.966
- REAL, PARAMETER :: MW_VAP = 18.0152
- REAL, PARAMETER :: A0 = 6.107799961
- REAL, PARAMETER :: A1 = 4.436518521e-01
- REAL, PARAMETER :: A2 = 1.428945805e-02
- REAL, PARAMETER :: A3 = 2.650648471e-04
- REAL, PARAMETER :: A4 = 3.031240396e-06
- REAL, PARAMETER :: A5 = 2.034080948e-08
- REAL, PARAMETER :: A6 = 6.136820929e-11
- REAL, PARAMETER :: ES0 = 6.1121
- REAL, PARAMETER :: C1 = 9.09718
- REAL, PARAMETER :: C2 = 3.56654
- REAL, PARAMETER :: C3 = 0.876793
- REAL, PARAMETER :: EIS = 6.1071
- REAL :: RHS
- REAL, PARAMETER :: TF = 273.16
- REAL :: TK
- REAL :: ES
- REAL :: QS
- REAL, PARAMETER :: EPS = 0.622
- REAL, PARAMETER :: SVP1 = 0.6112
- REAL, PARAMETER :: SVP2 = 17.67
- REAL, PARAMETER :: SVP3 = 29.65
- REAL, PARAMETER :: SVPT0 = 273.15
- CHARACTER (LEN=80) :: mess
- ! This subroutine computes mixing ratio (q, kg/kg) from basic variables
- ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
- ! The reference temperature (t_ref, C) is used to describe the temperature
- ! at which the liquid and ice phase change occurs.
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 0. ) , 100. )
- END DO
- END DO
- END DO
- IF ( wrt_liquid ) THEN
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! es is reduced by RH here to avoid problems in low-pressure cases
- if (t(i,k,j) .ne. 0.) then
- es=.01*rh(i,k,j)*svp1*10.*EXP(svp2*(t(i,k,j)-svpt0)/(t(i,k,j)-svp3))
- IF (es .ge. p(i,k,j)/100.)THEN
- q(i,k,j)=0.0
- WRITE(mess,*) 'Warning: vapor pressure exceeds total pressure, setting Qv to 0'
- CALL wrf_debug(0,mess)
- ELSE
- q(i,k,j)=eps*es/(p(i,k,j)/100.-es)
- ENDIF
- else
- q(i,k,j)=0.0
- WRITE(mess,*) 't(i,j,k) was 0 at ', i,j,k,', setting Qv to 0'
- CALL wrf_debug(0,mess)
- endif
- END DO
- END DO
- END DO
- ELSE
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- t1 = t(i,k,j) - 273.16
- ! Obviously dry.
- IF ( t1 .lt. -200. ) THEN
- q(i,k,j) = 0
- ELSE
- ! First compute the ambient vapor pressure of water
- ! Liquid phase t > 0 C
- IF ( t1 .GE. t_ref ) THEN
- ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
- ! Mixed phase -47 C < t < 0 C
- ELSE IF ( ( t1 .LT. t_ref ) .AND. ( t1 .GE. -47. ) ) THEN
- ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
- ! Ice phase t < -47 C
- ELSE IF ( t1 .LT. -47. ) THEN
- tk = t(i,k,j)
- rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
- c3 * (1. - tk / tf) + alog10(eis)
- ew = 10. ** rhs
- END IF
- ! Now sat vap pres obtained compute local vapor pressure
- ew = MAX ( ew , 0. ) * rh(i,k,j) * 0.01
- ! Now compute the specific humidity using the partial vapor
- ! pressures of water vapor (ew) and dry air (p-ew). The
- ! constants assume that the pressure is in hPa, so we divide
- ! the pressures by 100.
- q1 = mw_vap * ew
- q1 = q1 / (q1 + mw_air * (p(i,k,j)/100. - ew))
- q(i,k,j) = q1 / (1. - q1 )
- END IF
- END DO
- END DO
- END DO
- END IF
- ! For pressures above a defined level, reasonable Qv values should be
- ! a certain value or smaller. If they are larger than this, the input data
- ! probably had "missing" RH, and we filled in some values. This is an
- ! attempt to catch those. Also, set the minimum value for the entire
- ! domain that is above the selected pressure level.
-
- DO j = jts , MIN ( jde-1 , jte )
- DO k = kts , kte-1
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( p(i,k,j) .LT. qv_max_p_safe ) THEN
- IF ( q(i,k,j) .GT. qv_max_flag ) THEN
- q(i,k,j) = qv_max_value
- END IF
- END IF
- IF ( p(i,k,j) .LT. qv_min_p_safe ) THEN
- IF ( q(i,k,j) .LT. qv_min_flag ) THEN
- q(i,k,j) = qv_min_value
- END IF
- END IF
- END DO
- END DO
- END DO
- END SUBROUTINE rh_to_mxrat1
- !---------------------------------------------------------------------
- SUBROUTINE compute_eta ( znw , &
- eta_levels , max_eta , max_dz , &
- p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Compute eta levels, either using given values from the namelist (hardly
- ! a computation, yep, I know), or assuming a constant dz above the PBL,
- ! knowing p_top and the number of eta levels.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , INTENT(IN) :: max_dz
- REAL , INTENT(IN) :: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , tiso
- INTEGER , INTENT(IN) :: max_eta
- REAL , DIMENSION (max_eta) , INTENT(IN) :: eta_levels
- REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw
- ! Local vars
- INTEGER :: k
- REAL :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
- REAL , DIMENSION(kts:kte) :: dnw
- INTEGER , PARAMETER :: prac_levels = 17
- INTEGER :: loop , loop1
- REAL , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
- REAL , DIMENSION(kts:kte) :: alb , phb
- ! Gee, do the eta levels come in from the namelist?
- IF ( ABS(eta_levels(1)+1.) .GT. 0.0000001 ) THEN
- ! Check to see if the array is oriented OK, we can easily fix an upside down oops.
- IF ( ( ABS(eta_levels(1 )-1.) .LT. 0.0000001 ) .AND. &
- ( ABS(eta_levels(kde)-0.) .LT. 0.0000001 ) ) THEN
- DO k = kds+1 , kde-1
- znw(k) = eta_levels(k)
- END DO
- znw( 1) = 1.
- znw(kde) = 0.
- ELSE IF ( ( ABS(eta_levels(kde)-1.) .LT. 0.0000001 ) .AND. &
- ( ABS(eta_levels(1 )-0.) .LT. 0.0000001 ) ) THEN
- DO k = kds+1 , kde-1
- znw(k) = eta_levels(kde+1-k)
- END DO
- znw( 1) = 1.
- znw(kde) = 0.
- ELSE
- CALL wrf_error_fatal ( 'First eta level should be 1.0 and the last 0.0 in namelist' )
- END IF
- ! Check to see if the input full-level eta array is monotonic.
- DO k = kds , kde-1
- IF ( znw(k) .LE. znw(k+1) ) THEN
- PRINT *,'eta on full levels is not monotonic'
- PRINT *,'eta (',k,') = ',znw(k)
- PRINT *,'eta (',k+1,') = ',znw(k+1)
- CALL wrf_error_fatal ( 'Fix non-monotonic "eta_levels" in the namelist.input file' )
- END IF
- END DO
- ! Compute eta levels assuming a constant delta z above the PBL.
- ELSE
- ! Compute top of the atmosphere with some silly levels. We just want to
- ! integrate to get a reasonable value for ztop. We use the planned PBL-esque
- ! levels, and then just coarse resolution above that. We know p_top, and we
- ! have the base state vars.
- p_surf = p00
- znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
- 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
- DO k = 1 , prac_levels - 1
- znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5
- dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
- END DO
- DO k = 1, prac_levels-1
- pb = znu_prac(k)*(p_surf - p_top) + p_top
- temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
- ! temp = t00 + A*LOG(pb/p00)
- t_init = temp*(p00/pb)**(r_d/cp) - t0
- alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
- END DO
- ! Base state mu is defined as base state surface pressure minus p_top
- mub = p_surf - p_top
- ! Integrate base geopotential, starting at terrain elevation.
- phb(1) = 0.
- DO k = 2,prac_levels
- phb(k) = phb(k-1) - dnw_prac(k-1)*mub*alb(k-1)
- END DO
- ! So, now we know the model top in meters. Get the average depth above the PBL
- ! of each of the remaining levels. We are going for a constant delta z thickness.
- ztop = phb(prac_levels) / g
- ztop_pbl = phb(8 ) / g
- dz = ( ztop - ztop_pbl ) / REAL ( kde - 8 )
- ! Standard levels near the surface so no one gets in trouble.
- DO k = 1 , 8
- znw(k) = znw_prac(k)
- END DO
- ! Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
- ! Skamarock et al, NCAR TN 468. Use full levels, so
- ! use twice the thickness.
- DO k = 8, kte-1-2
- pb = znw(k) * (p_surf - p_top) + p_top
- temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
- ! temp = t00 + A*LOG(pb/p00)
- t_init = temp*(p00/pb)**(r_d/cp) - t0
- alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
- znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
- END DO
- znw(kte-2) = 0.000
- ! There is some iteration. We want the top level, ztop, to be
- ! consistent with the delta z, and we want the half level values
- ! to be consistent with the eta levels. The inner loop to 10 gets
- ! the eta levels very accurately, but has a residual at the top, due
- ! to dz changing. We reset dz five times, and then things seem OK.
- DO loop1 = 1 , 5
- DO loop = 1 , 10
- DO k = 8, kte-1-2
- pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
- temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
- ! temp = t00 + A*LOG(pb/p00)
- t_init = temp*(p00/pb)**(r_d/cp) - t0
- alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
- znw(k+1) = znw(k) - dz*g / ( mub*alb(k) )
- END DO
- IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
- print *,'Converged znw(kte) should be about 0.0 = ',znw(kte-2)
- END IF
- znw(kte-2) = 0.000
- END DO
- ! Here is where we check the eta levels values we just computed.
- DO k = 1, kde-1-2
- pb = (znw(k)+znw(k+1))*0.5 * (p_surf - p_top) + p_top
- temp = MAX ( tiso, t00 + A*LOG(pb/p00) )
- ! temp = t00 + A*LOG(pb/p00)
- t_init = temp*(p00/pb)**(r_d/cp) - t0
- alb(k) = (r_d/p1000mb)*(t_init+t0)*(pb/p1000mb)**cvpm
- END DO
- phb(1) = 0.
- DO k = 2,kde-2
- phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
- END DO
- ! Reset the model top and the dz, and iterate.
- ztop = phb(kde-2)/g
- ztop_pbl = phb(8)/g
- dz = ( ztop - ztop_pbl ) / REAL ( (kde-2) - 8 )
- END DO
- IF ( dz .GT. max_dz ) THEN
- print *,'z (m) = ',phb(1)/g
- do k = 2 ,kte-2
- print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
- end do
- print *,'dz (m) above fixed eta levels = ',dz
- print *,'namelist max_dz (m) = ',max_dz
- print *,'namelist p_top (Pa) = ',p_top
- CALL wrf_debug ( 0, 'You need one of three things:' )
- CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
- CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
- CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
- CALL wrf_debug ( 0, 'All are namelist options')
- CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
- END IF
- ! Add those 2 levels back into the middle, just above the 8 levels
- ! that semi define a boundary layer. After we open up the levels,
- ! then we just linearly interpolate in znw. So now levels 1-8 are
- ! specified as the fixed boundary layer levels given in this routine.
- ! The top levels, 12 through kte are those computed. The middle
- ! levels 9, 10, and 11 are equi-spaced in znw, and are each 1/2 the
- ! the znw thickness of levels 11 through 12.
- DO k = kte-2 , 9 , -1
- znw(k+2) = znw(k)
- END DO
- znw( 9) = 0.75 * znw( 8) + 0.25 * znw(12)
- znw(10) = 0.50 * znw( 8) + 0.50 * znw(12)
- znw(11) = 0.25 * znw( 8) + 0.75 * znw(12)
- END IF
- END SUBROUTINE compute_eta
- !---------------------------------------------------------------------
- SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Plow through each month, find the max, min values for each i,j.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
- ! Local vars
- INTEGER :: i , j , l
- REAL :: minner , maxxer
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- minner = field_in(i,1,j)
- maxxer = field_in(i,1,j)
- DO l = 2 , 12
- IF ( field_in(i,l,j) .LT. minner ) THEN
- minner = field_in(i,l,j)
- END IF
- IF ( field_in(i,l,j) .GT. maxxer ) THEN
- maxxer = field_in(i,l,j)
- END IF
- END DO
- field_min(i,j) = minner
- field_max(i,j) = maxxer
- END DO
- END DO
- END SUBROUTINE monthly_min_max
- !---------------------------------------------------------------------
- SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Linrarly in time interpolate data to a current valid time. The data is
- ! assumed to come in "monthly", valid at the 15th of every month.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- CHARACTER (LEN=24) , INTENT(IN) :: date_str
- REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
- REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
- ! Local vars
- INTEGER :: i , j , l
- INTEGER , DIMENSION(0:13) :: middle
- INTEGER :: target_julyr , target_julday , target_date
- INTEGER :: julyr , julday , int_month , month1 , month2
- REAL :: gmt
- CHARACTER (LEN=4) :: yr
- CHARACTER (LEN=2) :: mon , day15
- WRITE(day15,FMT='(I2.2)') 15
- DO l = 1 , 12
- WRITE(mon,FMT='(I2.2)') l
- CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
- middle(l) = julyr*1000 + julday
- END DO
- l = 0
- middle(l) = middle( 1) - 31
- l = 13
- middle(l) = middle(12) + 31
- CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
- target_date = target_julyr * 1000 + target_julday
- find_month : DO l = 0 , 12
- IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
- DO j = jts , MIN ( jde-1 , jte )
- DO i = its , MIN (ide-1 , ite )
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- int_month = l
- IF ( ( int_month .EQ. 0 ) .OR. ( int_month .EQ. 12 ) ) THEN
- month1 = 12
- month2 = 1
- ELSE
- month1 = int_month
- month2 = month1 + 1
- END IF
- field_out(i,j) = ( field_in(i,month2,j) * ( target_date - middle(l) ) + &
- field_in(i,month1,j) * ( middle(l+1) - target_date ) ) / &
- ( middle(l+1) - middle(l) )
- END DO
- END DO
- EXIT find_month
- END IF
- END DO find_month
- END SUBROUTINE monthly_interp_to_date
- !---------------------------------------------------------------------
- SUBROUTINE sfcprs (t, q, height, pslv, ter, avgsfct, p, &
- psfc, ez_method, &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Computes the surface pressure using the input height,
- ! temperature and q (already computed from relative
- ! humidity) on p surfaces. Sea level pressure is used
- ! to extrapolate a first guess.
- IMPLICIT NONE
- REAL, PARAMETER :: gamma = 6.5E-3
- REAL, PARAMETER :: pconst = 10000.0
- REAL, PARAMETER :: Rd = r_d
- REAL, PARAMETER :: TC = svpt0 + 17.5
- REAL, PARAMETER :: gammarg = gamma * Rd / g
- REAL, PARAMETER :: rov2 = Rd / 2.
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- LOGICAL , INTENT ( IN ) :: ez_method
- REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: pslv , ter, avgsfct
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
- INTEGER :: i
- INTEGER :: j
- INTEGER :: k
- INTEGER , DIMENSION (its:ite,jts:jte) :: k500 , k700 , k850
- LOGICAL :: l1
- LOGICAL :: l2
- LOGICAL :: l3
- LOGICAL :: OK
- REAL :: gamma78 ( its:ite,jts:jte )
- REAL :: gamma57 ( its:ite,jts:jte )
- REAL :: ht ( its:ite,jts:jte )
- REAL :: p1 ( its:ite,jts:jte )
- REAL :: t1 ( its:ite,jts:jte )
- REAL :: t500 ( its:ite,jts:jte )
- REAL :: t700 ( its:ite,jts:jte )
- REAL :: t850 ( its:ite,jts:jte )
- REAL :: tfixed ( its:ite,jts:jte )
- REAL :: tsfc ( its:ite,jts:jte )
- REAL :: tslv ( its:ite,jts:jte )
- ! We either compute the surface pressure from a time averaged surface temperature
- ! (what we will call the "easy way"), or we try to remove the diurnal impact on the
- ! surface temperature (what we will call the "other way"). Both are essentially
- ! corrections to a sea level pressure with a high-resolution topography field.
- IF ( ez_method ) THEN
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / avgsfct(i,j) ) ** ( - g / ( Rd * gamma ) )
- END DO
- END DO
- ELSE
- ! Find the locations of the 850, 700 and 500 mb levels.
- k850 = 0 ! find k at: P=850
- k700 = 0 ! P=700
- k500 = 0 ! P=500
- i = its
- j = jts
- DO k = kts+1 , kte
- IF (NINT(p(i,k,j)) .EQ. 85000) THEN
- k850(i,j) = k
- ELSE IF (NINT(p(i,k,j)) .EQ. 70000) THEN
- k700(i,j) = k
- ELSE IF (NINT(p(i,k,j)) .EQ. 50000) THEN
- k500(i,j) = k
- END IF
- END DO
- IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / t(i,1,j) ) ** ( - g / ( Rd * gamma ) )
- END DO
- END DO
- RETURN
- #if 0
- ! Possibly it is just that we have a generalized vertical coord, so we do not
- ! have the values exactly. Do a simple assignment to a close vertical level.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- DO k = kts+1 , kte-1
- IF ( ( p(i,k,j) - 85000. ) * ( p(i,k+1,j) - 85000. ) .LE. 0.0 ) THEN
- k850(i,j) = k
- END IF
- IF ( ( p(i,k,j) - 70000. ) * ( p(i,k+1,j) - 70000. ) .LE. 0.0 ) THEN
- k700(i,j) = k
- END IF
- IF ( ( p(i,k,j) - 50000. ) * ( p(i,k+1,j) - 50000. ) .LE. 0.0 ) THEN
- k500(i,j) = k
- END IF
- END DO
- END DO
- END DO
- ! If we *still* do not have the k levels, punt. I mean, we did try.
- OK = .TRUE.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
- OK = .FALSE.
- PRINT '(A)','(i,j) = ',i,j,' Error in finding p level for 850, 700 or 500 hPa.'
- DO K = kts+1 , kte
- PRINT '(A,I3,A,F10.2,A)','K = ',k,' PRESSURE = ',p(i,k,j),' Pa'
- END DO
- PRINT '(A)','Expected 850, 700, and 500 mb values, at least.'
- END IF
- END DO
- END DO
- IF ( .NOT. OK ) THEN
- CALL wrf_error_fatal ( 'wrong pressure levels' )
- END IF
- #endif
- ! We are here if the data is isobaric and we found the levels for 850, 700,
- ! and 500 mb right off the bat.
- ELSE
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- k850(i,j) = k850(its,jts)
- k700(i,j) = k700(its,jts)
- k500(i,j) = k500(its,jts)
- END DO
- END DO
- END IF
- ! The 850 hPa level of geopotential height is called something special.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ht(i,j) = height(i,k850(i,j),j)
- END DO
- END DO
- ! The variable ht is now -ter/ht(850 hPa). The plot thickens.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ht(i,j) = -ter(i,j) / ht(i,j)
- END DO
- END DO
- ! Make an isothermal assumption to get a first guess at the surface
- ! pressure. This is to tell us which levels to use for the lapse
- ! rates in a bit.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- psfc(i,j) = pslv(i,j) * (pslv(i,j) / p(i,k850(i,j),j)) ** ht(i,j)
- END DO
- END DO
- ! Get a pressure more than pconst Pa above the surface - p1. The
- ! p1 is the top of the level that we will use for our lapse rate
- ! computations.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
- p1(i,j) = 85000.
- ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0. ) THEN
- p1(i,j) = psfc(i,j) - pconst
- ELSE
- p1(i,j) = 50000.
- END IF
- END DO
- END DO
- ! Compute virtual temperatures for k850, k700, and k500 layers. Now
- ! you see why we wanted Q on pressure levels, it all is beginning
- ! to make sense.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- t850(i,j) = t(i,k850(i,j),j) * (1. + 0.608 * q(i,k850(i,j),j))
- t700(i,j) = t(i,k700(i,j),j) * (1. + 0.608 * q(i,k700(i,j),j))
- t500(i,j) = t(i,k500(i,j),j) * (1. + 0.608 * q(i,k500(i,j),j))
- END DO
- END DO
- ! Compute lapse rates between these three levels. These are
- ! environmental values for each (i,j).
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- gamma78(i,j) = ALOG(t850(i,j) / t700(i,j)) / ALOG (p(i,k850(i,j),j) / p(i,k700(i,j),j) )
- gamma57(i,j) = ALOG(t700(i,j) / t500(i,j)) / ALOG (p(i,k700(i,j),j) / p(i,k500(i,j),j) )
- END DO
- END DO
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
- t1(i,j) = t850(i,j)
- ELSE IF ( ( psfc(i,j) - 85000. ) .GE. 0. ) THEN
- t1(i,j) = t700(i,j) * (p1(i,j) / (p(i,k700(i,j),j))) ** gamma78(i,j)
- ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0.) THEN
- t1(i,j) = t500(i,j) * (p1(i,j) / (p(i,k500(i,j),j))) ** gamma57(i,j)
- ELSE
- t1(i,j) = t500(i,j)
- ENDIF
- END DO
- END DO
- ! From our temperature way up in the air, we extrapolate down to
- ! the sea level to get a guess at the sea level temperature.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- tslv(i,j) = t1(i,j) * (pslv(i,j) / p1(i,j)) ** gammarg
- END DO
- END DO
- ! The new surface temperature is computed from the with new sea level
- ! temperature, just using the elevation and a lapse rate. This lapse
- ! rate is -6.5 K/km.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- tsfc(i,j) = tslv(i,j) - gamma * ter(i,j)
- END DO
- END DO
- ! A small correction to the sea-level temperature, in case it is too warm.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- tfixed(i,j) = tc - 0.005 * (tsfc(i,j) - tc) ** 2
- END DO
- END DO
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- l1 = tslv(i,j) .LT. tc
- l2 = tsfc(i,j) .LE. tc
- l3 = .NOT. l1
- IF ( l2 .AND. l3 ) THEN
- tslv(i,j) = tc
- ELSE IF ( ( .NOT. l2 ) .AND. l3 ) THEN
- tslv(i,j) = tfixed(i,j)
- END IF
- END DO
- END DO
- ! Finally, we can get to the surface pressure.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- p1(i,j) = - ter(i,j) * g / ( rov2 * ( tsfc(i,j) + tslv(i,j) ) )
- psfc(i,j) = pslv(i,j) * EXP ( p1(i,j) )
- END DO
- END DO
- END IF
- ! Surface pressure and sea-level pressure are the same at sea level.
- ! DO j = jts , MIN(jde-1,jte)
- ! DO i = its , MIN(ide-1,ite)
- ! IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! IF ( ABS ( ter(i,j) ) .LT. 0.1 ) THEN
- ! psfc(i,j) = pslv(i,j)
- ! END IF
- ! END DO
- ! END DO
- END SUBROUTINE sfcprs
- !---------------------------------------------------------------------
- SUBROUTINE sfcprs2(t, q, height, psfc_in, ter, avgsfct, p, &
- psfc, ez_method, &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Computes the surface pressure using the input height,
- ! temperature and q (already computed from relative
- ! humidity) on p surfaces. Sea level pressure is used
- ! to extrapolate a first guess.
- IMPLICIT NONE
- REAL, PARAMETER :: Rd = r_d
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- LOGICAL , INTENT ( IN ) :: ez_method
- REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: psfc_in , ter, avgsfct
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
- INTEGER :: i
- INTEGER :: j
- INTEGER :: k
- REAL :: tv_sfc_avg , tv_sfc , del_z
- ! Compute the new surface pressure from the old surface pressure, and a
- ! known change in elevation at the surface.
- ! del_z = diff in surface topo, lo-res vs hi-res
- ! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
- IF ( ez_method ) THEN
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- tv_sfc_avg = avgsfct(i,j) * (1. + 0.608 * q(i,1,j))
- del_z = height(i,1,j) - ter(i,j)
- psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc_avg ) )
- END DO
- END DO
- ELSE
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- tv_sfc = t(i,1,j) * (1. + 0.608 * q(i,1,j))
- del_z = height(i,1,j) - ter(i,j)
- psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc ) )
- END DO
- END DO
- END IF
- END SUBROUTINE sfcprs2
- !---------------------------------------------------------------------
- SUBROUTINE sfcprs3( height , p , ter , slp , psfc , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- ! Computes the surface pressure by vertically interpolating
- ! linearly (or log) in z the pressure, to the targeted topography.
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: height, p
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: ter , slp
- REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
- INTEGER :: i
- INTEGER :: j
- INTEGER :: k
- LOGICAL :: found_loc
- REAL :: zl , zu , pl , pu , zm
- ! Loop over each grid point
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
- ! Special case where near the ocean level. Assume that the SLP is a good value.
- IF ( ter(i,j) .LT. 50 ) THEN
- psfc(i,j) = slp(i,j) + ( p(i,2,j)-p(i,3,j) ) / ( height(i,2,j)-height(i,3,j) ) * ter(i,j)
- CYCLE
- END IF
- ! Find the trapping levels
- found_loc = .FALSE.
- ! Normal sort of scenario - the model topography is somewhere between
- ! the height values of 1000 mb and the top of the model.
- found_k_loc : DO k = kts+1 , kte-2
- IF ( ( height(i,k ,j) .LE. ter(i,j) ) .AND. &
- ( height(i,k+1,j) .GT. ter(i,j) ) ) THEN
- zl = height(i,k ,j)
- zu = height(i,k+1,j)
- zm = ter(i,j)
- pl = p(i,k ,j)
- pu = p(i,k+1,j)
- psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
- found_loc = .TRUE.
- EXIT found_k_loc
- END IF
- END DO found_k_loc
- ! Interpolate betwixt slp and the first isobaric level above - this is probably the
- ! usual thing over the ocean.
- IF ( .NOT. found_loc ) THEN
- IF ( slp(i,j) .GE. p(i,2,j) ) THEN
- zl = 0.
- zu = height(i,3,j)
- zm = ter(i,j)
- pl = slp(i,j)
- pu = p(i,3,j)
- psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
- found_loc = .TRUE.
- ELSE
- found_slp_loc : DO k = kts+1 , kte-3
- IF ( ( slp(i,j) .GE. p(i,k+1,j) ) .AND. &
- ( slp(i,j) .LT. p(i,k ,j) ) ) THEN
- zl = 0.
- zu = height(i,k+1,j)
- zm = ter(i,j)
- pl = slp(i,j)
- pu = p(i,k+1,j)
- psfc(i,j) = EXP ( ( LOG(pl) * ( zm - zu ) + LOG(pu) * ( zl - zm ) ) / ( zl - zu ) )
- found_loc = .TRUE.
- EXIT found_slp_loc
- END IF
- END DO found_slp_loc
- END IF
- END IF
- ! Did we do what we wanted done.
- IF ( .NOT. found_loc ) THEN
- print *,'i,j = ',i,j
- print *,'p column = ',p(i,2:,j)
- print *,'z column = ',height(i,2:,j)
- print *,'model topo = ',ter(i,j)
- CALL wrf_error_fatal ( ' probs with sfc p computation ' )
- END IF
- END DO
- END DO
- END SUBROUTINE sfcprs3
- !---------------------------------------------------------------------
- SUBROUTINE filter_topo ( ht_in , xlat , msftx , fft_filter_lat , &
- ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte )
- IMPLICIT NONE
- INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
- ims , ime , jms , jme , kms , kme , &
- its , ite , jts , jte , kts , kte
- REAL , INTENT(IN) :: fft_filter_lat
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: ht_in
- REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: xlat , msftx
- ! Local vars
- INTEGER :: i , j , j_lat_pos , j_lat_neg
- INTEGER :: i_kicker , ik , i1, i2, i3, i4
- REAL :: length_scale , sum
- REAL , DIMENSION(its:ite,jts:jte) :: ht_out
- ! The filtering is a simple average on a latitude loop. Possibly a LONG list of
- ! numbers. We assume that ALL of the 2d arrays have been transposed so that
- ! each patch has the entire domain size of the i-dim local.
- IF ( ( its .NE. ids ) .OR. ( ite .NE. ide ) ) THEN
- CALL wrf_error_fatal ( 'filtering assumes all values on X' )
- END IF
- ! Starting at the south pole, we find where the
- ! grid distance is big enough, then go back a point. Continuing to the
- ! north pole, we find the first small grid distance. These are the
- ! computational latitude loops and the associated computational poles.
- j_lat_neg = 0
- j_lat_pos = jde + 1
- loop_neg : DO j = jts , MIN(jde-1,jte)
- IF ( xlat(its,j) .LT. 0.0 ) THEN
- IF ( ABS(xlat(its,j)) .LT. fft_filter_lat ) THEN
- j_lat_neg = j - 1
- EXIT loop_neg
- END IF
- END IF
- END DO loop_neg
- loop_pos : DO j = jts , MIN(jde-1,jte)
- IF ( xlat(its,j) .GT. 0.0 ) THEN
- IF ( xlat(its,j) .GE. fft_filter_lat ) THEN
- j_lat_pos = j
- EXIT loop_pos
- END IF
- END IF
- END DO loop_pos
- ! Set output values to initial input topo values for whole patch.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- ht_out(i,j) = ht_in(i,j)
- END DO
- END DO
- ! Filter the topo at the negative lats.
- DO j = j_lat_neg , jts , -1
- i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
- print *,'j = ' , j, ', kicker = ',i_kicker
- DO i = its , MIN(ide-1,ite)
- IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i+ik,j)
- END DO
- i1 = i - i_kicker + ide -1
- i2 = ide-1
- i3 = ids
- i4 = i-1
- DO ik = i1 , i2
- sum = sum + ht_in(ik,j)
- END DO
- DO ik = i3 , i4
- sum = sum + ht_in(ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i-ik,j)
- END DO
- i1 = i+1
- i2 = ide-1
- i3 = ids
- i4 = ids + ( i_kicker+i ) - ide
- DO ik = i1 , i2
- sum = sum + ht_in(ik,j)
- END DO
- DO ik = i3 , i4
- sum = sum + ht_in(ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- END IF
- END DO
- END DO
- ! Filter the topo at the positive lats.
- DO j = j_lat_pos , MIN(jde-1,jte)
- i_kicker = MIN( MAX ( NINT(msftx(its,j)) , 1 ) , (ide - ids) / 2 )
- print *,'j = ' , j, ', kicker = ',i_kicker
- DO i = its , MIN(ide-1,ite)
- IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i+ik,j) + ht_in(i-ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- ELSE IF ( ( i - i_kicker .LT. its ) .AND. ( i + i_kicker .LE. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i+ik,j)
- END DO
- i1 = i - i_kicker + ide -1
- i2 = ide-1
- i3 = ids
- i4 = i-1
- DO ik = i1 , i2
- sum = sum + ht_in(ik,j)
- END DO
- DO ik = i3 , i4
- sum = sum + ht_in(ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- ELSE IF ( ( i - i_kicker .GE. its ) .AND. ( i + i_kicker .GT. ide-1 ) ) THEN
- sum = 0.0
- DO ik = 1 , i_kicker
- sum = sum + ht_in(i-ik,j)
- END DO
- i1 = i+1
- i2 = ide-1
- i3 = ids
- i4 = ids + ( i_kicker+i ) - ide
- DO ik = i1 , i2
- sum = sum + ht_in(ik,j)
- END DO
- DO ik = i3 , i4
- sum = sum + ht_in(ik,j)
- END DO
- ht_out(i,j) = ( ht_in(i,j) + sum ) / REAL ( 2 * i_kicker + 1 )
- END IF
- END DO
- END DO
- ! Set output values to initial input topo values for whole patch.
- DO j = jts , MIN(jde-1,jte)
- DO i = its , MIN(ide-1,ite)
- ht_in(i,j) = ht_out(i,j)
- END DO
- END DO
- END SUBROUTINE filter_topo
- !---------------------------------------------------------------------
- SUBROUTINE init_module_initialize
- END SUBROUTINE init_module_initialize
- !---------------------------------------------------------------------
- END MODULE module_initialize_real
- #endif