/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
Large files files are truncated, but you can click here to view the full file
- !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. Tempo…
Large files files are truncated, but you can click here to view the full file