!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 hinterp_gribprep ! Module containing routines to process the gribprep output ! in hinterp USE map_utils USE wrf_metadata USE horiz_interp USE hinterp_setup USE domain_info USE gridded_data USE gribprep_io IMPLICIT NONE PRIVATE ! Definitions to set up linked list for ingest ! of 2D fields. TYPE twod REAL :: level REAL , DIMENSION(:,:) , POINTER :: field TYPE(twod) , POINTER :: next_slab END TYPE twod TYPE vars INTEGER :: num_levels CHARACTER (LEN=80) :: name TYPE(vars) , POINTER :: next_variable TYPE(twod) , POINTER :: slab END TYPE vars TYPE(vars) , TARGET :: head TYPE(vars) , POINTER :: current_var , new_var , old_var TYPE(twod) , POINTER :: current_slab , new_slab , old_slab LOGICAL :: initialized = .FALSE. , found = .FALSE. REAL , DIMENSION(:,:) , POINTER :: new_field ! Stuff for non-isobaric input fields. For now, we cannot tolerate ! isobaric and non-isobaric 3D fields in the same time period. LOGICAL :: have_iso3d, have_noniso3d, have_sfc, have_soilhgt PUBLIC ingest_gribprep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE ingest_gribprep(time,data_processed) ! This routine ingests the gribprep output files and interpolates the ! 2D slabs to the WRF grid. Meteorological fields are interpolated ! to the non-staggered grid, while fields that depend upon the ! land-water mask (e.g., SOIL*, etc.) are interpolated directly to ! the mass gridpoints, which may be different in the case of ! a staggered grid. ! ! 1. Input data is in gribprep format ! 2. Input data is one of the following projections: ! Cylindrical Equidistant (latlon grid) ! Mercator ! Lambert Conformal (tangent or secant) ! Polar Stereographic ! ! 3. Grid spacing must be same in both directions, except in the case ! of a latlon grid, even though the gribprep header has separate ! fields for dx and dy ! 4. Array must be ordered (i,j), where i is the East-West (x) direction ! and j is the North-South (y) direction. ! 5. Except in the case of a latlon grid, point (1,1) must correspond to ! the SW corner of the input data, and I is positive east and j is ! positive north. For latlon grids, the code looks at the sign ! of deltalat and deltalon to determine array order. ! IMPLICIT NONE CHARACTER(LEN=19),INTENT(IN) :: time LOGICAL, INTENT(OUT) :: data_processed REAL, DIMENSION(:,:), ALLOCATABLE :: array , & dum2d REAL :: xlo, & yla, & dlon, & dlat, & startlat, & startlon INTEGER :: ilo , jla , i , j , k, files_per_time LOGICAL :: initial_fg_read = .TRUE. LOGICAL :: found_any_files = .FALSE. INTEGER :: close_count , i_close , j_close REAL :: close_data REAL , DIMENSION(:,:), ALLOCATABLE :: landmask_fg LOGICAL :: landmask_fg_init, found_var REAL , DIMENSION (:,:), ALLOCATABLE :: xloc , yloc, xloc_t, yloc_t CHARACTER (LEN=24) :: hdate_input CHARACTER (LEN= 9) :: field_input CHARACTER (LEN=25) :: units_input CHARACTER (LEN=46) :: description_input LOGICAL :: constant_file TYPE (proj_info) :: proj_fg_new, proj_fg_old LOGICAL :: no_interp_needed LOGICAL :: same_as_prev INTEGER :: version CHARACTER(LEN=19) :: hdate REAL :: xfcst CHARACTER(LEN=32) :: source CHARACTER(LEN=8) :: field CHARACTER(LEN=24) :: units CHARACTER(LEN=46) :: description INTEGER :: llflag INTEGER :: idim, jdim REAL :: level CHARACTER(LEN=8) :: startloc REAL :: knownlat, knownlon, deltalat, deltalon REAL :: lat1, lon1, deltax, deltay, xlonc, truelat1, truelat2 INTEGER :: eof_status REAL, ALLOCATABLE :: u_temp(:,:), v_temp(:,:) LOGICAL :: rotate_fg_winds, rotate_wrf_winds INTEGER :: u_index, v_index INTEGER :: data_type REAL :: min_valid, max_valid, default_val INTEGER :: fg_unit LOGICAL :: already_opened INTEGER :: open_status INTEGER :: file_name_len CHARACTER(LEN=255) :: file_name INTEGER :: diff_files_per_single_time, time_index REAL :: val_mask LOGICAL :: make_landmask_fg INTEGER :: add_one_for_lsm REAL :: min_rh_in, max_rh_in ! For each time period that we ingest, set up the different ! list of names. This allows us to count how many 3d/2d fields ! we have encountered. IF (separate_lsm) THEN add_one_for_lsm = 1 ELSE add_one_for_lsm = 0 ENDIF all_names_sfc = ' ' all_names_up = ' ' num3d_fg = 0 num2d_fg = 0 num1d_fg = 0 have_iso3d = .false. have_noniso3d = .false. have_sfc = .false. ! We know how big the model grid should be, so we now ALLOCATE it ! as the correct size array. proj_fg_old%init = .false. initial_fg_read = .TRUE. IF ( initial_fg_read ) THEN ALLOCATE(dum2d(nx,ny)) ALLOCATE(xloc(nx,ny)) ALLOCATE(yloc(nx,ny)) ALLOCATE(xloc_t(nx,ny)) ALLOCATE(yloc_t(nx,ny)) initial_fg_read = .FALSE. END IF ! Loop over the different files that are needed for each ! time period. For example, this could include different places to get ! SST if it was not included in the file with the other meteorological ! variables. found_any_files = .FALSE. ! Determine if this processing is using initial conditions ! or lateral boundary conditions. find_tindex: DO time_index = 1,num_output_times IF (time_list(time_index) .EQ. time) EXIT find_tindex ENDDO find_tindex IF (time_index .LE. num_init_times) THEN files_per_time = num_files_per_single_time_init ELSE files_per_time = num_files_per_single_time_lbc ENDIF have_soilhgt = .false. diff_fg_files : DO diff_files_per_single_time = 1, & files_per_time + num_files_constants + add_one_for_lsm ! Get a UNIT number to OPEN. first_guess_unit : DO fg_unit = 10, 99 INQUIRE (UNIT=fg_unit, OPENED=already_opened) IF ( .NOT. already_opened) EXIT first_guess_unit END DO first_guess_unit ! F90 allows "as is" positioning, so make sure that does not ! slip in some how. REWIND ( UNIT=fg_unit ) ! Build the correct file name. This is inside a loop which will ! allow a search of all possible "roots" of the file name. file_name = ' ' // & ' ' IF ( diff_files_per_single_time .LE. num_files_constants ) THEN file_name_len = LEN_TRIM(& ADJUSTL(constants_names(diff_files_per_single_time))) file_name(1:file_name_len) = TRIM(ADJUSTL & (constants_names(diff_files_per_single_time))) constant_file = .TRUE. PRINT '(A,A)','Trying to OPEN constant file: ',TRIM(file_name) ELSEIF( (separate_lsm).AND. & (diff_files_per_single_time .LT. (num_files_constants+2))) THEN file_name_len= LEN_TRIM(ADJUSTL(lsm_prefix)) file_name(1:file_name_len) = & TRIM(ADJUSTL(lsm_prefix)) file_name(file_name_len+1:) = ':' // time(1:13) constant_file =.false. PRINT '(A,A)','Trying to OPEN LSM file: ',TRIM(file_name) ELSE IF (time_index .LE. num_init_times) THEN file_name_len = LEN_TRIM(ADJUSTL & (init_prefix & (diff_files_per_single_time-(num_files_constants+add_one_for_lsm)))) file_name(1:file_name_len) = & TRIM(ADJUSTL(init_prefix(& diff_files_per_single_time-(num_files_constants+add_one_for_lsm)))) file_name(file_name_len+1:) = ':' // time(1:13) constant_file = .FALSE. PRINT '(A,A)','Trying to OPEN init file: ',TRIM(file_name) ELSE file_name_len = LEN_TRIM(ADJUSTL & (lbc_prefix & (diff_files_per_single_time-(num_files_constants+add_one_for_lsm)))) file_name(1:file_name_len) = & TRIM(ADJUSTL(lbc_prefix(& diff_files_per_single_time-(num_files_constants+add_one_for_lsm)))) file_name(file_name_len+1:) = ':' // time(1:13) constant_file = .FALSE. PRINT '(A,A)','Trying to OPEN LBC file: ',TRIM(file_name) ENDIF END IF OPEN ( UNIT = fg_unit , & FILE = TRIM(file_name) , & STATUS = 'OLD' , & FORM = 'UNFORMATTED' , & IOSTAT = open_status ) IF ( open_status .NE. 0 ) THEN PRINT '(A,A)' , 'No file: ',TRIM(file_name) CYCLE diff_fg_files ELSE found_any_files = .TRUE. PRINT '(A,A)' , 'Successfully OPENed file: ', & TRIM(file_name) END IF ! This is the big loop over all of the data in this particular ! file. There is an outer loop over different files. landmask_fg_init = .false. this_fg_file_only : DO ! Read header data from FIRST GUESS file. READ ( UNIT=fg_unit , IOSTAT=eof_status) version ! Was this the end of the data set? IF ( eof_status .LT. 0 ) THEN PRINT '(A)','Found EOF in the FIRST GUESS file, exiting READ loop.' EXIT this_fg_file_only END IF ! Pull in another record of the header now that we know ! what the "format" will look like. This "hdate" gets the ! "hdate_input" stuff is so that the CHARACTER ! strings all are multiples of full bytes. ! We do not want to offend anyone. READ ( UNIT=fg_unit ) & hdate_input, xfcst, source, field_input, units_input, description_input, & level, idim, jdim, llflag hdate = hdate_input(1:19) field = field_input(1:8 ) units = units_input(1:24) description = description_input(1:46) ! Check the level flags to make sure we are not mixing/matching input ! vertical coordinate systems IF (level .GT. 200000.) THEN have_sfc = .true. ELSE IF ((level .LT. 200000.).AND.(level .GT. 200.)) THEN IF (have_noniso3d) THEN PRINT '(A,F7.0)', 'Apparent isobaric level found: ', level PRINT '(A8)', 'For field: ',field PRINT '(A)', 'However, we already had non-isobaric levels, and we cannot' PRINT '(A)', 'mix and match levels during one time period.' STOP 'iso_and_noniso_found' ELSE have_iso3d = .true. ENDIF ELSE IF (level .LE. 200) THEN IF (have_iso3d) THEN PRINT '(A,F7.0)', 'Apparent non-isobaric level found: ', level PRINT '(A8)', 'For field: ',field PRINT '(A)', 'However, we already had isobaric levels, and we cannot' PRINT '(A)', 'mix and match levels during one time period.' STOP 'iso_and_noniso_found' ELSE have_noniso3d = .true. ENDIF ENDIF ! We have the FG dimensions, so if we need the landmask_fg, we can allocate ! it and search for it IF (.NOT.landmask_fg_init) THEN ALLOCATE(landmask_fg(idim,jdim)) CALL get_gribprep_var('LANDSEA ', fg_unit, idim,jdim, & landmask_fg, found_var) IF (.NOT.found_var) THEN PRINT '(A)', 'No valid landmask found!' landmask_fg(:,:) = -1. make_landmask_fg = .true. ELSE make_landmask_fg = .false. WHERE(landmask_fg .GE. 0.5) landmask_fg = 1.0 WHERE(landmask_fg .LT. 0.5) landmask_fg = 0.0 print '(A,2F4.0)', 'Min/max landmask from FG: ', minval(landmask_fg), & maxval(landmask_fg) ENDIF landmask_fg_init = .true. landmask_fg = FLOAT(NINT(landmask_fg)) ! Prevents precision problems REWIND(fg_unit) CYCLE this_fg_file_only ENDIF ! We can test the llflag input as well. We can only deal with a couple of ! projections currently. Others will slide through, but be handled incorrectly ! without this test. The index values are defined in the library ! module_map_utils as public parameters, ! Current choices are: ! 0 - regular latitude/longitude projection (PROJ)LATLON) ! 1 - Mercator projection (PROJ_MERC) ! 3 - Lambert conformal projection (PROJ_LC) ! 5 - polar stereographic projection (PROJ_PS) ! IF ( (llflag .NE. PROJ_LATLON) .AND. & (llflag .NE. PROJ_MERC) .AND. & (llflag .NE. PROJ_LC) .AND. & (llflag .NE. PROJ_PS) ) THEN PRINT '(A,I3)','Unsupported projection code found in input file:',& llflag STOP 'Projection_type_unknown' END IF IF ( verbose ) THEN PRINT '(A)', 'INPUT DATA HEADER INFORMATION' PRINT '(A)', '-----------------------------' PRINT '(A,I2)','WRF pregrid version #: ', version PRINT '(A)', 'HDATE: '// hdate PRINT '(A,F7.2)', 'XFCST: ',xfcst PRINT '(A)', 'SOURCE: ' // source PRINT '(A)', 'FIELD: ' // field // ' in ' // units PRINT '(A)', 'DESCRIPTION: ' // description PRINT '(A,F10.1)','LEVEL: ',level PRINT '(A,I5,I5)', '(IDIM,JDIM): ', idim, jdim PRINT '(A,I2)', 'PROJECTION FLAG: ', llflag END IF ! Was this the time that we expected to find? We only care about this for first-guess ! fields that are not constant (the single input time would be assumed to be valid through ! the processing times requested). IF (( hdate(1:13).NE.time(1:13)).AND. & (.NOT. constant_file)) THEN PRINT '(A,A,A)','Expected to find date = ', & time(1:13),' in the input file.' PRINT '(A,A,A)','Instead, the date ',hdate(1:13), & ' was in the header.' STOP 'PROC_INGEST_GRIBPREP' ELSE IF ( constant_file ) THEN hdate(1:19) = time(1:19) IF ( verbose ) THEN PRINT '(A,A,A)','Constant file assumed valid for ', & time(1:13),'.' END IF END IF ! We know enough info from this READ statement to increment the ! counters for the 3d/2d data. This information is required on ! output in the header. If this is a "new" variable name, we ! increment either the "sfc" or the "up" counter. We can take care ! of duplicates at the end of this routine. IF ( (NINT(level) .EQ. 201300).OR.(NINT(level).EQ.200100 )) THEN found_name = .FALSE. find_name_sfc : DO max_fg_variables_index = 1 , max_fg_variables IF ( all_names_sfc(max_fg_variables_index)(1:8) .EQ. ' ' ) THEN all_names_sfc(max_fg_variables_index) = field // ' ' & // units // ' ' // description all_sources_sfc(max_fg_variables_index)(1:32) = source num2d_fg = num2d_fg + 1 found_name = .TRUE. EXIT find_name_sfc ELSE IF(field.EQ.all_names_sfc(max_fg_variables_index)(1:8))THEN found_name = .TRUE. EXIT find_name_sfc ELSE IF(field.NE.all_names_sfc(max_fg_variables_index)(1:8))THEN found_name = .FALSE. END IF END DO find_name_sfc ! If we did not find the name, something went wrong. IF ( .NOT. found_name ) THEN PRINT '(A)','Error in finding names for surface data.' STOP 'PROC_INGEST_FIRST_GUESS' END IF END IF ! Same search in the upper-air data. IF ( ( (NINT(level).EQ.100000).AND.(have_iso3d) ) .OR. & ( (NINT(level).EQ. 1) .AND. (have_noniso3d) ) ) THEN found_name = .FALSE. find_name_up : DO max_fg_variables_index = 1 , max_fg_variables IF ( all_names_up(max_fg_variables_index)(1:8) .EQ. ' ' ) THEN all_names_up(max_fg_variables_index) = field // ' ' & // units // ' ' // description all_sources_up(max_fg_variables_index)(1:32) = source num3d_fg = num3d_fg + 1 found_name = .TRUE. EXIT find_name_up ELSE IF(field.EQ. all_names_up(max_fg_variables_index)(1:8))THEN found_name = .TRUE. EXIT find_name_up ELSE IF(field.NE.all_names_up(max_fg_variables_index)(1:8))THEN found_name = .FALSE. END IF END DO find_name_up ! If we did not find the name, something went wrong. IF ( .NOT. found_name ) THEN PRINT '(A)','Error in finding names for upper-air data.' STOP 'search_num3d' END IF END IF ! Continue READs, if the first couple lines are there, ! the rest of the lines are supposed to be there as well, so no ! error traps. This methodology here is more flexible than the ! original MM5 REGRID, as it allows the input data ! set to be specified by either the center lat/lon or corner lat/lon. IF (llflag .EQ. PROJ_LATLON) THEN READ ( UNIT=fg_unit ) startloc,knownlat,knownlon, & deltalat, deltalon IF (knownlon .GT. 180) knownlon = knownlon - 360. IF (verbose) THEN PRINT '(A)', 'CYLINDRICAL EQUIDISTANT projection information:' PRINT '(A)', '-- Location of known lat/lon point: ' // startloc PRINT '(A,2F10.2)', '-- Known lat/lon: ',knownlat,knownlon PRINT '(A,2F10.2)', '-- Dlat/Dlon: ', deltalat, deltalon ENDIF ! The startlon value needs to be consistent with the way that the ! program produces the gridded longitude field for the output ! domain. The longitudes will be between (inclusively) ! -180 and +180. IF (startloc .NE. 'CENTER ') THEN ! We must be on one of the corners already, so this is easy startlat = knownlat startlon = knownlon ELSE ! We are on the center point, so compute the origin point startlat = knownlat - (FLOAT(jdim)-1.)*0.5*deltalat startlon = knownlon - (FLOAT(idim)-1.)*0.5*deltalon IF (verbose ) THEN PRINT '(A,2F10.2)', 'Computed lat/lon of (1,1)):', & startlat,startlon END IF END IF IF ( ( startlon .GT. 180 ) .AND. ( startlon .LE. 360 ) ) THEN startlon = startlon - 360 END IF ! Set up proj_fg_new CALL map_set(PROJ_LATLON,startlat,startlon,0.,deltalon,deltalat,& 0.,idim,jdim,proj_fg_new) ELSE IF ( llflag .EQ. PROJ_MERC) THEN READ(UNIT=fg_unit)startloc,lat1,lon1,deltax,deltay, & truelat1 IF (lon1 .gt. 180.) lon1 = lon1-360. IF (deltax .NE. deltay) THEN PRINT '(A)', 'WARNING: Mercator projection has different dx/dy!' PRINT '(A)', ' This is not yet supported, but I will assume' PRINT '(A)', ' that dy = the value set for dx and continue.' deltay = deltax ENDIF IF (verbose) THEN PRINT '(A)', 'MERCATOR projection information:' PRINT '(A)', '-- Location of known lat/lon point: ' // startloc PRINT '(A,2F10.2)', '-- Known lat/lon: ',lat1,lon1 PRINT '(A,2F10.2)', '-- Dx/Dy: ', deltax, deltay PRINT '(A,F10.2)', '-- Truelat: ', truelat1 ENDIF ! Initialize the proj_fg_new structure CALL map_set(PROJ_MERC,lat1,lon1,deltax*1000.,0.,truelat1,0., & idim, jdim, proj_fg_new) IF ( startloc .EQ. 'CENTER') THEN ! Proj structure assumes SW corner for lat1/lon2, so lets ! use ij_to_latlon to compute SW corner point and re-init ! the structure CALL ij_to_latlon(proj_fg_new,-FLOAT(idim+1)/2.+2.,-FLOAT(jdim+1)/2.+2., & lat1,lon1) CALL map_set(PROJ_MERC,lat1,lon1,deltax*1000.,0.,truelat1,0., & idim, jdim, proj_fg_new) ENDIF ELSE IF ( llflag .EQ. PROJ_LC) THEN READ ( UNIT=fg_unit ) startloc, lat1, lon1, & deltax, deltay, xlonc, truelat1, truelat2 IF (lon1 .GT. 180.) lon1 = lon1 - 360. IF (xlonc .GT. 180.) xlonc = xlonc - 360. IF (deltax .NE. deltay) THEN PRINT '(A)', 'WARNING: Different dx/dy values in LC projection.' PRINT '(A)', ' This is not yet supported, but I will continue' PRINT '(A)', ' by assuming correct value is dx.' deltay = deltax ENDIF IF (verbose) THEN PRINT '(A)', 'LAMBERT CONFORMAL projection information:' PRINT '(A)', '-- Location of known lat/lon point: ' // startloc PRINT '(A,2F10.2)', '-- Known lat/lon: ',lat1,lon1 PRINT '(A,2F10.2)', '-- Dx/Dy: ', deltax, deltay PRINT '(A,F10.2)', '-- Standard Lon: ', xlonc PRINT '(A,2F10.2)', '-- Truelats: ', truelat1, truelat2 ENDIF ! Initialize projection structure CALL map_set(PROJ_LC,lat1,lon1,deltax*1000.,xlonc,truelat1,truelat2, & idim, jdim, proj_fg_new) IF ( startloc .EQ. 'CENTER') THEN ! Proj structure assumes SW corner for lat1/lon2, so lets ! use ij_to_latlon to compute SW corner point and re-init ! the structure CALL ij_to_latlon(proj_fg_new,-FLOAT(idim+1)/2.+2.,-FLOAT(jdim+1)/2.+2., & lat1,lon1) CALL map_set(PROJ_LC,lat1,lon1,deltax*1000.,xlonc,truelat1,truelat2, & idim, jdim, proj_fg_new) ENDIF ELSE IF ( llflag .EQ. PROJ_PS ) THEN READ ( UNIT=fg_unit ) startloc, lat1, lon1, deltax, deltay, & xlonc, truelat1 IF (lon1 .GT. 180.) lon1 = lon1 - 360. IF (xlonc .GT. 180.) xlonc = xlonc - 360. IF (deltax .NE. deltay) THEN PRINT '(A)', 'WARNING: Different dx/dy values in PSprojection.' PRINT '(A)', ' This is not yet supported, but I will continue' PRINT '(A)', ' by assuming correct value is dx.' deltay = deltax ENDIF IF (verbose) THEN PRINT '(A)', 'POLAR STEREOGRAPHIC projection information:' PRINT '(A)', '-- Location of known lat/lon point: ' // startloc PRINT '(A,2F10.2)', '-- Known lat/lon: ',lat1,lon1 PRINT '(A,2F10.2)', '-- Dx/Dy: ', deltax, deltay PRINT '(A,F10.2)', '-- Standard Lon: ', xlonc PRINT '(A,2F10.2)', '-- Truelats: ', truelat1, truelat2 ENDIF ! Initialize projection structure CALL map_set(PROJ_PS,lat1,lon1,deltax*1000.,xlonc,truelat1,0., & idim, jdim, proj_fg_new) IF ( startloc .EQ. 'CENTER') THEN ! Proj structure assumes SW corner for lat1/lon2, so lets ! use ij_to_latlon to compute SW corner point and re-init ! the structure CALL ij_to_latlon(proj_fg_new,-FLOAT(idim+1)/2.+2.,-FLOAT(jdim+1)/2.+2., & lat1,lon1) CALL map_set(PROJ_PS,lat1,lon1,deltax*1000.,xlonc,truelat1,0., & idim, jdim, proj_fg_new) ENDIF ENDIF ! We now have a projection structure for the first guess data ! as well as the WRF domain. If this is not our first time ! through this point, we also know about the projection of ! the previous first guess field. Lets check to see if ! we need to do any interpolation at all by comparing our ! WRF domain projection to the FG projection. same_as_prev = .false. IF (proj_fg_old%init) THEN CALL compare_projections(proj_fg_new, proj_fg_old, same_as_prev) ENDIF IF (same_as_prev) THEN IF (verbose) & PRINT '(A)', 'Using same projection as previous field.' ELSE PRINT '(A)', 'Initializing new projection for this input field.' ! Check to see if interp is needed. If the input domain ! is identical to WRF domain, then we only need to extract ! data gridpoint for gridpoing CALL compare_projections(proj_fg_new, proj_wrf,no_interp_needed) IF (no_interp_needed) THEN PRINT '(A)', 'Input data already matches WRF domain!' DO j = 1, ny DO i = 1, nx xloc(i,j) = FLOAT(i) yloc(i,j) = FLOAT(j) ENDDO ENDDO DO j = 1, ny-1 DO i = 1, nx-1 xloc_t(i,j) = FLOAT(i) + 0.5 ! Arakawa C mass grid yloc_t(i,j) = FLOAT(j) + 0.5 ! Arakawa C mass grid ENDDO ENDDO ! The mass point grid does not actually used the top row ! or right column, so set these values to -1. as a flag ! for the interpolator xloc_t(nx,:) = -1. xloc_t(:,ny) = -1. yloc_t(:,ny) = -1. yloc_t(nx,:) = -1. ELSE DO j = 1, ny DO i = 1, nx !mp-nmm-BLS IF (proj_wrf%code .NE. PROJ_ROTLAT) THEN CALL latlon_to_ij(proj_fg_new,latitude(i,j,n_ind), & longitude(i,j,n_ind), xloc(i,j),yloc(i,j)) ELSE CALL latlon_to_ij(proj_fg_new,latitude(i,j,1), & longitude(i,j,1), xloc(i,j),yloc(i,j)) ENDIF !mp-nmm-BLS IF ( (NINT(xloc(i,j)).LT.0).OR.(NINT(xloc(i,j)).GT.idim+1) .OR. & (NINT(yloc(i,j)).LT.0).OR.(NINT(yloc(i,j)).GT.jdim+1) ) THEN PRINT '(A)', 'WRF MET Domain outside of FG domain!' PRINT '(A,2F5.1)', 'XLOC/YLOC(i,j) = ', xloc(i,j),yloc(i,j) STOP 'PROC_INGEST_GRIBPREP' ENDIF ! For non-wrapping data sets, limit the ! xloc/yloc IF (.NOT. proj_fg_new%cyclic) THEN xloc(i,j) = MAX(xloc(i,j),1.) xloc(i,j) = MIN(xloc(i,j),FLOAT(idim)) yloc(i,j) = MAX(yloc(i,j),1.) yloc(i,j) = MIN(yloc(i,j),FLOAT(jdim)) ENDIF ENDDO ENDDO print *, 'Min/Max XLOC (SRC) = ',MINVAL(xloc),MAXVAL(xloc) print *, 'Min/Max YLOC (SRC) = ',MINVAL(yloc),MAXVAL(yloc) DO j = 1, ny-1 DO i = 1,nx-1 CALL latlon_to_ij(proj_fg_new,latitude(i,j,t_ind), & longitude(i,j,t_ind), xloc_t(i,j),yloc_t(i,j)) ! We do not need to check for being on domain in this case ! because the used points of the t_grid are completely contained ! within the non-staggered grid, which was checked above. ENDDO ENDDO xloc_t(nx,:) = -1. xloc_t(:,ny) = -1. yloc_t(:,ny) = -1. yloc_t(nx,:) = -1. ENDIF ENDIF proj_fg_old = proj_fg_new ! We know how big the input array will be, so we ALLOCATE the ! correct size array. Since the array size can change, this is ALLOCATED ! and DEALLOCATEd for each READ statement. ALLOCATE (array(idim,jdim)) ! Read the 2d array. READ ( UNIT=fg_unit ) array ! Though -1e30 is a lovely flag value, to be sure, use zero instead. WHERE ( array .LT. -1.e20 ) array = 0 END WHERE ! Some quicky print statements of the input data, just the middle column, ! to bolster our confidence that this is being handled in a reasonable fashion. IF (verbose) THEN PRINT '(A)', 'Some values from the middle column:' PRINT '(8F10.1)',array(idim/2,1:jdim:5) END IF ! If this level is "higher up" than (physically above) the ptop defined in the ! NAMELIST, go back to the start of this READ loop. The data ! for this entry has been completely ingested, and would now be ! overwritten by the next entry (if this is the case where the data is higher ! up than the requested ptop). Note that this is after the ! EOF tests, so no infinite loop happens. IF ((level .LT. ptop_in_Pa ).AND.(have_iso3d)) THEN IF ( verbose) THEN PRINT '(A,F7.0,A,F7.0,A,A,A)','The NAMELIST ptop variable (',ptop_in_Pa, & ' Pa) requires that this level (',level,' Pa) of ', & field,' be omitted.' END IF ! We have to DEALLOCATE the array since we will skip that part of the ! loop with this CYCLE. DEALLOCATE (array) ! Get to the bottom of the loop and get the next data from the file. CYCLE this_fg_file_only END IF ! If this is the "LANDSEA" field and we need it, then we should ! have already gotten it. IF (field .EQ. 'LANDSEA ')THEN DEALLOCATE(array) CYCLE this_fg_file_only ENDIF ! Some models, like the RUC, originally have QVAPOR and ! VPTMP as state variables, but grib_prep derives T and RH ! from those. If these unnecessary vars are still in the ! file, skip them. IF ((field(1:6).EQ.'QVAPOR').OR.(field(1:5).EQ.'VPTMP'))THEN DEALLOCATE(array) CYCLE this_fg_file_only ENDIF ! If the background model has SOILCAT an/or VEGCAT variables ! provided, but we are not running LSM_HINTERP_METHOD = 0, ! then we can skip them. Otherwise, we keep them and ! process them using nearest neighbor (always!) IF ((field(1:6).EQ.'VEGCAT').OR.(field(1:7).EQ.'SOILCAT'))THEN IF (lsm_hinterp_method .NE. 0) THEN DEALLOCATE(array) CYCLE this_fg_file_only ENDIF ENDIF ! If we have already obtained a SOILHGT for this time, ! but we found it again, then skip IF (field(1:8) .EQ. 'SOILHGT ') THEN IF (have_soilhgt) THEN DEALLOCATE(array) CYCLE this_fg_file_only ELSE have_soilhgt = .true. ENDIF ENDIF ! Fix for AFWA AGRMET units IF (field(1:6).EQ.'CANOPY') THEN array = array*1000. ! Convert from meters to kg/m^2 field(1:8) = 'CANWAT ' ENDIF ! Fix RUC snow units IF (field(1:7) .EQ. 'SNOWRUC') THEN array = array*1000. ! convert from meters to k/m^2 field(1:8) = 'SNOW ' ENDIF ! Clean out the dum2d array before we use it to prevent problems ! with masked fields. BLS 12/12/2002 dum2d(:,:) = 0. ! Determine whether or not this is a masked field. IF ( ( field(1:2) .EQ. 'ST') .OR. & ( field(1:2) .EQ. 'SM') .OR. & ( field(1:2) .EQ. 'SW') .OR. & ( field(1:5) .EQ. 'SOILW') .OR. & ( field(1:5) .EQ. 'SOILM') .OR. & ( field(1:5) .EQ. 'SOILT') .OR. & ! wrfhelp: move SST and SEAICE out of masked interpolation ! ( field(1:3) .EQ. 'SEA' ) .OR. & ! ( field(1:3) .EQ. 'SST' ).OR. & ! ( field(1:8) .EQ. 'SKINTEMP' ).OR. & ( field(1:4) .EQ. 'SNOW' ) .OR. & ( field(1:6) .EQ. 'CANWAT' ) .OR. & ( field(1:8) .EQ. 'WEASD ' ))THEN ! Hopefully, we already have a true land/water mask ! from this FG source, but if not, create a dummy ! with negative values so the masked interpolater ! knows to try and guess whether it is a land or ! water value IF (.NOT. ALLOCATED(landmask_fg)) THEN ALLOCATE(landmask_fg(idim,jdim)) landmask_fg = -1. ENDIF IF (field(1:5).EQ. 'SNOWC') THEN min_valid = 0. max_valid = 1. default_val = 0. val_mask = 1. ELSE IF (field(1:3).EQ. 'SEA') THEN min_valid = 0. max_valid = 1. default_val = 0. val_mask = 0. ELSE IF((field(1:2).EQ.'SM').OR.(field(1:5).EQ.'SOILM').OR.& (field(1:2).EQ.'SW').OR.(field(1:5).EQ.'SOILW')) THEN ! Changed min_valid from 0.001 min_valid = 0.0000 max_valid = 0.9999 default_val = 0.30 val_mask = 1.0 ELSE IF((field(1:2).EQ.'ST').OR.(field(1:5).EQ.'SOILT'))THEN min_valid = 100. max_valid = 400. default_val = 280. val_mask = 1.0 ELSE IF(field(1:8).EQ.'SKINTEMP')THEN min_valid = 100. max_valid = 400. default_val = 280. val_mask = 2.0 ELSE IF(field(1:3).EQ.'SST')THEN min_valid = 270. max_valid = 350. default_val = 290. val_mask = 0. ELSE IF(field(1:6).EQ.'CANWAT')THEN min_valid = 0. max_valid = 10.0 default_val = 0.06 val_mask = 1. ELSE IF(field(1:8).EQ.'SOILHGT ')THEN min_valid = -1000. max_valid = 10000. default_val = 0.00 val_mask = 1. ELSE min_valid = 0. max_valid = 1000. default_val = 0. val_mask = 2.0 ENDIF PRINT '(2A)', 'Doing masked interpolation for ', field ! If this is a "flag" field, like snow cover or sea ice, we ! should do the nearest neighbor method. Otherwise, use ! the linear method. IF ((field(1:6) .EQ. 'SEAICE').OR.(field(1:8).EQ.'SNOWCOVR'))THEN CALL interpolate_masked_val(idim, jdim, landmask_fg, array, & nx-1, ny-1, landmask_wrf(1:nx-1,1:ny-1), & dum2d(1:nx-1,1:ny-1), xloc_t(1:nx-1,1:ny-1), & yloc_t(1:nx-1,1:ny-1), make_landmask_fg, & min_valid, max_valid, & default_val, val_mask, METHOD_NEAREST) ELSE CALL interpolate_masked_val(idim, jdim, landmask_fg, array, & nx-1, ny-1, landmask_wrf(1:nx-1,1:ny-1), & dum2d(1:nx-1,1:ny-1), xloc_t(1:nx-1,1:ny-1), & yloc_t(1:nx-1,1:ny-1), make_landmask_fg, & min_valid, max_valid, & default_val, val_mask, interp_method_lsm) ENDIF ! Fill right row and upper column just for cleanliness. dum2d(:,ny) = dum2d(:,ny-1) dum2d(nx,:) = dum2d(nx-1,:) ELSEIF ((field(1:6).EQ.'VEGCAT').OR.(field(1:7).EQ.'SOILCAT'))THEN ! We are running lsm_hinterp_method 0 and background model has ! landuse and/or soil category data that needs to be passed ! through. We must put the values from the background onto the ! mass points, and always use nearest neighbor! ! Make sure we do not have any fractional values array = FLOAT(NINT(array)) val_mask = 1. ! AGRMET data was bitmasked in GRIB to only have land points, ! so the water points end up being zero. Correct that here ! by setting it to the USGS wate value (16) IF (field(1:8).EQ. 'VEGCAT ') THEN WHERE(array .EQ. 0.) array = 16. ! Make sure landmask_fg is consistent landmask_fg(:,:) = 1. WHERE(array .EQ. 16) landmask_fg = 0. min_valid = 1. max_valid = 24. default_val = 9. ENDIF ! Same issue for soil cateogry, which is supposed to use 14 ! for water points IF (field(1:8).EQ.'SOILCAT ') THEN WHERE(array .EQ. 0.) array = 14 ! Make sure landmask_fg is consistent landmask_fg(:,:) = 1. WHERE(array .EQ. 14) landmask_fg = 0. min_valid = 1. max_valid = 16. default_val = 8. ENDIF PRINT *,'Passing the background ',field(1:8), ' through to output.' CALL interpolate_masked_val(idim, jdim, landmask_fg, array, & nx-1, ny-1, landmask_wrf(1:nx-1,1:ny-1), & dum2d(1:nx-1,1:ny-1), xloc_t(1:nx-1,1:ny-1), & yloc_t(1:nx-1,1:ny-1), .false., & min_valid, max_valid, & default_val, val_mask, METHOD_NEAREST) ELSE ! This is presumed to be a standard,non-masked meteorological field, ! for which we interpolate to the non-staggered grid. ! In the case of hinterp_method = 2, it is possible for the ! interpolated result to have values less than the original ! minimum value and/or values greater than the original ! maximum. This may cause a problem with relative humidity, ! so before we interpolate, capture the min/max values ! and use those to reset the output values IF (field(1:8) .EQ. 'RH ') THEN min_rh_in = MINVAL(array) max_rh_in = MAXVAL(array) ENDIF CALL interpolate_standard(idim, jdim, array, & nx, ny, xloc, yloc, interp_method, & dum2d) ! If we were RH and hinterp_method gt 1, then apply correction for ! min and max values IF ((field(1:8) .EQ. 'RH ').AND.(interp_method .GT. METHOD_LINEAR))THEN PRINT *, 'Limiting min/max RH to original input min/max ', min_rh_in,max_rh_in WHERE (dum2d .GT. max_rh_in) dum2d = max_rh_in WHERE (dum2d .LT. min_rh_in) dum2d = min_rh_in ENDIF ! If we just processed the HGT variable at the surface, then allocate ! and generate the soilhgt_bg data array using the same masked interpolation ! procedures as the LSM IF ((field(1:8).EQ.'HGT ').AND.(NINT(level).EQ.200100))THEN IF (.NOT.ALLOCATED(soilhgt_bg)) ALLOCATE(soilhgt_bg(nx,ny)) soilhgt_bg(:,:)=-9999. print *, 'Interpolating soil height with method: ', interp_method_lsm CALL interpolate_masked_val(idim,jdim,landmask_fg,array, & nx-1, ny-1, landmask_wrf(1:nx-1,1:ny-1), & soilhgt_bg(1:nx-1,1:ny-1), xloc_t(1:nx-1,1:ny-1), & yloc_t(1:nx-1,1:ny-1), make_landmask_fg, & -1000., 20000., & 1., 1., interp_method_lsm) soilhgt_bg(:,ny) =soilhgt_bg(:,ny-1) soilhgt_bg(nx,:) =soilhgt_bg(nx-1,:) ENDIF ENDIF ! Ensure positive definite where appropriate SELECT CASE(field(1:8)) CASE('QLIQUID ') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('QICE ') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('QRAIN ') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('QGRAUPEL') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('QSNOW ') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('QNI ') WHERE(dum2d .LT. 0.) dum2d = 0. CASE('RH ') WHERE(dum2d .LT. 1.) dum2d = 1. CASE('WEASD ') WHERE(dum2d .LT. 0.) dum2d = 0. field(1:8) = 'SNOW ' CASE('CANWAT ') WHERE(dum2d .LT. 0.) dum2d = 0. END SELECT IF ( verbose ) THEN PRINT '(//,A,1X,I8)',field,NINT(level) DO j=ny,1,-10 PRINT '(10f9.1)',dum2d(1:nx:10,j) END DO END IF ! Store this 2d data on the model grid. CALL proc_grid_store ( field // ' ' // units // ' ' // description , level , dum2d , nx , ny ) ! Since the input data can have different sizes for each level ! (such as with the surface compared to the upper-air levels), ! we need to DEALLOCATE and ALLOCATE the input array for each ! slice. DEALLOCATE (array) END DO this_fg_file_only IF (ALLOCATED(landmask_fg)) DEALLOCATE(landmask_fg) ! CLOSE this file. We can reuse this UNIT again if required. CLOSE (UNIT=fg_unit) IF (ALLOCATED(landmask_fg)) DEALLOCATE(landmask_fg) PRINT '(A,A/)' , 'CLOSE file: ',TRIM(file_name) END DO diff_fg_files ! Check to see if we have at least OPENed up a single file. This is ! to provide a gentle and controlled exit from the program. IF ( .NOT. found_any_files ) THEN PRINT '(A)','Could not find any files to OPEN.' PRINT '(A)','You may want to check the filenames and the dates in the NAMELIST file.' data_processed = .false. ELSE data_processed = .true. ! We have ingested all of the first guess data for this time. ! We can now go through the list of stored names and determine ! how many 3d/2d variables there are. The only replications ! are if there is a field in the 3d and in the 2d array (such ! as temperature). We remove the name from the 2d array and ! decrement the 2d counter. check_2d : DO fg_variables_sfc_index = 1 , max_fg_variables ! If this is blank, EXIT this loop. IF ( all_names_sfc(fg_variables_sfc_index)(1:8) .EQ. ' ' ) THEN EXIT check_2d END IF check_3d : DO fg_variables_up_index = 1 , max_fg_variables ! If this is blank, EXIT this loop. IF ( all_names_up(fg_variables_up_index)(1:8) .EQ. ' ' ) THEN EXIT check_3d END IF ! Are the fields the same? IF ( all_names_up(fg_variables_up_index) .EQ. all_names_sfc(fg_variables_sfc_index) ) THEN IF ( verbose ) THEN PRINT '(A,A,A)','Removing ',all_names_sfc(fg_variables_sfc_index)(1:8),' from the 2d list.' END IF all_names_sfc(fg_variables_sfc_index)(1:8) = 'XXXXXXXX' num2d_fg = num2d_fg - 1 END IF END DO check_3d END DO check_2d PRINT '(A,I4,A,I4,A)','Found ',num3d_fg,' 3d arrays and ',num2d_fg,' 2d arrays for this time.' ! Figure out how many levels there are in the 3d data. This can ! be different for different variables, so pick the maximum of all of the numbers. current_var => head find_max_levels : DO WHILE ( ASSOCIATED(current_var) ) number_of_original_levels = MAX ( number_of_original_levels , current_var%num_levels ) current_var => current_var%next_variable END DO find_max_levels IF ( verbose ) THEN PRINT '(A,I4,A)','Found a maximum of ',number_of_original_levels,' for the 3d data.' END IF ! Zap the ALLOCATEd space for the model domain values. DEALLOCATE(dum2d) DEALLOCATE(xloc) DEALLOCATE(yloc) IF (ALLOCATED(xloc_t)) DEALLOCATE(xloc_t) IF (ALLOCATED(yloc_t)) DEALLOCATE(yloc_t) ! Convert the linked list to a more convenient array-based representation ! (this used to be in the hinterp.F main driver) IF ( .NOT. have_noniso3d) THEN CALL proc_list_to_array ELSE CALL proc_list_to_array_noniso ENDIF ! Free up some space CALL proc_zap_space_list ! Update the domain metadata dom_meta%zdim = number_of_original_levels + num_new_levels ! Rotate the winds if needed IF (.NOT. no_interp_needed) THEN ! The FG and WRF grids are not the same projection, so we ! might need to rotate them. The rotation involves ! converting the FG grid-relative winds to true winds ! (if they are not already), and then rotating the true ! winds to the WRF grid-relative IF ( (proj_fg_new%code .NE. PROJ_LATLON).AND. & (proj_fg_new%code .NE. PROJ_MERC) ) THEN rotate_fg_winds = .true. ELSE rotate_fg_winds = .false. ENDIF ! If the WRF projection is mercator, then true winds ! are also grid-relative IF (proj_wrf%code .NE. PROJ_MERC) THEN rotate_wrf_winds = .true. ELSE rotate_wrf_winds = .false. ENDIF IF ((rotate_fg_winds).OR.(rotate_wrf_winds)) THEN ALLOCATE(u_temp(nx,ny)) ALLOCATE(v_temp(nx,ny)) ! Search through the 3d data array to get the index ! of u and v DO i = 1 , fg_variables_up_index IF ( all_names_up(i)(1:8) .EQ. 'U ' ) THEN u_index = i ELSE IF ( all_names_up(i)(1:8) .EQ. 'V ' ) THEN v_index = i END IF END DO !mp-BLS ! For ROTLAT (NMM) stuff, want to use angles read in from ! static.wrfsi.rotlat file !mp-BLS DO k = 1, dom_meta%zdim DO j = 1, ny DO i = 1, nx IF (rotate_fg_winds) THEN !mp-BLS IF (proj_wrf%code .NE. PROJ_ROTLAT) THEN CALL gridwind_to_truewind(longitude(i,j,n_ind), & proj_fg_new,all_3d(i,j,k,u_index), & all_3d(i,j,k,v_index),u_temp(i,j),v_temp(i,j)) ELSE CALL gridwind_to_truewind(longitude(i,j,2), & proj_fg_new,all_3d(i,j,k,u_index), & all_3d(i,j,k,v_index),u_temp(i,j),v_temp(i,j)) ENDIF all_3d(i,j,k,u_index) = u_temp(i,j) all_3d(i,j,k,v_index) = v_temp(i,j) ENDIF IF (rotate_wrf_winds) THEN IF (proj_wrf%code .NE. PROJ_ROTLAT) THEN CALL truewind_to_gridwind(longitude(i,j,n_ind), & proj_wrf,all_3d(i,j,k,u_index), & all_3d(i,j,k,v_index),u_temp(i,j),v_temp(i,j)) ELSE u_temp(i,j)=all_3d(i,j,k,u_index)*cosalp(i,j) - & all_3d(i,j,k,v_index)*sinalp(i,j) v_temp(i,j)=all_3d(i,j,k,u_index)*sinalp(i,j) + & all_3d(i,j,k,v_index)*cosalp(i,j) ENDIF all_3d(i,j,k,u_index) = u_temp(i,j) all_3d(i,j,k,v_index) = v_temp(i,j) ENDIF ENDDO ENDDO ENDDO DEALLOCATE(u_temp) DEALLOCATE(v_temp) ENDIF ENDIF ENDIF END SUBROUTINE ingest_gribprep !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE proc_grid_store ( name , level , field , ix , jx ) ! Stores a 2D array into the linked list IMPLICIT NONE CHARACTER (LEN=80) , INTENT(IN) :: name REAL , INTENT(IN) :: level INTEGER , INTENT(IN) :: ix,jx REAL , DIMENSION(ix,jx) , INTENT(IN) :: field LOGICAL :: bad = .FALSE. initialized_TF : IF ( .NOT. initialized ) THEN IF ( verbose ) THEN PRINT '(A,A)', 'Uninitialized ',name END IF initialized = .TRUE. ! Initialize the head of the linked list. head%num_levels=0 NULLIFY(head%next_variable) NULLIFY(head%slab) ! We are going to index around with the current pointer. current_var => head ! Store the name of this variable. Since this is the first ! one, there is no search to see what order they go in. current_var%name = name ! ALLOCATE space for a new data level. This does not contain ! the gridded 2d data itself. ALLOCATE (new_slab) NULLIFY (new_slab%next_slab) current_var%slab => new_slab ! Store which level this data is valid for. current_var%slab%level = level ! Allocate the space needed to store the 2d data. ALLOCATE (new_field(ix,jx)) current_var%slab%field => new_field ! Store the input gridded 2d data in this space. Again, ! since this is the first store, we do not need to search ! for the correct vertical level. current_var%slab%field = field ! Initialize the number of 2d fields stored. The first one ! is always the easy one. current_var%num_levels = 1 ELSE IF ( initialized ) THEN IF ( verbose ) THEN PRINT '(A,A)', 'Initialized ',name END IF ! Check to see where this name can fit in. We do not ! need to order these alphabetically. Typically, there ! are going to be less than 30 different names, so a ! linear search is not a problem. If this name does not ! already exist, then ALLOCATE space for this new variable. current_var => head found = .FALSE. find_right_name : DO WHILE ( ASSOCIATED(current_var) ) IF ( current_var%name(1:8) .EQ. name(1:8) ) THEN found = .TRUE. EXIT find_right_name END IF old_var => current_var current_var => current_var%next_variable END DO find_right_name ! If we need to add this name into the linked list, we ! now know where to add it in (this is always just at ! the end of the list, as explained above). new_name_to_make : IF ( .NOT. found ) THEN IF ( verbose ) THEN PRINT '(A,A)', 'Name not found: ',name END IF ! ALLOCATE space for this new variable. ALLOCATE(new_var) new_var%num_levels = 0 NULLIFY (new_var%next_variable) NULLIFY (new_var%slab) ! Include the new variable in the linked list (the one ! before this has to know how to get to him). old_var%next_variable => new_var ! Put the current pointer on the new variable. current_var => new_var ! Store the name of this variable. current_var%name = name ELSE IF ( found ) THEN IF ( verbose ) THEN PRINT '(A,A)', 'Name found: ',name END IF END IF new_name_to_make ! Increment the levels counter. current_var%num_levels = current_var%num_levels + 1 ! Find where to stick the vertical level. This is ordered such that the ! first level is closest to the ground (i.e., the largest pressure). ! If there is only one level, then no search. ! This section has been modified to account for non-isobaric levels, which ! are assumed to be numbered with lower values being closer to the surface ! (opposite of isobaric levels) first_level_TF : IF ( current_var%num_levels .EQ. 1 ) THEN IF ( verbose ) THEN PRINT '(A,A,1X,I8)', 'First level for this variable: ',name,NINT(level) END IF ! Make the space for the new level type. Have the pointer ! index to it. ALLOCATE (new_slab) NULLIFY (new_slab%next_slab) current_var%slab => new_slab ! Now, make space for the 2d data. ALLOCATE (new_field(ix,jx)) current_var%slab%field => new_field ! Load up the 2d data into the holder. current_var%slab%level = level current_var%slab%field = field ! We need to sort the levels since we already have at least ! one of them. ELSE IF ( verbose ) THEN PRINT '(A,A,1X,I8)', 'Other levels exist for this variable: ',name,NINT(level) END IF found = .FALSE. current_slab => current_var%slab old_slab => current_var%slab ! Does the new level go before all levels? This is only checked ! once (outside of the following DO loop). bigger_than_first_level: IF ( ((have_iso3d).AND.(level.GT.current_slab%level)).OR. & ( (have_noniso3d).AND.(level.LT.current_slab%level) ) ) THEN IF ( verbose ) THEN PRINT '(A,A,1X,2I8)', 'This level goes before the others: ',& name,NINT(level),NINT(current_slab%level) END IF ! ALLOCATE space for the slab type and the 2d data. ALLOCATE (new_slab) NULLIFY (new_slab%next_slab) ! Make space for the 2d data, load up the space. ALLOCATE (new_field(ix,jx)) new_slab%field => new_field ! Have the new slab point to the rest of the list. new_slab%next_slab => current_slab ! Have the top of the variable list point to this slab. current_var%slab => new_slab ! Fill in all of the data that we have for this slab. new_slab%level = level new_field = field ! foo this does not work on DEC, so using the above kludge ! new_slab%field = field ! The data has been placed in the linked list successfully, so we ! set a flag so that we do not go into the following loop. found = .TRUE. END IF bigger_than_first_level ! Either the new level goes in between two levels or at the end of ! the list. Check these two choices. A DO loop is required as we ! traverse all of the way to the bottom. We do not go in this ! loop if we have successfully already placed the level above. find_right_level : DO WHILE ( ASSOCIATED(current_slab) .AND. ( .NOT. found) ) ! The new level goes at the end of the list. vertical_search : IF ( ( .NOT. ASSOCIATED (current_slab%next_slab) ) .AND. & ( ((level .LT. current_slab%level).AND.(have_iso3d)).OR. & ((level .GT. current_slab%level).AND.(have_noniso3d)) ) ) THEN IF ( verbose ) THEN PRINT '(A,A,1X,2I8)', 'This level goes after all the others: ',name,NINT(level),NINT(current_slab%level) END IF ! ALLOCATE space for the slab type and the 2d data. ALLOCATE (new_slab) NULLIFY (new_slab%next_slab) ! Make space for the 2d data, load up the space. ALLOCATE (new_field(ix,jx)) new_slab%field => new_field ! Fill in all of the data that we have for this slab. new_slab%level = level new_field = field ! new_slab%field = field ! Have the previous slab in the list point to this slab. current_slab%next_slab => new_slab ! EXIT the DO loop, the data has been placed in the linked ! list successfully. found = .TRUE. EXIT find_right_level ! The new level goes between two levels. ELSE IF ( ( level .GT. current_slab%level ) .AND. & ( level .LT. old_slab%level ) ) THEN IF ( verbose) THEN PRINT '(A,A,1X,3I8)', 'This level goes between two other levels: ', & name,NINT(level),NINT(old_slab%level),NINT(current_slab%level) END IF ! ALLOCATE space for the slab type and the 2d data. ALLOCATE (new_slab) NULLIFY (new_slab%next_slab) ! Make space for the 2d data, load up the space. ALLOCATE (new_field(ix,jx)) new_slab%field => new_field ! Have the new slab point to the rest of the list. new_slab%next_slab => current_slab ! Have the previous slab in the list point to this slab. old_slab%next_slab => new_slab ! Fill in all of the data that we have for this slab. new_slab%level = level new_field = field ! foo this does not work on DEC, so using the above kludge ! new_slab%field = field ! EXIT the DO loop, the data has been placed in the linked ! list successfully. found = .TRUE. EXIT find_right_level ! Oops, we found a duplicate. This is the same level with the ! same name, and since this is a single file, at the same time. ! This could happen, for example, if there was an SST data segment ! in the first guess file, as well as in a separate SST file. In any ! case, disregard the second occurrence. ELSE IF ( ABS(level-old_slab%level) .LT. 1 ) THEN ! Inform folks of the situation. PRINT '(A,A,A,F10.1,A)','A duplicate variable/level was detected: ',name(1:8),' at level ',level,' Pa.' PRINT '(A)','The first field/level data values are retained, & &and all subsequent occurrences are disregarded.' ! Decrement the levels counter to get rid of this "bad" value. current_var%num_levels = current_var%num_levels - 1 ! EXIT the loop - effectively, leave this routine. Neither ! a new slab (for the new level) nor a new 2d field (for the ! actual (IX,JX) field) is generated. We are exiting this ! routine without any changes to the linked list. bad = .TRUE. EXIT find_right_level END IF vertical_search IF ( verbose ) THEN PRINT '(A)', 'Have not yet found where the level goes.' END IF ! We have not found where the new level goes, so ! increment the pointers. old_slab => current_slab current_slab => current_slab%next_slab ! Just in case the duplicate is at the end of the linked list, we must check ! against this level now. In any case, disregard the second occurrence. IF ( ABS(level-old_slab%level) .LT. 1 ) THEN ! Inform folks of the situation. PRINT '(A,A,A,F10.1,A)','A duplicate variable/level was detected: ',name(1:8),' at level ',level,' Pa.' PRINT '(A)','The first field/level data values are retained, & &and all subsequent occurrences are disregarded.' ! Decrement the levels counter to get rid of this "bad" value. current_var%num_levels = current_var%num_levels - 1 ! EXIT the loop - effectively, leave this routine. Neither ! a new slab (for the new level) nor a new 2d field (for the ! actual (JX,IX) field) is generated. We are exiting this ! routine without any changes to the linked list. bad = .TRUE. EXIT find_right_level END IF END DO find_right_level IF ( ( .NOT. found ) .AND. ( .NOT. bad ) ) THEN PRINT '(A)', 'Error in finding correct level in data store.' STOP 'data_level_store' END IF END IF first_level_TF END IF initialized_TF ! Re-initialize the "bad" value flag to .FALSE. We want to ! assume that the incoming data is .NOT. bad. The bad flag is ! only set when we find a duplicate of the same variable ! at the same vertical level for the same valid time. bad = .FALSE. END SUBROUTINE proc_grid_store !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE proc_zap_space_list IMPLICIT NONE ! Point to the head with an incrementing pointer, and ! one that will stay put. current_var => head new_var => head ! WHILE the pointer is at something that is ASSOCIATED, ! drop down into its levels. each_variable : DO WHILE ( ASSOCIATED ( new_var ) ) ! Hold onto where we are pointing since this is not a ! doubly linked list. old_var => new_var ! Initialize the various levels pointer. new_slab => new_var%slab ! Check each one of the levels. each_level : DO WHILE ( ASSOCIATED ( new_slab ) ) ! Again, hold onto where we are as this is not ! doubly linked. old_slab => new_slab ! Check for the next level. new_slab => new_slab%next_slab ! DEALLOCATE what we can. DEALLOCATE ( old_slab%field ) DEALLOCATE ( old_slab ) END DO each_level ! Increment to the next variable. new_var => new_var%next_variable ! If this is not the head of the list (which is a ! TARGET - not ALLOCATABLE), zap this variable. IF ( current_var%name .NE. old_var%name ) THEN IF ( ASSOCIATED(old_var) ) THEN DEALLOCATE ( old_var ) ENDIF END IF END DO each_variable ! Clean up the head to an unused an pristine state. head%name = ' ' head%num_levels = 0 NULLIFY ( head%next_variable ) NULLIFY ( head%slab ) ! Finally, set the linked list initialized flag to FALSE. initialized = .FALSE. END SUBROUTINE proc_zap_space_list !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE proc_list_to_array ! ! This routine takes the linked list of 2D slabs combines ! them appropriately into 2D and 3D arrays to populate ! the all_3d and all_2d structures defined in the ! module_gridded_data. ! ! HISTORY ! ======= ! Nov 2000 - Brent L. Shaw, NOAA/FSL ! (Adapted from NCAR MM5v3 REGRIDDER Routine) ! Apr 2001 - Brent L. Shaw, NOAA/FSL ! Mods to support improved code structure of hinterp IMPLICIT NONE ! These variables are used for the storage and vertical ! interpolation of the 3d data. INTEGER :: fg_vertical_index , new_levels_index REAL , DIMENSION(nx,ny) :: slab_bottom, slab_top, slab_middle REAL :: pressure_bottom, pressure_top, pressure_middle LOGICAL :: vertically_interpolate = .FALSE. LOGICAL :: any_surface_data ! Flags used for searching for vertical coordinate limiting cases. INTEGER :: num_same_press , nnl LOGICAL :: need_ptop_level ! Initialize the two name arrays (both for surface and upper-air). all_names_sfc = ' ' all_names_up = ' ' fg_variables_sfc_index = 0 fg_variables_up_index = 0 new_levels_index = 0 ! Find out how many levels are in the 3d data. It is possible ! that some fields do not have as many levels as other fields ! (such as RH, there may not be a surface value). Loop through ! all of the variable names and find the maximum of the number ! of levels. any_surface_data = .FALSE. current_var => head number_of_original_levels = 0 find_max_levels : DO WHILE ( ASSOCIATED(current_var) ) number_of_original_levels = MAX ( number_of_original_levels , current_var%num_levels ) IF ( verbose) THEN PRINT '(A,I8,A,A)','Found ',current_var%num_levels,' levels of data for variable ',current_var%name END IF IF ( ( current_var%num_levels .GT. 1 ) .AND. ( .NOT. any_surface_data ) ) THEN IF ( ABS(current_var%slab%level-200100.) .LT. 1 ) THEN any_surface_data = .TRUE. END IF END IF current_var => current_var%next_variable END DO find_max_levels IF ( verbose ) THEN PRINT '(A,I8,A)','Found ',number_of_original_levels,' as maximum number of original vertical levels.' END IF !!Next section commented out.. MOD BLS-20030306 ! If there is no surface data, we need to let the folks know. ! IF ( .NOT. any_surface_data ) THEN ! PRINT '(A)',' ' ! PRINT '(A)','****************************************************************************' ! PRINT '(A)','****************************************************************************' ! PRINT '(A)','There is no surface data in your data set.' ! PRINT '(A)','The 1000 hPa level data is being used as the surface data for all fields.' ! PRINT '(A)','****************************************************************************' ! PRINT '(A)','****************************************************************************' ! PRINT '(A)',' ' ! END IF ! Correct the pressure level data for a few things. 1) If we have a ! level and the same level is requested as a new level, decrement the ! number of original levels; 2) make sure that ptop is one of the ! levels in the final data (it cannot be a new level). current_var => head num_same_press = 0 need_ptop_level = .TRUE. ! Loop over all of the variables. find_p_peculiarities : DO WHILE ( ASSOCIATED(current_var) ) ! Find the variable which has the most vertical levels. If there are ! several with the same number, we just find the first one. IF ( number_of_original_levels .EQ. current_var%num_levels ) THEN ! Loop over each first-guess level for this variable. current_slab => current_var%slab first_guess_p_loop : DO WHILE ( ASSOCIATED ( current_slab ) ) ! Loop over each requested pressure level from the namelist. find_dups : DO nnl = 1 , num_new_levels ! Do we have a requested level that already exists? IF ( ABS ( current_slab%level - new_levels(nnl) ) .LT. 1 ) THEN num_same_press = num_same_press + 1 IF ( verbose ) THEN PRINT '(A,I8,A)','Found duplicate pressure, one from namelist and the other & &from the first guess ',NINT(current_slab%level),' Pa.' END IF END IF END DO find_dups ! Another thing, did we find a first-guess level for ptop? IF ( ABS ( current_slab%level - ptop_in_Pa ) .LT. 1 ) THEN need_ptop_level = .FALSE. IF ( verbose ) THEN PRINT '(A)','Found first-guess level for ptop.' END IF END IF current_slab => current_slab%next_slab END DO first_guess_p_loop EXIT find_p_peculiarities ELSE current_var => current_var%next_variable END IF END DO find_p_peculiarities ! If there is no surface data for the 3d arrays, we need to add an extra level. ! Next 3 lines commented out MOD BLS-20030306 !IF ( .NOT. any_surface_data ) THEN ! number_of_original_levels = number_of_original_levels + 1 !END IF ! If there are any replicated levels, reduce the number of original levels. number_of_original_levels = number_of_original_levels - num_same_press ! If ptop is not in the first-guess levels, we need to stop. IF ( need_ptop_level ) THEN PRINT '(A,I8,A)','The requested PTOP variable in the namelist is ',NINT(ptop_in_Pa),' Pa.' PRINT '(A)','This must be in the first-guess field levels.' STOP 'change_ptop' END IF ! The vertical array of pressure is ALLOCATEd. This is the sum ! of the data from the first guess file and the new levels from the ! NAMELIST. Note that the surface level is included in the count ! from the first guess file. ALLOCATE (output_levels(number_of_original_levels+num_new_levels) ) ! ALLOCATE the space required for the 3d data. One array for ! each variable. The vertical dimension is the sum of the ! input levels (from the first guess fields themselves) and the ! requested new levels (from the NAMELIST). ALLOCATE ( all_3d(nx,ny,number_of_original_levels+num_new_levels,num3d_fg) ) ! ALLOCATE the space for the 2d arrays. !ALLOCATE ( all_2d(nx,ny,num2d_fg*3) ) ! COMMENTED OUT BLS-20030306 ALLOCATE ( all_2d(nx,ny,num2d_fg) ) ! REPLACED ABOVE BLS-20030306 ! Initialize some of the required counters to zero. fg_variables_sfc_index = 0 fg_variables_up_index = 0 fg_vertical_index = 0 ! We are going to loop over all of the variables, so start ! by pointing to the first one. We loop until we are no ! longer pointing to an ASSOCIATEd variable. current_var => head search_all_variables : DO WHILE ( ASSOCIATED(current_var) ) ! If there is only a single level, this is 2d data. As ! is seen below, the 2d data is a piece of cake to store. twod_or_threed : IF ( current_var%num_levels .EQ. 1 ) THEN ! Find where to stick this data. Fill an array of the ! name, then store the 2d data in the same location of ! the ALLOCATEd 2d array. Increment the index of found ! data, then load the data, then save the name of the data. ! This is easy. fg_variables_sfc_index = fg_variables_sfc_index + 1 all_2d(:,:,fg_variables_sfc_index) = current_var%slab%field all_names_sfc(fg_variables_sfc_index) = current_var%name IF ( verbose ) THEN PRINT '(A,I8,A,A)','Store data in arrays: Found 2d data index #',fg_variables_sfc_index, & ' named ',all_names_sfc(fg_variables_sfc_index)(1:8) END IF ! If there is more than a single level, this is a 3d array. As ! is judged by the following rest of the IF test, this is a bit ! more involved. There are vertical levels to find, possible ! interpolations and maybe missing surface values. Fasten your ! seat belt. ELSE IF ( current_var%num_levels .GT. 1 ) THEN ! Increment the index of found upper-air data. This is a ! count of the total number of 3d fields in the data set. fg_variables_up_index = fg_variables_up_index + 1 ! Save the name of this 3d array. all_names_up(fg_variables_up_index) = current_var%name ! Point to the first level of this variable. The data is ! stored so that the largest pressure levels come first (ie, ! the surface). The flag value for the surface pressure is ! 200100 Pa. current_slab => current_var%slab ! Fill in the surface value into the first level of ! this variable. Save the accompanying pressure level ! that defines the surface. All of the data in the linked list ! are stored in 2d slices. We build the 3d array with one slice ! at a time, then move to the next variable name. fg_vertical_index = 1 all_3d(:,:,fg_vertical_index,fg_variables_up_index) = current_slab%field ! Because this may be a level other than the surface, we force ! the pressure level to be defined as the surface. The first ! level in the array has to be the surface if this is a 3d array. output_levels(fg_vertical_index) = 200100 IF ( verbose ) THEN PRINT '(A,I8,A,A,A,I8)','Store data in arrays: Found 3d data index #',fg_variables_up_index, & ' named ',all_names_up(fg_variables_up_index)(1:8),' pressure = ', & output_levels(fg_vertical_index) END IF ! If this is the actual 200100 Pa level, this was indeed ! the surface data that we just loaded. In that case we ! can increment the level. Otherwise, we keep the level ! the same and reuse the data at the next level. This could ! happen, for example, if there was no surface RH field. IF ( NINT(current_slab%level) .EQ. 200100 ) THEN current_slab => current_slab%next_slab END IF ! Initializations required for vertical interpolation. slab_bottom = current_slab%field pressure_bottom = current_slab%level IF ( num_new_levels .GT. 0 ) THEN vertically_interpolate = .TRUE. END IF ! Now that we have the surface out of the way, loop ! through all of the other levels. We have to check for the ! possibility of needing to do a vertical interpolation. For ! each level (interpolated or not), we will increment the ! vertical index, store the 3d datas 2d slice, and save the ! pressure level for this slice. get_all_levels : DO WHILE ( ASSOCIATED(current_slab) ) ! Increment vertical location of the 3d array level. This ! is effectively the vertical index in the 3d array. Since we ! are in this loop, we know that there is another data level ! to process. fg_vertical_index = fg_vertical_index + 1 ! Check for vertical interpolation. This is a test that surrounds a ! loop through all of the vertical levels to be interpolated to. ! The only "oops" is if a user tries to interpolate between the ! surface (200100 Pa) and 100000 Pa (for example). This would be seen ! as a valid computation. Hopefully, no one tries. do_vertical_interpolation : IF ( vertically_interpolate ) THEN ! To vertically interpolate, we need the surrounding pressure ! levels. We saved the previous pressure (pressure_bottom) at ! the end of this loop. Here is the current pressure (pressure_top). pressure_top = current_slab%level ! Loop over all possible new levels to interpolate to. check_for_new_levels : DO new_levels_index = 1 , num_new_levels ! The candidate new pressure level. pressure_middle = new_levels(new_levels_index) ! If the new vertical level is between the previous level and the ! current level, then perform the vertical interpolation linear in ! pressure. IF ( ( pressure_bottom .GT. pressure_middle ) .AND. & ( pressure_middle .GT. pressure_top ) ) THEN slab_top = current_slab%field slab_middle = ( (pressure_bottom-pressure_middle)*slab_top + & (pressure_middle-pressure_top )*slab_bottom ) / & (pressure_bottom-pressure_top ) ! A successful venture through a weighted mean. Now, save the ! level for this pressure, the actual interpolated data, and finally ! increment the vertical index to reflect the new level. output_levels(fg_vertical_index) = NINT(pressure_middle) all_3d(:,:,fg_vertical_index,fg_variables_up_index) = slab_middle IF ( verbose ) THEN PRINT '(A,I8,A,A,A,I8)','Store data in arrays: Found 3d data index #', & fg_variables_up_index, & ' named ',all_names_up(fg_variables_up_index)(1:8), & ' pressure = ', output_levels(fg_vertical_index) END IF fg_vertical_index = fg_vertical_index + 1 ! Insure that the data that was requested as new levels was not available ! in the original data set. For example, if we ask for 60000 Pa as a new level, but ! the first guess data already provides it, there will be problems in the count of ! the vertical levels. We have to check each variable since we do not have to have ! the same number of levels for any given field. ELSE IF ( ( ABS ( pressure_bottom - pressure_middle ) .LT. 1 ) .OR. & ( ABS ( pressure_middle - pressure_top ) .LT. 1 ) ) THEN IF ( verbose ) THEN PRINT '(A)','Found a duplicated vertical level in the requested new ' PRINT '(A)','levels and in the first guess field.' END IF CYCLE check_for_new_levels END IF END DO check_for_new_levels END IF do_vertical_interpolation ! Store the 3d data and vertical location for this level. Whether or ! not there was a vertical interpolation, we still have the current level ! in the linked list to extract. Save the pressure level and this slice ! of the 3d data. output_levels(fg_vertical_index) = NINT(current_slab%level) all_3d(:,:,fg_vertical_index,fg_variables_up_index) = current_slab%field IF ( verbose ) THEN PRINT '(A,I8,A,A,A,I8)','Store data in arrays: Found 3d data index #', & fg_variables_up_index, & ' named ',all_names_up(fg_variables_up_index)(1:8), & ' pressure = ', output_levels(fg_vertical_index) END IF ! If we are doing vertical interpolations, well, then it is a requirement ! to save things for the next go round. This is the "lower" level that ! will be used to check the trap of the middle level. IF ( vertically_interpolate ) THEN pressure_bottom = current_slab%level slab_bottom = current_slab%field END IF ! Point to the next slab, which in this case, is point to the ! next vertical level. current_slab => current_slab%next_slab END DO get_all_levels END IF twod_or_threed ! We have done all of the levels for this variable. Point to ! the next variable. current_var => current_var%next_variable END DO search_all_variables END SUBROUTINE proc_list_to_array !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE proc_list_to_array_noniso ! ! This routine takes the linked list of 2D slabs combines ! them appropriately into 2D and 3D arrays to populate ! the all_3d and all_2d structures defined in the ! module_gridded_data for non-isobaric data. ! ! HISTORY ! ======= ! Apr 2002 - Brent L. Shaw, NOAA/FSL ! New routine to support non-isobaric levels IMPLICIT NONE INTEGER :: fg_vertical_index LOGICAL :: any_surface_data LOGICAL :: have_press3d ! Initialize the two name arrays (both for surface and upper-air). all_names_sfc = ' ' all_names_up = ' ' fg_variables_sfc_index = 0 fg_variables_up_index = 0 have_press3d = .false. ! Find out how many levels are in the 3d data. It is possible ! that some fields do not have as many levels as other fields ! (such as RH, there may not be a surface value). Loop through ! all of the variable names and find the maximum of the number ! of levels. any_surface_data = .FALSE. current_var => head number_of_original_levels = 0 find_max_levels : DO WHILE ( ASSOCIATED(current_var) ) number_of_original_levels = MAX ( number_of_original_levels , current_var%num_levels ) IF ( verbose) THEN PRINT '(A,I8,A,A)','Found ',current_var%num_levels,' levels of data for variable ',current_var%name END IF IF ( ( current_var%num_levels .GT. 1 ) .AND. ( .NOT. any_surface_data ) ) THEN IF ( ABS(current_var%slab%level-200100.) .LT. 1 ) THEN any_surface_data = .TRUE. PRINT '(A)', 'We are processing non-isorbaric data. For now, this' PRINT '(A)', 'program assumes surface data is the lowest level in the 3d' PRINT '(A)', 'data, and that surface fields will not be provided in' PRINT '(A)', 'a level labeled 200100. Pa. For now, you need to' PRINT '(A)', 'remove such entries from you Vtable for this data source and' PRINT '(A)', 'rerun grib_prep.' STOP 'sfc_level_found_in_nonisobaric_data' END IF END IF IF ((current_var%num_levels .GT. 1) .AND. & (current_var%name(1:8).EQ.'PRESSURE')) THEN have_press3d = .true. ENDIF current_var => current_var%next_variable END DO find_max_levels IF (.NOT. have_press3d) THEN PRINT '(A)', 'No 3D pressure array found for this non-isobaric data set.' PRINT '(A)', 'Add a PRESSURE entry to your Vtable and re-run gribprep.' STOP 'No_3d_pressure' ENDIF IF ( verbose ) THEN PRINT '(A,I8,A)','Found ',number_of_original_levels,' as maximum number of original vertical levels.' END IF ! Correct the pressure level data for a few things. 1) If we have a ! level and the same level is requested as a new level, decrement the ! number of original levels current_var => head ! The vertical array of output_levels is ALLOCATEd. This is the sum ! of the data from the first guess file. ALLOCATE (output_levels(number_of_original_levels) ) ! ALLOCATE the space required for the 3d data. One array for ! each variable. The vertical dimension is the number of the ! input levels (from the first guess fields themselves). ALLOCATE ( all_3d(nx,ny,number_of_original_levels,num3d_fg) ) ! ALLOCATE the space for the 2d arrays. ! ALLOCATE ( all_2d(nx,ny,num2d_fg*3) ) ! COMMENTED OUT BLS-20030306 ALLOCATE ( all_2d(nx,ny,num2d_fg) ) ! REPLACES ABOVE BLS-20030306 ! Initialize some of the required counters to zero. fg_variables_sfc_index = 0 fg_variables_up_index = 0 fg_vertical_index = 0 ! We are going to loop over all of the variables, so start ! by pointing to the first one. We loop until we are no ! longer pointing to an ASSOCIATEd variable. current_var => head search_all_variables : DO WHILE ( ASSOCIATED(current_var) ) ! If there is only a single level, this is 2d data. As ! is seen below, the 2d data is a piece of cake to store. twod_or_threed : IF ( current_var%num_levels .EQ. 1 ) THEN ! Find where to stick this data. Fill an array of the ! name, then store the 2d data in the same location of ! the ALLOCATEd 2d array. Increment the index of found ! data, then load the data, then save the name of the data. ! This is easy. fg_variables_sfc_index = fg_variables_sfc_index + 1 all_2d(:,:,fg_variables_sfc_index) = current_var%slab%field all_names_sfc(fg_variables_sfc_index) = current_var%name IF ( verbose ) THEN PRINT '(A,I8,A,A)','Store data in arrays: Found 2d data index #',fg_variables_sfc_index, & ' named ',all_names_sfc(fg_variables_sfc_index)(1:8) END IF ! If there is more than a single level, this is a 3d array. As ! is judged by the following rest of the IF test, this is a bit ! more involved. There are vertical levels to find, possible ! interpolations and maybe missing surface values. Fasten your ! seat belt. ELSE IF ( current_var%num_levels .GT. 1 ) THEN ! Increment the index of found upper-air data. This is a ! count of the total number of 3d fields in the data set. fg_variables_up_index = fg_variables_up_index + 1 ! Save the name of this 3d array. all_names_up(fg_variables_up_index) = current_var%name ! Point to the first level of this variable. The data is ! stored so that the lowest level values come first (ie, ! the surface). The flag value for the surface pressure is ! 200100 Pa. current_slab => current_var%slab ! Fill in the surface value into the first level of ! this variable. Save the accompanying pressure level ! that defines the surface. All of the data in the linked list ! are stored in 2d slices. We build the 3d array with one slice ! at a time, then move to the next variable name. fg_vertical_index = 1 all_3d(:,:,fg_vertical_index,fg_variables_up_index) = current_slab%field output_levels(fg_vertical_index) = current_slab%level IF ( verbose ) THEN PRINT '(A,I8,A,A,A,I8)','Store data in arrays: Found 3d data index #',fg_variables_up_index, & ' named ',all_names_up(fg_variables_up_index)(1:8),' level # = ', & output_levels(fg_vertical_index) END IF ! Now that we have the surface out of the way, loop ! through all of the other levels. We have to check for the ! possibility of needing to do a vertical interpolation. For ! each level (interpolated or not), we will increment the ! vertical index, store the 3d datas 2d slice, and save the ! pressure level for this slice. ! Unlike the isobaric version, we will go ahead and increment ! the slab right now, so we dont end up with level 1 in both ! level 1 and 2. Bug fix by BLS on 20030815 current_slab => current_slab%next_slab get_all_levels : DO WHILE ( ASSOCIATED(current_slab) ) ! Increment vertical location of the 3d array level. This ! is effectively the vertical index in the 3d array. Since we ! are in this loop, we know that there is another data level ! to process. fg_vertical_index = fg_vertical_index + 1 ! Store the 3d data and vertical location for this level. Whether or ! not there was a vertical interpolation, we still have the current level ! in the linked list to extract. Save the pressure level and this slice ! of the 3d data. output_levels(fg_vertical_index) = NINT(current_slab%level) all_3d(:,:,fg_vertical_index,fg_variables_up_index) = current_slab%field IF ( verbose ) THEN PRINT '(A,I8,A,A,A,I8)','Store data in arrays: Found 3d data index #', & fg_variables_up_index, & ' named ',all_names_up(fg_variables_up_index)(1:8), & ' level # = ', output_levels(fg_vertical_index) END IF ! Point to the next slab, which in this case, is point to the ! next vertical level. current_slab => current_slab%next_slab END DO get_all_levels END IF twod_or_threed ! We have done all of the levels for this variable. Point to ! the next variable. current_var => current_var%next_variable END DO search_all_variables END SUBROUTINE proc_list_to_array_noniso !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END MODULE hinterp_gribprep