!dis !dis Open Source License/Disclaimer, Forecast Systems Laboratory !dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305 !dis !dis This software is distributed under the Open Source Definition, !dis which may be found at http://www.opensource.org/osd.html. !dis !dis In particular, redistribution and use in source and binary forms, !dis with or without modification, are permitted provided that the !dis following conditions are met: !dis !dis - Redistributions of source code must retain this notice, this !dis list of conditions and the following disclaimer. !dis !dis - Redistributions in binary form must provide access to this !dis notice, this list of conditions and the following disclaimer, and !dis the underlying source code. !dis !dis - All modifications to this software must be clearly documented, !dis and are solely the responsibility of the agent making the !dis modifications. !dis !dis - If significant modifications or enhancements are made to this !dis software, the FSL Software Policy Manager !dis (softwaremgr@fsl.noaa.gov) should be notified. !dis !dis THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN !dis AND ARE FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES !dis GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND !dis AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS !dis OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME !dis NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND !dis DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS. !dis !dis MODULE vinterp_etap ! Module that handles the vertical interpolation of basic state variables ! to the Eta-P coordinate. ! ! EtaP = Pdry - Ptop Pdry = Dry pressure (Pa) ! ---------- Ptop = Pressure at model top (Pa) ! mu mu = Psfc - Pvapor_bottom - Ptop ! ! General Procedure: ! ! 1. Compute downward integrated vapor-pressure for each gridpoint ! ( Pvapor) ! 2. Subtract the Pvapor value from the pressure at each point to ! get the dry pressure (Pdry) at each point ! 3. Compute mu (dry surface pressure minus top pressure) ! 4. Compute EtaP value using newly compute dry pressure for each ! point ! 5. Interpolate grids to desired EtaP levels. ! 6. Diagnose WRF state variables ! 7. Horizontally stagger variables as required and output. ! ! HISTORY ! Feb 2002 - Original version -- B. Shaw, NOAA/FSL USE wrfsi_io USE wrfsi_static USE wrf_metadata USE vinterp_domain USE vinterp_setup USE date_pack USE grid_utils USE physical_constants USE diagnostic_vars USE vinterp_utils USE WriteField IMPLICIT NONE PRIVATE ! Variables shared among contained subroutines REAL, ALLOCATABLE :: etap_src(:,:,:) REAL, ALLOCATABLE :: mu(:,:) INTEGER, ALLOCATABLE :: sfc_k_ind(:,:) CHARACTER(LEN=8) :: processed_var_list(200) INTEGER :: num_processed LOGICAL :: use_sfc LOGICAL :: skip_lowest TYPE(wrfvar_metadata) :: var_meta_out PUBLIC etap_driver CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE etap_driver IMPLICIT NONE CALL setup_etap CALL interp_state_etap CALL process_others IF (ALLOCATED(etap_src)) DEALLOCATE(etap_src) IF (ALLOCATED(mu)) DEALLOCATE(mu) IF (ALLOCATED(sfc_k_ind)) DEALLOCATE(sfc_k_ind) print '(A,I2)','Finished EtaP driver for ', dom_meta%id RETURN END SUBROUTINE etap_driver !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE setup_etap IMPLICIT NONE character(LEN=4) fileType real, allocatable, dimension(:) :: level_1d integer :: kk num_processed = 0 processed_var_list(:) = ' ' dom_out = dom_meta dom_out%zdim = output_nz fileType = setup_info%output_file_type if(fileType == 'BIN' .or. fileType == 'BOTH') then CALL init_new_output_file(dom_out,setup_info%output_prefix, & setup_info%current_date) if(setup_info%output_file_type == 'BIN') DryRun=.false. endif if(fileType == 'WRF' .or. fileType == 'BOTH') then CALL init_new_wrf_output_file(dom_out,setup_info%wrf_output_prefix, & setup_info%current_date) endif ! Output the Full EtaP level values var_meta_out%name = etapfull_id var_meta_out%units = 'Dimensionless' var_meta_out%description = 'EtaP value on vertical full levels (mass)' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=1 var_meta_out%dim_val(1) = output_nz var_meta_out%dim_val(2) = 0 var_meta_out%dim_val(3) = 0 var_meta_out%dim_val(4) = 0 var_meta_out%dim_val(5) = 0 var_meta_out%dim_desc(1) = 'VERT' var_meta_out%dim_desc(2) = ' ' var_meta_out%dim_desc(3) = ' ' var_meta_out%dim_desc(4) = ' ' var_meta_out%dim_desc(5) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 0 var_meta_out%start_index(3) = 0 var_meta_out%start_index(4) = 0 var_meta_out%start_index(5) = 0 var_meta_out%stop_index(1) = output_nz var_meta_out%stop_index(2) = 0 var_meta_out%stop_index(3) = 0 var_meta_out%stop_index(4) = 0 var_meta_out%stop_index(5) = 0 var_meta_out%h_stagger_index = 0 var_meta_out%v_stagger_index = zstag_full_index var_meta_out%array_order = '+Z ' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'User defined levels generated by WRFSI/vinterp' var_meta_out%field_time_type = 'CONSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time num_processed = num_processed + 1 processed_var_list(num_processed) = etapfull_id allocate (level_1d(output_nz)) do kk = 1, output_nz level_1d(kk) = output_levels(kk,zstag_full_index) end do call write_field ( var_meta_out%name,level_1d,var_meta_out) ! call write_field ( var_meta_out%name,output_levels(:,zstag_full_index),var_meta_out) deallocate (level_1d) RETURN END SUBROUTINE setup_etap !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE interp_state_etap ! This routine computes the value of etap for the input data. IMPLICIT NONE INTEGER :: i,j,k REAL, ALLOCATABLE :: psfc_bg(:,:) LOGICAL :: psfc_bg_present LOGICAL :: found_k REAL :: deltap REAL :: wgt0 INTEGER :: adjustz_cnt REAL, ALLOCATABLE :: press3d_in(:,:,:) REAL, ALLOCATABLE :: t3d_in(:,:,:) REAL, ALLOCATABLE :: rh3d_in(:,:,:) REAL, ALLOCATABLE :: z3d_in(:,:,:) REAL, ALLOCATABLE :: t3d_out(:,:,:) REAL, ALLOCATABLE :: rh3d_out(:,:,:) REAL, ALLOCATABLE :: theta3d_in(:,:,:) REAL, ALLOCATABLE :: theta3d_out(:,:,:) REAL, ALLOCATABLE :: qv3d_in(:,:,:) REAL, ALLOCATABLE :: qv3d_out(:,:,:) REAL, ALLOCATABLE :: uv3d_in(:,:,:) REAL, ALLOCATABLE :: u3d_out(:,:,:) REAL, ALLOCATABLE :: v3d_out(:,:,:) REAL, ALLOCATABLE :: dum2ds(:,:) ! Initialize some variables psfc_bg_present = .false. ! Allocate space for the etap_src, mu, and temp/RH/pressure arrays. ! The input data ! has already been horizontally interpolated to the WRF non-staggered ! grid (A-grid), so the dimensions found in the dom_meta record ! are used for allocation. ASSUMPTION HERE: If surface fields ! are provided, they will be the first slab in the 3D array, regardless ! of vertical coordinate and whether or not there are 3D levels ! "below ground". This requires a sort on a column by column basis. The ! only exception is surface pressure, which is carried as a separate ! 2d field, which then allows us to move the surface slab to the ! appropriate position in each column. More to follow below... IF (ALLOCATED(etap_src)) DEALLOCATE(etap_src) ALLOCATE(etap_src (dom_meta%xdim,dom_meta%ydim,dom_meta%zdim )) IF (ALLOCATED(mu)) DEALLOCATE (mu) ALLOCATE(mu(dom_meta%xdim,dom_meta%ydim)) IF (ALLOCATED(press3d_in)) DEALLOCATE(press3d_in) ALLOCATE(press3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim )) IF (ALLOCATED(t3d_in)) DEALLOCATE(t3d_in) ALLOCATE(t3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim )) IF (ALLOCATED(rh3d_in)) DEALLOCATE(rh3d_in) ALLOCATE(rh3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim )) IF (ALLOCATED(z3d_in)) DEALLOCATE(z3d_in) ALLOCATE(z3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim )) IF (ALLOCATED(psfc_bg)) DEALLOCATE(psfc_bg) ALLOCATE(psfc_bg(dom_meta%xdim,dom_meta%ydim)) ! Get the temperature, RH, and surface pressure fields. ! Determine whether or not input data is isobaric. IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 3D temperature field...' CALL get_variable(setup_info%input_prefix, t_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN t3d_in = real_array(:,:,:,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = t_id ELSE PRINT *,'Problem getting ', t_id DEALLOCATE(real_array) STOP 'No_T3D_data' ENDIF if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(t3d_in),MAXVAL(t3d_in) IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 3D height field...' CALL get_variable(setup_info%input_prefix, height_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN z3d_in = real_array(:,:,:,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = height_id ELSE PRINT *,'Problem getting ', height_id STOP 'No_Z3D_data' ENDIF if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(z3d_in), & MAXVAL(z3d_in) IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 3D RH field...' CALL get_variable(setup_info%input_prefix, rh_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN rh3d_in = real_array(:,:,:,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = rh_id ELSE PRINT *,'Problem getting ', rh_id STOP 'No_RH3D_data' ENDIF if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(rh3d_in), & MAXVAL(rh3d_in) IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 2D PSFC_BG field...' CALL get_variable(setup_info%input_prefix, press_sfc_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN psfc_bg = real_array(:,:,1,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = press_sfc_id psfc_bg_present = .true. IF(.not.Dryrun)PRINT '(A,2F12.1)','Min/Max values = ',MINVAL(psfc_bg), & MAXVAL(psfc_bg) ELSE IF(.not.Dryrun)PRINT *,'No background surface pressure:', press_sfc_id IF(.not.Dryrun)PRINT *,'Not a big deal...we will make do without!' DEALLOCATE(psfc_bg) ENDIF ! Get the pressures IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting pressure array' CALL get_variable(setup_info%input_prefix, pressure_id, dom_meta%id, & setup_info%current_date,status) IF (status .NE. 0) THEN PRINT *, 'No pressure data found using: ',pressure_id STOP 'Need_pressure_data' ELSE num_processed = num_processed + 1 processed_var_list(num_processed) = pressure_id IF (var_info%ndim .eq. 1) THEN ! Isobaric input: ! Create a 3D array from the 1D aray IF(.not.Dryrun)PRINT *,"Creating 3D pressure array from 1D" IF(.not.Dryrun)PRINT *,'Levels found = ', real_array(:,1,1,1,1) DO j = 1,dom_meta%ydim DO i = 1,dom_meta%xdim press3d_in(i,j,:) = real_array(:,1,1,1,1) ENDDO ENDDO ELSE IF (var_info%ndim .eq. 3) THEN IF(.not.Dryrun)PRINT *, 'Pressure array already 3D' press3d_in = real_array(:,:,:,1,1) IF(.not.Dryrun)PRINT '(A,2F10.1)', 'Min/Max pressures: ',MINVAL(press3d_in), & MAXVAL(press3d_in) ELSE PRINT *, 'Pressure array not 1D or 3d; ndim = ', var_info%ndim STOP 'Help_I_am_not_ready_for_this' ENDIF ENDIF ! Do we have surface data? If so, it should be the bottom layer ! of all arrays. We require hinterp have a 200100 Pa level as ! the first level as a placeholder for surface data. If the bottom ! pressure level is not set to this value, then we assume no surface ! data is present in any of the 3D arrays and no sorting is required. If ! surface data is indicated as present, then we require hinterp to ! fill the lowest slab with a "missing" value for any of the other ! 3D variables for which no surface value was available, just so ! we know. IF (MAXVAL(press3d_in(:,:,1)).LT.200000.) THEN IF(.not.Dryrun)PRINT *, 'Appears as if no surface data is present.' IF(.not.Dryrun)PRINT *, 'Setting use_sfc = false and skipping array sorting.' use_sfc = .false. skip_lowest = .false. ELSE ! Allocate the array to keep track of the k index at which we ! insert the surface data IF (ALLOCATED(sfc_k_ind)) DEALLOCATE(sfc_k_ind) ALLOCATE (sfc_k_ind(dom_meta%xdim,dom_meta%ydim)) ! We must have both the surface height field AND the surface ! pressure field to use the surface data. Later we can add code ! to work if just one of those is available. IF ( (MINVAL(z3d_in(:,:,1)) .EQ. rmissing).OR. & (.NOT. psfc_bg_present) ) THEN if(.not.DryRun) then PRINT *, 'Surface data is unusable due to lack of both' PRINT *, 'a surface height and pressure value. Both are required' PRINT *, 'for now, but we will enhance the code later. For this' PRINT *, 'run, surface fields will be ignored.' endif skip_lowest = .true. use_sfc = .false. DEALLOCATE(sfc_k_ind) ELSE ! Sort by pressure and move height values. IF(.not.Dryrun)PRINT *, 'Sorting using PSFC' use_sfc = .true. skip_lowest = .false. ! Loop over horizontal domain DO j = 1,dom_meta%ydim DO i = 1,dom_meta%xdim found_k = .false. ! See if surface is already in the right place as the ! bottom slab of the array. IF (psfc_bg(i,j) .GE. press3d_in(i,j,2)) THEN sfc_k_ind(i,j) = 1 press3d_in(i,j,1) = psfc_bg(i,j) found_k = .true. ELSE ! Find the correct position find_sfc_k: DO k=2,dom_meta%zdim-1 IF ( (psfc_bg(i,j).LT.press3d_in(i,j,k)) .AND. & (psfc_bg(i,j).GE.press3d_in(i,j,k+1)) ) THEN sfc_k_ind(i,j) = k ! Shift all pressures below down 1 level, removing ! bottom-most pressure value press3d_in(i,j,1:k-1) = press3d_in(i,j,2:k) ! Insert surface value at correct level press3d_in(i,j,k) = psfc_bg(i,j) found_k = .true. EXIT find_sfc_k ENDIF ENDDO find_sfc_k ENDIF ! Make sure we successfully found insertion point! IF (.NOT. found_k) THEN PRINT *, 'Could not find where to insert surface pressure.' PRINT *, 'I/J = ', i,j PRINT *, 'Surface pressure = ', psfc_bg(i,j) PRINT *, 'Column pressures = ', press3d_in(i,j,:) STOP 'problem_sorting_pressure' ENDIF ENDDO ENDDO ! Make the height array consistent with the pressure array. IT ! is a special case, because heights must be monotonically increasing ! as pressure decreases. First, just call the shift_sfc_slab ! routine, regardless of whether or not the values are filled ! with rmissing, because our sanity check will correct missing values ! as well. CALL shift_sfc_slab(z3d_in) ! Here is the sanity check... adjustz_cnt = 0 DO j = 1,dom_meta%ydim DO i = 1,dom_meta%xdim k = sfc_k_ind(i,j) IF ( k .EQ. 1) THEN IF ( z3d_in(i,j,k).GT.z3d_in(i,j,k+1) ) THEN ! PRINT *, 'Height problem at i/j/k/k+1',i,j,k,k+1, 'Values =',& ! z3d_in(i,j,k), z3d_in(i,j,k+1) ! If we are within 10 mb, do a simple correction, otherwise ! stop with an error deltap = press3d_in(i,j,k) - press3d_in(i,j,k+1) IF ((deltap .LT. 1000.).AND.(deltap.GE.0)) THEN ! PRINT *, 'Applying simple correction.' adjustz_cnt = adjustz_cnt + 1 z3d_in(i,j,k) = z3d_in(i,j,k+1) - deltap/10. ELSE PRINT *, 'Cannot correct. Need to debug...' print *, 'i,j,k = ', i,j,k print *, 'z3d(i,j,k) z3d(i,j,k+1) = ', & z3d_in(i,j,k),z3d_in(i,j,k+1) print *, 'press(k), press(k+1) = ', press3d_in(i,j,k), & press3d_in(i,j,k+1) STOP 'bad_surface_height_vs_press' ENDIF ENDIF ELSE IF ( (z3d_in(i,j,k).LT.z3d_in(i,j,k-1)) .OR. & (z3d_in(i,j,k).GT.z3d_in(i,j,k+1)) ) THEN adjustz_cnt = adjustz_cnt+1 ! PRINT *, 'Bad height at i/j/k = ', i,j,k ! PRINT *, 'Z0/Z/Z1 = ', z3d_in(i,j,k-1:k+1) ! PRINT *, 'Replacing with interpolated value' wgt0 = 1.0 - ( ALOG( press3d_in(i,j,k-1)/press3d_in(i,j,k))/& ALOG( press3d_in(i,j,k-1)/press3d_in(i,j,k+1))) z3d_in(i,j,k) = wgt0*z3d_in(i,j,k-1) + & (1.-wgt0)*z3d_in(i,j,k+1) ! PRINT *, 'New Z0/Z/Z1 = ',z3d_in(i,j,k-1:k+1) ENDIF ENDIF ENDDO ENDDO IF(.not.Dryrun)PRINT *, 'Number of columns where Zsfc had to be adjusted: ', & adjustz_cnt ! ELSE ! Sort by height and derive surface pressure ! COMPLETE THIS WORK LATER. FOR NOW, JUST MENTION ! NOT SUPPORTED AND SET USE_SFC = .FALSE. ! PRINT *, 'Sorting by ZSFC and deriving PSFC' ! use_sfc = .true. ! skip_lowest = .false. ENDIF ENDIF ! Now, if required, shift the surface levels of temperature ! and RH into their appropriate position in each column. IF (use_sfc) THEN IF (MINVAL(t3d_in(:,:,1)).NE.rmissing) THEN CALL shift_sfc_slab(t3d_in) ELSE PRINT *, 'Sfc Temperature field missing.' STOP 'need_to_add_support' ENDIF IF (MINVAL(rh3d_in(:,:,1)).NE.rmissing) THEN CALL shift_sfc_slab(rh3d_in) ELSE PRINT *, 'Sfc RH field missing.' STOP 'need_to_add_support' ENDIF ENDIF ! Generate the 3d array of Eta values on the background gridpoints. ! Subroutine also returns the 2D mu array, which is computed using ! the WRF topography to ensure a pressure field consistent with the ! forecast model. CALL compute_eta_3d(dom_meta%xdim, dom_meta%ydim, dom_meta%zdim, & press3d_in, t3d_in, z3d_in, rh3d_in, dom_meta%top_level, & terrain_hgt_n, etap_src, mu ) ! Diagnose potential temperature (theta) and Qv ALLOCATE (theta3d_in (dom_out%xdim,dom_out%ydim,dom_meta%zdim)) ALLOCATE(qv3d_in(dom_out%xdim,dom_out%ydim,dom_meta%zdim)) DO j = 1, dom_meta%ydim DO i = 1, dom_meta%xdim DO k = 1, dom_meta%zdim theta3d_in(i,j,k) = compute_theta(t3d_in(i,j,k),press3d_in(i,j,k)) qv3d_in(i,j,k) = compute_vapor_mixing_ratio(t3d_in(i,j,k), & press3d_in(i,j,k),rh3d_in(i,j,k),.true.) ENDDO ENDDO ENDDO ! At this point, we can free up some memory that is no longer needed DEALLOCATE(z3d_in) DEALLOCATE(t3d_in) DEALLOCATE(rh3d_in) DEALLOCATE(press3d_in) IF (ALLOCATED(psfc_bg)) DEALLOCATE (psfc_bg) ! Allocate space for output state thermdynamic variables ALLOCATE (qv3d_out(dom_out%xdim,dom_out%ydim,dom_out%zdim)) ALLOCATE (theta3d_out(dom_out%xdim,dom_out%ydim,dom_out%zdim)) CALL interp_eta2eta_lin(dom_meta%xdim, dom_meta%ydim, & dom_meta%zdim, dom_out%zdim, & etap_src,output_levels(:,zstag_half_index), & mu,dom_out%top_level,theta3d_in,theta3d_out,'THETA ') DEALLOCATE(theta3d_in) CALL interp_eta2eta_lin(dom_meta%xdim, dom_meta%ydim, & dom_meta%zdim, dom_out%zdim, & etap_src,output_levels(:,zstag_half_index), & mu,dom_out%top_level,qv3d_in,qv3d_out,'QVAPOR ') DEALLOCATE(qv3d_in) WHERE(qv3d_out .LT. 0.) qv3d_out = 0.0 ! Begin writing out data here so we can free up memory as it becomes ! possible. Note that staggering will be handled in output routines ! contained in this module CALL output_basic_thermo(theta3d_out,qv3d_out) DEALLOCATE(theta3d_out) DEALLOCATE(qv3d_out) ALLOCATE (u3d_out(dom_out%xdim,dom_out%ydim,dom_out%zdim)) ALLOCATE (v3d_out(dom_out%xdim,dom_out%ydim,dom_out%zdim)) ALLOCATE (uv3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim)) ! Process the u wind field IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 3D U field...' CALL get_variable(setup_info%input_prefix, u_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN uv3d_in = real_array(:,:,:,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = u_id ELSE PRINT *,'Problem getting ', u_id STOP 'No_U3D_data' ENDIF IF(.not.Dryrun)PRINT '(A,2F9.1)','Min/Max values = ',MINVAL(uv3d_in), & MAXVAL(uv3d_in) IF (use_sfc) THEN IF (MINVAL(uv3d_in(:,:,1)).NE.rmissing) THEN CALL shift_sfc_slab(uv3d_in) ELSE PRINT *, 'Sfc U field missing.' STOP 'need_to_add_support' ENDIF ENDIF CALL interp_eta2eta_lin(dom_meta%xdim, dom_meta%ydim, & dom_meta%zdim, dom_out%zdim, & etap_src,output_levels(:,zstag_half_index), & mu,dom_out%top_level,uv3d_in,u3d_out,'U ') ! Process the v wind field IF(.not.Dryrun)PRINT '(A,A)', 'INTERP_STATE_ETAP: Getting 3D V field...' CALL get_variable(setup_info%input_prefix, v_id, dom_meta%id, & setup_info%current_date,status) IF (status.EQ.0) THEN uv3d_in = real_array(:,:,:,1,1) num_processed = num_processed + 1 processed_var_list(num_processed) = v_id ELSE PRINT *,'Problem getting ', v_id STOP 'No_V3D_data' ENDIF IF(.not.Dryrun)PRINT '(A,2F9.1)','Min/Max values = ',MINVAL(uv3d_in), & MAXVAL(uv3d_in) IF (use_sfc) THEN IF (MINVAL(uv3d_in(:,:,1)).NE.rmissing) THEN CALL shift_sfc_slab(uv3d_in) ELSE PRINT *, 'Sfc V field missing.' STOP 'need_to_add_support' ENDIF ENDIF CALL interp_eta2eta_lin(dom_meta%xdim, dom_meta%ydim, & dom_meta%zdim, dom_out%zdim, & etap_src,output_levels(:,zstag_half_index), & mu,dom_out%top_level,uv3d_in,v3d_out,'V ') DEALLOCATE (uv3d_in) ! Output the wind field (staggering done in output routine) CALL output_basic_momentum(u3d_out,v3d_out) DEALLOCATE(u3d_out) DEALLOCATE(v3d_out) RETURN END SUBROUTINE interp_state_etap !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE shift_sfc_slab(data3d) IMPLICIT NONE REAL, INTENT(INOUT) :: data3d(:,:,:) REAL :: sfcval INTEGER :: i,j DO j = 1 , dom_meta%ydim DO i = 1 , dom_meta%xdim IF (sfc_k_ind(i,j) .GT. 1) THEN sfcval = data3d(i,j,1) data3d(i,j,1:sfc_k_ind(i,j)-1) = data3d(i,j,2:sfc_k_ind(i,j)) data3d(i,j,sfc_k_ind(i,j)) = sfcval ENDIF ENDDO ENDDO RETURN END SUBROUTINE shift_sfc_slab !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE output_basic_thermo(theta,qv) IMPLICIT NONE REAL, INTENT(INOUT) :: theta(:,:,:) REAL, INTENT(INOUT) :: qv(:,:,:) REAL, ALLOCATABLE :: dum3ds(:,:,:) REAL, ALLOCATABLE :: mu_out(:,:) ALLOCATE(mu_out(dom_meta%xdim,dom_meta%ydim)) IF (setup_info%output_vars .LE. 2) THEN ! Stagger these fields to T grid CALL arakawa_c_n2t(mu,dom_meta%xdim,dom_meta%ydim,1, & mu_out) ALLOCATE(dum3ds(dom_meta%xdim,dom_meta%ydim,dom_out%zdim)) CALL arakawa_c_n2t(theta,dom_meta%xdim,dom_meta%ydim,dom_out%zdim, & dum3ds) theta = dum3ds CALL arakawa_c_n2t(qv,dom_meta%xdim,dom_meta%ydim,dom_out%zdim, & dum3ds) qv = dum3ds var_meta_out%h_stagger_index = t_ind DEALLOCATE (dum3ds) ELSE mu_out = mu var_meta_out%h_stagger_index = n_ind ENDIF ! Output mu var_meta_out%name = mumass_id var_meta_out%units = 'Pa ' var_meta_out%description = 'Surface scaled dry pressure ' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=2 var_meta_out%dim_val(1) = dom_meta%xdim var_meta_out%dim_val(2) = dom_meta%ydim var_meta_out%dim_val(3:var_maxdims) = 1 var_meta_out%dim_desc(1) = 'E-W' var_meta_out%dim_desc(2) = 'N-S' var_meta_out%dim_desc(3:var_maxdims) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 1 var_meta_out%start_index(3:var_maxdims) = 0 var_meta_out%stop_index(1) = dom_meta%xdim var_meta_out%stop_index(2) = dom_meta%ydim var_meta_out%stop_index(3:var_maxdims) = 0 var_meta_out%v_stagger_index = 0 var_meta_out%array_order = '+X+Y ' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'Vertical derivation by SI' var_meta_out%field_time_type = 'INSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time call write_field ( var_meta_out%name, mu_out,var_meta_out) DEALLOCATE(mu_out) ! Output theta var_meta_out%name = theta_id var_meta_out%units = 'K ' var_meta_out%description = 'Potential temperature ' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=3 var_meta_out%dim_val(1) = dom_meta%xdim var_meta_out%dim_val(2) = dom_meta%ydim var_meta_out%dim_val(3) = dom_out%zdim var_meta_out%dim_val(4:var_maxdims) = 1 var_meta_out%dim_desc(1) = 'E-W' var_meta_out%dim_desc(2) = 'N-S' var_meta_out%dim_desc(3) = 'VERT' var_meta_out%dim_desc(4:var_maxdims) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 1 var_meta_out%start_index(3) = 1 var_meta_out%start_index(4:var_maxdims) = 0 var_meta_out%stop_index(1) = dom_meta%xdim var_meta_out%stop_index(2) = dom_meta%ydim var_meta_out%stop_index(3) = dom_out%zdim var_meta_out%stop_index(4:var_maxdims) = 0 var_meta_out%v_stagger_index = zstag_half_index var_meta_out%array_order = '+X+Y+Z' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'Vertical interpolation by SI' var_meta_out%field_time_type = 'INSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time call write_field ( var_meta_out%name, theta, var_meta_out, 'z ' ) ! Output qv var_meta_out%name = qvapor_id var_meta_out%units = 'kg kg{-1}' var_meta_out%description = 'Water vapor mixing ratio ' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=3 var_meta_out%dim_val(1) = dom_meta%xdim var_meta_out%dim_val(2) = dom_meta%ydim var_meta_out%dim_val(3) = dom_out%zdim var_meta_out%dim_val(4:var_maxdims) = 1 var_meta_out%dim_desc(1) = 'E-W' var_meta_out%dim_desc(2) = 'N-S' var_meta_out%dim_desc(3) = 'VERT' var_meta_out%dim_desc(4:var_maxdims) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 1 var_meta_out%start_index(3) = 1 var_meta_out%start_index(4:var_maxdims) = 0 var_meta_out%stop_index(1) = dom_meta%xdim var_meta_out%stop_index(2) = dom_meta%ydim var_meta_out%stop_index(3) = dom_out%zdim var_meta_out%stop_index(4:var_maxdims) = 0 var_meta_out%v_stagger_index = zstag_half_index var_meta_out%array_order = '+X+Y+Z' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'Vertical interpolation by SI' var_meta_out%field_time_type = 'INSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time call write_field ( var_meta_out%name, qv, var_meta_out, 'z ' ) RETURN END SUBROUTINE output_basic_thermo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE output_basic_momentum(u3d,v3d) IMPLICIT NONE REAL,INTENT(INOUT) :: u3d(:,:,:) REAL,INTENT(INOUT) :: v3d(:,:,:) REAL, ALLOCATABLE :: dum3ds(:,:,:) INTEGER :: ustag,vstag IF (setup_info%output_vars .LE. 2) THEN ! Stagger these fields to U and V grids ALLOCATE(dum3ds(dom_meta%xdim,dom_meta%ydim,dom_out%zdim)) CALL arakawa_c_n2u(u3d,dom_meta%xdim,dom_meta%ydim,dom_out%zdim, & dum3ds) u3d = dum3ds CALL arakawa_c_n2v(v3d,dom_meta%xdim,dom_meta%ydim,dom_out%zdim, & dum3ds) v3d = dum3ds DEALLOCATE (dum3ds) ustag = u_ind vstag = v_ind ELSE ustag = n_ind vstag = n_ind ENDIF ! Output u3d var_meta_out%name = u_id var_meta_out%units = 'm s{-1}' var_meta_out%description = 'Grid-relative u-component of wind ' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=3 var_meta_out%dim_val(1) = dom_meta%xdim var_meta_out%dim_val(2) = dom_meta%ydim var_meta_out%dim_val(3) = dom_out%zdim var_meta_out%dim_val(4:var_maxdims) = 1 var_meta_out%dim_desc(1) = 'E-W' var_meta_out%dim_desc(2) = 'N-S' var_meta_out%dim_desc(3) = 'VERT' var_meta_out%dim_desc(4:var_maxdims) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 1 var_meta_out%start_index(3) = 1 var_meta_out%start_index(4:var_maxdims) = 0 var_meta_out%stop_index(1) = dom_meta%xdim var_meta_out%stop_index(2) = dom_meta%ydim var_meta_out%stop_index(3) = dom_out%zdim var_meta_out%stop_index(4:var_maxdims) = 0 var_meta_out%h_stagger_index = ustag var_meta_out%v_stagger_index = zstag_half_index var_meta_out%array_order = '+X+Y+Z' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'Vertical interpolation by SI' var_meta_out%field_time_type = 'INSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time call write_field ( var_meta_out%name, u3d, var_meta_out, 'z ' ) ! Output u3d var_meta_out%name = v_id var_meta_out%units = 'm s{-1}' var_meta_out%description = 'Grid-relative v-component of wind ' var_meta_out%domain_id = dom_meta%id var_meta_out%ndim=3 var_meta_out%dim_val(1) = dom_meta%xdim var_meta_out%dim_val(2) = dom_meta%ydim var_meta_out%dim_val(3) = dom_out%zdim var_meta_out%dim_val(4:var_maxdims) = 1 var_meta_out%dim_desc(1) = 'E-W' var_meta_out%dim_desc(2) = 'N-S' var_meta_out%dim_desc(3) = 'VERT' var_meta_out%dim_desc(4:var_maxdims) = ' ' var_meta_out%start_index(1) = 1 var_meta_out%start_index(2) = 1 var_meta_out%start_index(3) = 1 var_meta_out%start_index(4:var_maxdims) = 0 var_meta_out%stop_index(1) = dom_meta%xdim var_meta_out%stop_index(2) = dom_meta%ydim var_meta_out%stop_index(3) = dom_out%zdim var_meta_out%stop_index(4:var_maxdims) = 0 var_meta_out%h_stagger_index = vstag var_meta_out%v_stagger_index = zstag_half_index var_meta_out%array_order = '+X+Y+Z' var_meta_out%field_type = 'REAL' var_meta_out%field_source_prog = 'SI' var_meta_out%source_desc = 'Vertical interpolation by SI' var_meta_out%field_time_type = 'INSTANT' var_meta_out%vt_date_start = dom_meta%vt_date var_meta_out%vt_time_start = dom_meta%vt_time var_meta_out%vt_date_stop = dom_meta%vt_date var_meta_out%vt_time_stop = dom_meta%vt_time call write_field ( var_meta_out%name,v3d, var_meta_out, 'z ' ) RETURN END SUBROUTINE output_basic_momentum !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 SUBROUTINE process_others ! This routine scans the file looking for any variables that have not been ! processed yet. It looks at the dimensions, etc. to determine if the ! variable needs to be vertically interpolated and/or destaggered into ! multiple output fields. If not, it just passes them on, metadata included ! into the new output file. IMPLICIT NONE INTEGER :: v,s,k REAL, ALLOCATABLE :: dum3d_in(:,:,:),dum3d(:,:,:),dum3ds(:,:,:) REAL, ALLOCATABLE :: dum2d(:,:),dum2ds(:,:) CHARACTER (LEN=1) :: stag_string LOGICAL :: processed ! First, reset the file pointer to the beginning of the input file by ! calling get_domain_metadata. The API should have something like ! get_first_variable to eventually replace this kludge method. CALL get_domain_metadata(dom_meta%id, setup_info%current_date, & setup_info%input_prefix) IF(.not.Dryrun)PRINT *, 'Here is what we have processed so far: ' dumplist: DO k = 1, 200 IF(.not.Dryrun)PRINT *, processed_var_list(k) IF (processed_var_list(k) .EQ. ' ') EXIT dumplist ENDDO dumplist ! Now, we can loop through calling read_next_variable and using the ! status to determine if we have reached the end of the file (module_wrf_io ! sets the status=2 when EOF is reached) status = 0 CALL read_next_variable filescan_loop: DO WHILE (status .EQ. 0) ! Compare the variable name with those in the processed_var list processed = .false. check_var_loop: DO v = 1,num_processed IF (var_info%name .EQ. processed_var_list(v))THEN processed = .true. EXIT check_var_loop ENDIF ENDDO check_var_loop IF (processed) THEN CALL read_next_variable CYCLE filescan_loop ENDIF IF(.not.Dryrun)PRINT '(2A)', 'Other Variable: ', var_info%name ! We have a variable that has not been processed, so figure out ! what we need to do with it. IF ((var_info%ndim .EQ. 3).AND.(var_info%array_order.EQ.'+X+Y+Z '))THEN ALLOCATE (dum3d (dom_out%xdim,dom_out%ydim,dom_out%zdim)) ! if(.not.DryRun) then ALLOCATE (dum3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim)) dum3d_in = real_array(:,:,:,1,1) IF (use_sfc) THEN CALL shift_sfc_slab(dum3d_in) ENDIF CALL interp_eta2eta_lin(dom_meta%xdim, dom_meta%ydim, & dom_meta%zdim, dom_out%zdim, & etap_src,output_levels(:,zstag_half_index), & mu,dom_out%top_level,dum3d_in,dum3d,var_info%name) var_meta_out = var_info var_meta_out%dim_val(3) = output_nz var_meta_out%stop_index(3) = output_nz IF ( (var_meta_out%h_stagger_index .EQ. n_ind).AND. & (setup_info%output_vars .LE. 2) )THEN ALLOCATE(dum3ds(dom_out%xdim,dom_out%ydim,dom_out%zdim)) IF ( (var_meta_out%name(1:1) .NE. 'U') .AND. & (var_meta_out%name(1:1) .NE. 'V')) THEN CALL arakawa_c_n2t(dum3d,dom_out%xdim,dom_out%ydim, & dom_out%zdim,dum3ds) var_meta_out%h_stagger_index = t_ind ELSE IF (var_meta_out%name(1:1) .EQ. 'U') THEN CALL arakawa_c_n2u(dum3d,dom_out%xdim,dom_out%ydim, & dom_out%zdim,dum3ds) var_meta_out%h_stagger_index = u_ind ELSE IF (var_meta_out%name(1:1) .EQ. 'V') THEN CALL arakawa_c_n2v(dum3d,dom_out%xdim,dom_out%ydim, & dom_out%zdim,dum3ds) var_meta_out%h_stagger_index = v_ind ENDIF dum3d = dum3ds DEALLOCATE(dum3ds) ENDIF DEALLOCATE(dum3d_in) ! endif ! .not.DryRun call write_field ( var_info%name, dum3d, var_meta_out, 'z ' ) DEALLOCATE(dum3d) ELSE IF ((var_info%ndim .EQ. 3).AND. & (var_info%array_order .EQ. '+X+Y+S '))THEN ! Staggered 2D variable stag_output_loop: DO s=1,var_info%dim_val(3) var_meta_out = var_info ! Replace the last character of the variable name with stagger index WRITE (stag_string, '(I1)') s var_meta_out%name(8:8) = stag_string var_meta_out%description = TRIM(var_meta_out%description) // & ' on stagger number ' // stag_string var_meta_out%ndim = 2 var_meta_out%dim_val(3) = 0 var_meta_out%stop_index(3) = 0 var_meta_out%h_stagger_index = s var_meta_out%array_order = '+X+Y ' var_meta_out%dim_desc(3) = ' ' call write_field ( var_info%name, real_array(:,:,s,1,1), var_meta_out) ENDDO stag_output_loop ELSE IF ( var_info%ndim .EQ. 2) THEN ALLOCATE(dum2d(dom_out%xdim,dom_out%ydim)) dum2d = real_array(:,:,1,1,1) ! Pass this variable right on through var_meta_out = var_info ! If this variable is background landusec or soilctop, ! we need to set the flags to let our static processor ! know we already have it IF (var_info%name(1:6) .EQ. 'VEGCAT') have_bg_landusec = .true. IF (var_info%name(1:7) .EQ. 'SOILCAT') have_bg_soilctop = .true. IF ((var_meta_out%h_stagger_index .EQ. n_ind).AND. & (setup_info%output_vars .LE. 2) ) THEN ALLOCATE (dum2ds(dom_out%xdim,dom_out%ydim)) IF ( (var_meta_out%name(1:1) .NE. 'U') .AND. & (var_meta_out%name(1:1) .NE. 'V')) THEN CALL arakawa_c_n2t(dum2d,dom_out%xdim,dom_out%ydim,1,dum2ds) var_meta_out%h_stagger_index = t_ind ELSE IF (var_meta_out%name(1:1) .EQ. 'U') THEN CALL arakawa_c_n2u(dum2d,dom_out%xdim,dom_out%ydim,1,dum2ds) var_meta_out%h_stagger_index = u_ind ELSE IF (var_meta_out%name(1:1) .EQ. 'V') THEN CALL arakawa_c_n2v(dum2d,dom_out%xdim,dom_out%ydim,1,dum2ds) var_meta_out%h_stagger_index = v_ind ENDIF dum2d = dum2ds DEALLOCATE (dum2ds) ENDIF call write_field ( var_info%name, dum2d , var_meta_out) DEALLOCATE(dum2d) ELSE IF ( var_info%ndim .EQ. 1) THEN ! Pass this variable right on through var_meta_out = var_info call write_field ( var_info%name, real_array(:,1,1,1,1), var_meta_out) ELSE IF (var_info%ndim .EQ. 0) THEN ! Pass this variable right on through var_meta_out = var_info call write_field ( var_info%name, real_array(1,1,1,1,1), var_meta_out) ELSE IF(.not.Dryrun)PRINT '(A)', 'Unsure how to process this variable...omitting from output:' IF(.not.Dryrun)PRINT '(3A,I1,2A)', 'NAME: ', var_info%name, ' NDIM: ', var_info%ndim,& ' ARRAY ORDER: ', var_info%array_order ENDIF CALL read_next_variable ENDDO filescan_loop IF(.not.Dryrun)PRINT '(A)', 'ETAP_PROCESS_OTHERS: Reached end of file.' RETURN END SUBROUTINE process_others !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE vinterp_etap