SUBROUTINE proc_ingest_first_guess USE constants USE diags USE first_guess_data USE file_data USE gauss USE gridded_data USE header_data USE horiz_interp USE link_list_info USE map_utils USE namelist_info USE terrain_data IMPLICIT NONE TYPE(proj_info) :: proj REAL , DIMENSION(:,:) , ALLOCATABLE :: array , array_backwards , & dum2d , longitude , latitude , tempgrd , landuse , terrain REAL :: xlo , yla , dlon , dlat , known_i , known_j INTEGER :: ilo , jla , i , j , icd , icd_old LOGICAL :: initial_fg_read = .TRUE. LOGICAL :: found_any_files = .FALSE. INTEGER :: close_count , i_close , j_close REAL :: close_data REAL :: point_x ! x location of point with unknown latitude, usually center REAL , DIMENSION (:,:), ALLOCATABLE :: xloc , yloc REAL :: angle LOGICAL :: converged , similar , wrap REAL :: xlonc_mod CHARACTER (LEN=24) :: hdate_input CHARACTER (LEN= 9) :: field_input CHARACTER (LEN=25) :: units_input CHARACTER (LEN=46) :: description_input LOGICAL :: constant_file ! additions D.Ovens (ovens@uw.edu) to handle WRF ungrib output character (len=8) :: startloc character (len=32) :: map_source logical :: is_wind_grid_rel real :: earth_radius ! end additions ! 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. all_names_sfc = ' ' all_names_up = ' ' num3d_fg = 0 num2d_fg = 0 num1d_fg = 0 ! Initialize these known_i and known_j values. If this is a projected ! data set, these will be set. If this is just your regular run of the ! mill lat/lon spacing, this is all we get. known_i = 1 known_j = 1 ! We know how big the model grid should be, so we now ALLOCATE it ! as the correct size array. Note that this is done once per program ! run. initial_fg_read = .TRUE. IF ( initial_fg_read ) THEN ALLOCATE(dum2d(jxf,ixf)) ALLOCATE(landuse(jxf,ixf)) ALLOCATE(tempgrd(jxf,ixf)) ALLOCATE(terrain(jxf,ixf)) ALLOCATE(latitude(jxf,ixf)) ALLOCATE(longitude(jxf,ixf)) ALLOCATE(xloc(jxf,ixf)) ALLOCATE(yloc(jxf,ixf)) 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. diff_fg_files : DO diff_files_per_single_time = 1 , rec3%max_files_per_single_time + & rec3%max_files_constants ! Get a UNIT number to OPEN. first_guess_file = -1 first_guess_unit : DO file_index = file_index_min , file_index_max INQUIRE (UNIT=file_index, OPENED=already_opened) IF ( .NOT. already_opened) THEN first_guess_file = file_index EXIT first_guess_unit END IF END DO first_guess_unit IF ( first_guess_file .EQ. -1 ) THEN PRINT '(A)','Could not find a spare unit for the FIRST GUESS data.' print '(A)','first_guess_file_unit' call exit(1) END IF ! F90 allows "as is" positioning, so make sure that does not ! slip in some how. REWIND ( UNIT=first_guess_file ) ! 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. rec3%max_files_constants ) THEN file_name_len = LEN_TRIM(rec3%constants_full_name(diff_files_per_single_time)) file_name(1:file_name_len) = TRIM(rec3%constants_full_name(diff_files_per_single_time)) IF (file_name(1:1) .EQ. ' ' ) THEN PRINT '(A)','Your constant file names may have "spaces" at the beginning. Remove them.' PRINT '(A)','zap_the_spaces_in_const_filename' call exit(1) END IF constant_file = .TRUE. IF ( rec4%print_file ) THEN PRINT '(A,A)','Trying to OPEN constant file: ',TRIM(file_name) END IF ELSE file_name_len = LEN_TRIM(rec3%root(diff_files_per_single_time-rec3%max_files_constants)) file_name(1:file_name_len) = TRIM(rec3%root(diff_files_per_single_time-rec3%max_files_constants)) IF (file_name(1:1) .EQ. ' ' ) THEN PRINT '(A)','Your file names may have "spaces" at the beginning. Remove them.' print '(A)', 'zap_the_spaces_in_filename' call exit(1) END IF IF ( (rec1%time_increment_second/60)*60 .NE. rec1%time_increment_second ) THEN file_name(file_name_len+1:) = ':' // rec1%current_date(1:19) ELSE IF ( (rec1%time_increment_second/3600)*3600 .NE. rec1%time_increment_second ) THEN file_name(file_name_len+1:) = ':' // rec1%current_date(1:16) ELSE file_name(file_name_len+1:) = ':' // rec1%current_date(1:13) END IF constant_file = .FALSE. IF ( rec4%print_file ) THEN PRINT '(A,A)','Trying to OPEN first-guess file: ',TRIM(file_name) END IF END IF OPEN ( UNIT = first_guess_file , & FILE = TRIM(file_name) , & STATUS = 'OLD' , & FORM = 'UNFORMATTED' , & IOSTAT = open_status ) IF ( open_status .NE. 0 ) THEN PRINT '(A,A)' , 'Error OPENing first_guess file: ',TRIM(file_name) CYCLE diff_fg_files ELSE found_any_files = .TRUE. PRINT '(A,A)' , 'Successfully OPENed first_guess file: ',TRIM(file_name) END IF ! These are some initializations to help speed up the interpolation ! if the domain requires the similar latitude longitude computation ! for each level, for each variable. idim_old = -1 jdim_old = -1 llflag_old = -1 lat1_old = -1 lon1_old = -1 deltax_old = -1 deltay_old = -1 xlonc_old = -1 truelat1_old = -1 truelat2_old = -1 icd_old = -1 icd = -2 ! This is the big loop over all of the data in this particular ! file. There is an outer loop over different files. this_fg_file_only : DO ! Read header data from FIRST GUESS file. READ ( UNIT=first_guess_file , 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. IF ( version .EQ. 1 ) THEN PRINT '(A)','This is an old format (IVERSION=1), the character strings may be incorrect.' READ ( UNIT=first_guess_file ) & hdate_input, field_input, units_input, description_input, level, idim, jdim, llflag xfcst = 0. hdate = hdate_input(1:19) field = field_input(1:8 ) units = field_input(9:9 ) // units_input(1:15) // ' ' description = units_input(16:25) // description_input(1:36) ELSE IF ( version .EQ. 2 ) THEN READ ( UNIT=first_guess_file ) & hdate_input, xfcst, field_input, units_input, description_input, level, idim, jdim, llflag xfcst = xfcst * 60. hdate = hdate_input(1:19) field = field_input(1:8 ) units = field_input(9:9 ) // units_input(1:15) // ' ' description = units_input(16:25) // description_input(1:36) ELSE IF ( version .EQ. 3 ) THEN READ ( UNIT=first_guess_file ) & hdate_input, xfcst, field_input, units_input, description_input, level, idim, jdim, llflag xfcst = xfcst * 60. hdate = hdate_input(1:19) field = field_input(1:8 ) units = units_input(1:24) description = description_input(1:46) ELSE IF ( version .EQ. 5 ) THEN READ ( UNIT=first_guess_file ) & hdate_input, xfcst, map_source, field_input, units_input, description_input, level, idim, jdim, llflag xfcst = xfcst * 60. hdate = hdate_input(1:19) field = field_input(1:8 ) units = units_input(1:24) description = description_input(1:46) ELSE PRINT '(A)','This program fully supports IVERSION=3; less so on IVERSION=2; and IVERSION=1 is old.' PRINT '(A,I4,A)','You tried version ',version,'. No can do.' PRINT '(A)','wrong_version_of_input' call exit(1) END IF ! 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. ! Current choices are: ! 0 - regular latitude/longitude projection ! 1 - Mercator projection ! 3 - Lambert conformal projection ! 4 - Gaussian projection ! 5 - polar stereographic projection IF ( ( llflag .NE. 0 ) .AND. & ( llflag .NE. 1 ) .AND. & ( llflag .NE. 3 ) .AND. & ( llflag .NE. 4 ) .AND. & ( llflag .NE. 5 ) ) THEN PRINT '(A)','This program can only handle data with specific projections.' PRINT '(A,I8,A)','The input data uses projection ',llflag,', which is not acceptable.' PRINT '(A)','Projection_type_unknown' call exit(1) END IF IF ( rec4%print_echo ) THEN PRINT '(A,1X,F10.3,1X,A,1X,I8,2I4,I2)' , hdate, xfcst, field, NINT(level), idim, jdim, 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. rec1%current_date(1:13) ) .AND. ( .NOT. constant_file ) ) THEN PRINT '(A,A,A)','Expected to find date = ',rec1%current_date(1:13),' in the input file.' PRINT '(A,A,A)','Instead, the date ',hdate(1:13),' was in the header.' PRINT '(A)','Wrong_date_in_first_guess_header' call exit(1) ELSE IF ( constant_file ) THEN hdate(1:19) = rec1%current_date(1:19) IF ( rec4%print_echo ) THEN PRINT '(A,A,A)','Constant file assumed valid for ',rec1%current_date(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 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 didn't find the name, something went wrong. IF ( .NOT. found_name ) THEN PRINT '(A)','Error in finding names for 3d/2d search for surface data.' PRINT '(A)','search_num2d' call exit(1) END IF END IF ! Same search in the upper-air data. IF ( NINT(level) .EQ. 100000 ) 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 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 didn't find the name, something went wrong. IF ( .NOT. found_name ) THEN PRINT '(A)','Error in finding names for 3d/2d search for upper-air data.' PRINT '(A)','search_num3d' call exit(1) 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. If this is a domain that ! is not constant latitude/longitude, then we need to eventually compute an ! additional important value - namely the center latitude. IF ( llflag .EQ. 0 ) THEN if (version.eq.5) then READ ( UNIT=first_guess_file ) startloc, startlat, startlon, deltalat, deltalon, earth_radius else READ ( UNIT=first_guess_file ) startlat, startlon, 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 ( ( startlon .GT. 180 ) .AND. ( startlon .LE. 360 ) ) THEN startlon = startlon - 360 END IF IF ( rec4%print_echo ) THEN PRINT '(4F10.4)', startlat, startlon, deltalat, deltalon END IF IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN wind_llflag = llflag END IF known_i = 1 known_j = 1 ELSE IF ( llflag .EQ. 1 ) THEN if (version.eq.5) then READ ( UNIT=first_guess_file ) startloc, lat1, lon1, deltax, deltay, truelat1, earth_radius else READ ( UNIT=first_guess_file ) lat1, lon1, deltax, deltay, truelat1 endif IF ( rec4%print_echo ) THEN PRINT '(6F10.4)', lat1, lon1, deltax, deltay, truelat1 END IF cone_fac = 0.0 truelat2 = -99. angle = 360. * deltax * (idim-1) / & ( 2. * pi * earth_radius_m/1000. * COS(truelat1*radians_per_degree) ) lon2 = lon1 + angle IF ( lon2 .GT. 180 ) THEN lon2 = lon2 - 360 END IF IF ( lon1 .LT. lon2 ) THEN xlonc = ( lon2 + lon1 ) / 2. ELSE xlonc = ( ( 360 + lon2 ) + lon1 ) / 2. END IF IF ( xlonc .GT. 180 ) THEN xlonc = xlonc - 360 END IF ! Hey, was this the same as last time? Can we avoid an expensive computation? ! What we want to find out is if this is the same projection information, we ! do not care about the specifics of the field or level. IF ( ( idim_old .EQ. idim ) .AND. ( jdim_old .EQ. jdim ) .AND. ( llflag_old .EQ. llflag ) .AND. & ( lat1_old .EQ. lat1 ) .AND. ( lon1_old .EQ. lon1 ) .AND. & ( deltax_old .EQ. deltax ) .AND. ( deltay_old .EQ. deltay ) .AND. ( xlonc_old .EQ. xlonc ) .AND. & ( truelat1_old .EQ. truelat1 ) ) THEN similar = .TRUE. ELSE similar = .FALSE. END IF ! There are a few values that we need to compute from the available data. IF ( .NOT. similar ) THEN IF ( lon1 .LT. -180. ) lon1 = lon1 + 360. IF ( lon1 .GT. 180. ) lon1 = lon1 - 360. IF ( xlonc .LT. -180. ) xlonc = xlonc + 360. IF ( xlonc .GT. 180. ) xlonc = xlonc - 360. IF ( ( deltax .GT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = 1. known_j = 1. ELSE IF ( ( deltax .GT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = 1. known_j = jdim ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = idim known_j = 1. ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = idim known_j = jdim END IF if(nint(known_i*known_j).gt.1) then stop 'fix_this_mercator' endif CALL map_set(proj_merc,lat1,lon1,1.,1.,deltax*1000.,xlonc,truelat1,truelat2,proj) END IF IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN wind_llflag = llflag wind_xlonc = xlonc wind_xlatc = xlatc wind_cone_fac = cone_fac END IF ELSE IF ( llflag .EQ. 3 ) THEN if (version.eq.5) then READ ( UNIT=first_guess_file ) startloc, lat1, lon1, deltax, deltay, xlonc, truelat1, truelat2, earth_radius else READ ( UNIT=first_guess_file ) lat1, lon1, deltax, deltay, xlonc, truelat1, truelat2 endif IF ( rec4%print_echo ) THEN PRINT '(7F10.4)', lat1, lon1, deltax, deltay, xlonc, truelat1, truelat2 END IF ! Hey, was this the same as last time? Can we avoid an expensive computation? ! What we want to find out is if this is the same projection information, we ! do not care about the specifics of the field or level. IF ( ( idim_old .EQ. idim ) .AND. ( jdim_old .EQ. jdim ) .AND. ( llflag_old .EQ. llflag ) .AND. & ( lat1_old .EQ. lat1 ) .AND. ( lon1_old .EQ. lon1 ) .AND. & ( deltax_old .EQ. deltax ) .AND. ( deltay_old .EQ. deltay ) .AND. ( xlonc_old .EQ. xlonc ) .AND. & ( truelat1_old .EQ. truelat1 ) .AND. ( truelat2_old .EQ. truelat2 ) ) THEN similar = .TRUE. ELSE similar = .FALSE. END IF ! There are a few values that we need to compute from the available data. IF ( .NOT. similar ) THEN IF ( lon1 .LT. -180. ) lon1 = lon1 + 360. IF ( lon1 .GT. 180. ) lon1 = lon1 - 360. IF ( xlonc .LT. -180. ) xlonc = xlonc + 360. IF ( xlonc .GT. 180. ) xlonc = xlonc - 360. IF ( ( deltax .GT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = 1. known_j = 1. ELSE IF ( ( deltax .GT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = 1. known_j = jdim ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = idim known_j = 1. ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = idim known_j = jdim END IF if(nint(known_i*known_j).gt.1) then stop 'fix_this_lambert' endif CALL map_set(proj_lc,lat1,lon1,1.,1.,deltax*1000.,xlonc,truelat1,truelat2,proj) END IF IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN wind_llflag = llflag wind_xlonc = xlonc wind_xlatc = xlatc wind_cone_fac = proj%cone END IF ELSE IF ( llflag .EQ. 4 ) THEN if (version.eq.5) then READ ( UNIT=first_guess_file ) startloc, startlat, startlon, deltalon , nlat, earth_radius else READ ( UNIT=first_guess_file ) startlat, startlon, deltalon , nlat 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 ( ( startlon .GT. 180 ) .AND. ( startlon .LE. 360 ) ) THEN startlon = startlon - 360 END IF IF ( rec4%print_echo ) THEN PRINT '(4F10.4)', startlat, startlon, deltalon , nlat END IF IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN wind_llflag = llflag END IF known_i = 1 known_j = 1 ELSE IF ( llflag .EQ. 5 ) THEN if (version.eq.5) then READ ( UNIT=first_guess_file ) startloc, lat1, lon1, deltax, deltay, xlonc, truelat1, earth_radius READ ( UNIT=first_guess_file ) is_wind_grid_rel else READ ( UNIT=first_guess_file ) lat1, lon1, deltax, deltay, xlonc, truelat1 endif ! If this is a *special* AFWA domain of seaice, fix the center longitude to ! what is considered an MM5 standard. The specifics: it is a se ice field, ! it is polar stereographic, it has a negatice true lat, the (1,1) location ! is in the upper left hand corner. Based on these, swap the central longitude. ! This replaces the mod sent to AFWA 23 April 2003. IF ( ( field(1:6) .EQ. 'SEAICE' ) .AND. & ( truelat1 .LT. 0. ) .AND. & ( deltax .GT. 0. ) .AND. & ( deltay .LT. 0. ) ) THEN PRINT *,'AFWA OLD LONC=',xlonc xlonc = xlonc - 180. END IF IF ( rec4%print_echo ) THEN PRINT '(6F10.4)', lat1, lon1, deltax, deltay, xlonc, truelat1 END IF cone_fac = 1.0 truelat2 = -99. ! Hey, was this the same as last time? Can we avoid an expensive computation? ! What we want to find out is if this is the same projection information, we ! do not care about the specifics of the field or level. IF ( ( idim_old .EQ. idim ) .AND. ( jdim_old .EQ. jdim ) .AND. ( llflag_old .EQ. llflag ) .AND. & ( lat1_old .EQ. lat1 ) .AND. ( lon1_old .EQ. lon1 ) .AND. & ( deltax_old .EQ. deltax ) .AND. ( deltay_old .EQ. deltay ) .AND. ( xlonc_old .EQ. xlonc ) .AND. & ( truelat1_old .EQ. truelat1 ) ) THEN similar = .TRUE. ELSE similar = .FALSE. END IF ! There are a few values that we need to compute from the available data. IF ( .NOT. similar ) THEN IF ( lon1 .LT. -180. ) lon1 = lon1 + 360. IF ( lon1 .GT. 180. ) lon1 = lon1 - 360. IF ( xlonc .LT. -180. ) xlonc = xlonc + 360. IF ( xlonc .GT. 180. ) xlonc = xlonc - 360. IF ( ( deltax .GT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = 1. known_j = 1. ELSE IF ( ( deltax .GT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = 1. known_j = jdim ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .GT. 0 ) ) THEN known_i = idim known_j = 1. ELSE IF ( ( deltax .LT. 0 ) .AND. ( deltay .LT. 0 ) ) THEN known_i = idim known_j = jdim END IF CALL map_set(proj_ps,lat1,lon1,known_i,known_j,ABS(deltax*1000.),xlonc,truelat1,truelat2,proj) END IF IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN wind_llflag = llflag wind_xlonc = xlonc wind_xlatc = xlatc wind_cone_fac = cone_fac END IF END IF ! 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. IF ( ( idim .NE. idim_old ) .OR. ( jdim .EQ. jdim_old ) ) THEN IF ( ALLOCATED ( array ) ) THEN DEALLOCATE ( array ) END IF ALLOCATE (array(idim,jdim)) END IF ! Read the 2d array. READ ( UNIT=first_guess_file ) array ! There are data sets that set HUGE POSITIVE flag values for masked fields. Our ! humble convention, of course, is to instead have HUGE NEGATIVE flags. WHERE ( array .GT. 1.E19 ) array = -1.E20 END WHERE ! If the projected data is flipped upside down or side to side, then we'll ! just handle it here expeditiously. I cannot fathom changing indices ! below in each of the interpolation blocks just to handle this situation. ! First, this is only if the data is projected (polar only right now), and ! second it is only if the known_i or known_j values are not = 1. IF ( ( llflag .EQ. 5 ) .AND. & ( ( NINT(known_i) .NE. 1 ) .OR. ( NINT(known_j) .NE. 1 ) ) ) THEN IF ( ALLOCATED ( array_backwards ) ) THEN DEALLOCATE ( array_backwards ) END IF ALLOCATE (array_backwards(idim,jdim)) ! There are four possible origin locations that may be chosen for ! domain, let's call them A, B, C, and D. If the origin is set as ! point A, then we do nothing (we are not even in this IF test!). ! We handle each of the three remaining flips. ! C D ! ----------------------------- ! | | ! | | ! | | ! | | ! | | ! | | ! | | ! | | ! | | ! | | ! | | ! ----------------------------- ! A B ! Flip B : known_i = dim, known_j = 1 IF ( ( NINT(known_i) .EQ. idim ) .AND. ( NINT(known_j) .EQ. 1 ) ) THEN DO j = 1 , jdim DO i = 1 , idim array_backwards(idim-i+1,j) = array(i,j) END DO END DO ! Flip C : known_i = 1, known_j = jdim ELSE IF ( ( NINT(known_i) .EQ. 1 ) .AND. ( NINT(known_j) .EQ. jdim ) ) THEN DO j = 1 , jdim DO i = 1 , idim array_backwards(i,jdim-j+1) = array(i,j) END DO END DO ! Flip D : known_i = dim, known_j = jdim ELSE IF ( ( NINT(known_i) .EQ. idim ) .AND. ( NINT(known_j) .EQ. jdim ) ) THEN DO j = 1 , jdim DO i = 1 , idim array_backwards(idim-i+1,jdim-j+1) = array(i,j) END DO END DO END IF ! Put the flipped data back in the array with which we wish to play. array = array_backwards END IF ! 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 ( rec4%print_echo ) THEN 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. rec2%ptop ) THEN IF ( rec4%print_interp) THEN PRINT '(A,F7.0,A,F7.0,A,A,A)','The NAMELIST ptop variable (',rec2%ptop, & ' Pa) requires that this level (',level,' Pa) of ', & field,' be omitted.' END IF ! 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 a surface pressure field, we do not want it. The sea ! level pressure is important, the surface pressure is a function ! of the terrain elevation of the first guess data. Besides being ! fairly useless, if we end up with two PRESSURE levels (sea level ! and surface), the darn program thinks that PRESSURE is a 3d field. IF ( ( NINT(level) .EQ. 200100 ) .AND. ( field .EQ. 'PRESSURE' ) ) THEN PRINT '(A)','Not processing surface pressure as a 2d array.' ! Get to the bottom of the loop and get the next data from the file. CYCLE this_fg_file_only END IF ! Version 2 of the masked tests: ! Some of the first guess data sets have masked data. If this ! is a masked field, interpolating has problems if one of the ! 4 surrounding points is a flag value (typically zero - real helpful). ! So, check to see if this is a known masked field, and if so, ! check to see if there are some flag values, and if so, then do some ! simple things to see how we get the data. The easiest case is when ! all 4 surrounding points are real data (tests #1 and #2), we just use ! the real data (surprise!). When there are no valid data values in the ! surrounding 4 points (tests #3 and #4), then we look to expand the ! search to 16 points. Failing that, we assign a fixed value (15 C for ! soil temperature, 0 for both soil moisture and sea ice). The last tests ! are when there are some (not 4 and not 0) values to be used for the ! masked interpolation (tests #5 and #6). In this case we use as much ! data as is possible. The last two tests (#7 and #8) are the null data ! cases for when we have a land grid point with a water based field (or ! a water point with a land based field). These values are set to 0. ! Version 3 of the tests: ! If there is a location with sea ice, then we have to allow there to be ! a soil temperature. ! Version 4 of the tests: ! Allow the TEMPGRD, if it exists, to help with islands. If TEMPGRD is ! not available, use the function of the latitude for the soil temps. IF ( ( field(1:5) .EQ. 'SOILT' ) .OR. & ( field(1:7) .EQ. 'SOILHGT' ) .OR. & ( field(1:5) .EQ. 'SOILM' ) .OR. & ( field(1:5) .EQ. 'SOILW' ) .OR. & ( field(1:6) .EQ. 'MONALB' ) .OR. & ( field(1:8) .EQ. 'ALBSNOMX' ) .OR. & ( field(1:3) .EQ. 'SEA' ) .OR. & ( field(1:7) .EQ. 'CANOPYM' ) .OR. & ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. & ( field(1:5) .EQ. 'WEASD' ) ) THEN ! Store the latitude, longitude, elevation, land use and the deep soil ! temperature from the TERRAIN file into some more attractively named local ! variables. latitude = ter_data(latitude_cross_index)%array longitude = ter_data(longitude_cross_index)%array terrain = ter_data(terrain_index)%array IF ( tempgrd_index .GT. 0 ) THEN tempgrd = ter_data(tempgrd_index)%array END IF landuse = ter_data(land_use_index)%array icd_old = icd icd = 1 ! For the regular lat/lon grids, the longitude field needs to be ! massaged to allow correct horizontal interpolation across ! Greenwich and the date line. IF ( llflag .EQ. 0 )THEN IF ( deltalon .LT. 0 ) THEN WHERE ( longitude .GT. startlon ) longitude = longitude - 360. END WHERE ELSE IF ( deltalon .GT. 0 ) THEN WHERE ( longitude .LT. startlon ) longitude = longitude + 360. END WHERE END IF END IF ! Loop through each of the model grid points. ! pjj/cray !$OMP parallel do default(shared) & !$OMP private(i,j,yla,xlo,ilo,jla,dlon,dlat,xlonc_mod,deltalat,deltalon,close_count,close_data) mask_i_loop : DO i = 1, ixf mask_j_loop : DO j = 1, jxf yla = latitude(j,i) xlo = longitude(j,i) ! Make the interpolation specific to the first guess input projection. IF ( llflag .EQ. 0 ) THEN ilo = MAX ( MIN(FLOOR((xlo-startlon)/deltalon) + 1, idim-1) , 1 ) jla = MAX ( MIN(FLOOR((yla-startlat)/deltalat) + 1, jdim-1) , 1 ) dlon = xlo - ( startlon + (ilo-1) * deltalon ) IF ( ABS(dlon) .GT. ABS(deltalon) ) THEN PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) dlon = 1. END IF dlat = yla - ( startlat + (jla-1) * deltalat ) ELSE IF ( ( llflag .EQ. 1 ) .OR. ( llflag .EQ. 3 ) .OR. ( llflag .EQ. 5 ) ) THEN IF (xlonc .GT. 180 ) THEN xlonc_mod = xlonc - 360. ELSE IF (xlonc .LT. -180 ) THEN xlonc_mod = xlonc + 360. ELSE xlonc_mod = xlonc END IF IF ( ( .NOT. similar ) .OR. ( icd .NE. icd_old ) ) THEN IF ( ( llflag .EQ. 1 ) .OR. & ( llflag .EQ. 3 ) .OR. & ( llflag .EQ. 5 ) ) THEN CALL latlon_to_ij(proj,yla,xlo,xloc(j,i) , yloc(j,i)) END IF IF ( ( NINT(xloc(j,i)) .GT. idim ) .OR. & ( NINT(yloc(j,i)) .GT. jdim ) .OR. & ( NINT(xloc(j,i)) .LT. 1 ) .OR. & ( NINT(yloc(j,i)) .LT. 1 ) ) THEN PRINT '(A)','The chosen domain from the terrain data is not completely contained & &within the input first guess domain.' PRINT '(A,I8,A,A,A)','This is the ',NINT(level),' Pa level of ',field(1:8),'.' PRINT '(A)','Model_domain_outside_first_guess_mask' call exit(1) END IF END IF ilo = MIN(FLOOR(xloc(j,i)) , idim-1) jla = MIN(FLOOR(yloc(j,i)) , jdim-1) dlon = xloc(j,i) - ilo dlat = yloc(j,i) - jla deltalat = 1. deltalon = 1. END IF ! Test #1a: No flag values, this is over land and this is a land field (specifically, ! this is the soil temperature or the soil height). IF ( ( array(ilo ,jla ) .GT. big_negative ) .AND. & ( array(ilo+1,jla ) .GT. big_negative ) .AND. & ( array(ilo ,jla+1) .GT. big_negative ) .AND. & ( array(ilo+1,jla+1) .GT. big_negative ) .AND. & ! ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILT' ) .OR. ( field(1:7) .EQ. 'SOILHGT' ) .OR. & ( field(1:7) .EQ. 'CANOPYM' ) .OR. ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) .OR. & ( field(1:6) .EQ. 'MONALB' ) .OR. ( field(1:8) .EQ. 'ALBSNOMX' ) ) ) THEN IF ( ( ilo .GE. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= ( ( deltalon-dlon ) * ( ( deltalat-dlat ) * array ( ilo , jla ) & + dlat * array ( ilo , jla+1 ) ) & + dlon * ( ( deltalat-dlat ) * array ( ilo+1 , jla ) & + dlat * array ( ilo+1 , jla+1 ) ) ) & / ( deltalat * deltalon ) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= array ( 1 , jla ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .GE. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( ilo , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( 1 , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) END IF ! Test #1b: No flag values, this is over land and this is a land field (specifically, this is ! the soil moisture). A check is set to make sure that no value larger than 0.5 gets into ! the average. ELSE IF ( ( array(ilo ,jla ) .GT. big_negative ) .AND. & ( array(ilo ,jla ) .LE. 0.5 ) .AND. & ( array(ilo+1,jla ) .GT. big_negative ) .AND. & ( array(ilo+1,jla ) .LE. 0.5 ) .AND. & ( array(ilo ,jla+1) .GT. big_negative ) .AND. & ( array(ilo ,jla+1) .LE. 0.5 ) .AND. & ( array(ilo+1,jla+1) .GT. big_negative ) .AND. & ( array(ilo+1,jla+1) .LE. 0.5 ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILM' ) .OR. ( field(1:5) .EQ. 'SOILW' ) ) ) THEN IF ( ( ilo .GE. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= ( ( deltalon-dlon ) * ( ( deltalat-dlat ) * array ( ilo , jla ) & + dlat * array ( ilo , jla+1 ) ) & + dlon * ( ( deltalat-dlat ) * array ( ilo+1 , jla ) & + dlat * array ( ilo+1 , jla+1 ) ) ) & / ( deltalat * deltalon ) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= array ( 1 , jla ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .GE. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( ilo , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( 1 , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) END IF ! Test #1c: This is where the soil moisture is set to 1.0 over the water. ELSE IF ( ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) .AND. & ( ( field(1:5) .EQ. 'SOILM' ) .OR. ( field(1:5) .EQ. 'SOILW' ) ) ) THEN dum2d(j,i) = 1. ! Test #2: No flag values, this is over water and this is a water field. ELSE IF ( ( array(ilo ,jla ) .GT. big_negative ) .AND. & ( array(ilo+1,jla ) .GT. big_negative ) .AND. & ( array(ilo ,jla+1) .GT. big_negative ) .AND. & ( array(ilo+1,jla+1) .GT. big_negative ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) .AND. & ( field(1:6) .EQ. 'SEAICE' ) ) THEN IF ( ( ilo .GE. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= ( ( deltalon-dlon ) * ( ( deltalat-dlat ) * array ( ilo , jla ) & + dlat * array ( ilo , jla+1 ) ) & + dlon * ( ( deltalat-dlat ) * array ( ilo+1 , jla ) & + dlat * array ( ilo+1 , jla+1 ) ) ) & / ( deltalat * deltalon ) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= array ( 1 , jla ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .GE. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( ilo , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( 1 , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) END IF ! Test #3: All flag values, this is over land and this is a land field. ELSE IF ( ( array(ilo ,jla ) .LT. big_negative ) .AND. & ( array(ilo+1,jla ) .LT. big_negative ) .AND. & ( array(ilo ,jla+1) .LT. big_negative ) .AND. & ( array(ilo+1,jla+1) .LT. big_negative ) .AND. & ! ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILT' ) .OR. ( field(1:7) .EQ. 'SOILHGT' ) .OR. & ( field(1:6) .EQ. 'MONALB' ) .OR. ( field(1:8) .EQ. 'ALBSNOMX' ) .OR. & ( field(1:7) .EQ. 'CANOPYM' ) .OR. ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) ) ) THEN ! Check to see how many values of real data are within the surrounding ! 16 point square. close_count = 0 close_data = 0 DO j_close = MAX (1,jla-1) , MIN(jdim,jla+2) DO i_close = MAX (1,ilo-1) , MIN(idim,ilo+2) IF ( array(i_close,j_close) .GT. big_negative ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) END IF END DO END DO ! Test #3a: There is someone close to use as data. Use the average of all ! of the available (yet distant) data. ! But wait, if this is the albedo field, and we are over water, who CARES ! how many points are valid, we set this to a default value of 8%. IF ( ( close_count .GT. 0 ) .AND. ( field(1:6) .EQ. 'MONALB' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 8. ELSE IF ( ( close_count .GT. 0 ) .AND. ( field(1:8) .EQ. 'ALBSNOMX' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 8. ELSE IF ( ( close_count .GT. 0 ) .AND. & ( ( field(1:7) .EQ. 'CANOPYM' ) .OR. ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 0. ELSE IF ( close_count .GT. 0 ) THEN dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Island-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Test #3b: There is nobody close to this grid point, and it is soil temperature. ! There is a predetermined value of soil temperature that is used which is a ! function of latitude. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:5) .EQ. 'SOILT' ) ) THEN IF ( ( landuse(j,i) .NE. bhi(23,1) ) .AND. ( tempgrd_index .GT. 0 ) ) THEN dum2d(j,i) = tempgrd(j,i) + 0.0065 * terrain(j,i) ELSE CALL missing_soil ( yla , dum2d(j,i) , terrain(j,i) ) dum2d(j,i) = dum2d(j,i) + 0.0065 * terrain(j,i) END IF ! Test #3c: There is nobody close to this grid point, and it is the soil ! height. Guess what, our pre-determined value here is just the terrain ! elevation. Why do this? The soil temps are going to be adjusted later based ! on the difference in elevation between the soil height and the model's elevation. ! If we make them identical in this island location, then the adjustment is zero. ! Since the annual-mean, deep-soil temperature is already adjusted to the ! elevation, we are set. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:7) .EQ. 'SOILHGT' ) ) THEN dum2d(j,i) = terrain(j,i) ! This is an island type situation, and we are setting the maximum albedo. ! Very few sub-grid scale islands are covered with snow, so let's just set ! a value of 16% for each month and a max of 60%. If we are over water, ! set the albedo to 8%. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:6) .EQ. 'MONALB' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) ) THEN dum2d(j,i) = 16. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:6) .EQ. 'MONALB' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 8. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:8) .EQ. 'ALBSNOMX' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) ) THEN dum2d(j,i) = 60. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:8) .EQ. 'ALBSNOMX' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 8. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:7) .EQ. 'LANDSEA' ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) ) THEN dum2d(j,i) = 1. ELSE IF ( ( close_count .EQ. 0 ) .AND. & ( ( field(1:7) .EQ. 'CANOPYM' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) ) THEN dum2d(j,i) = 0. ELSE IF ( ( close_count .EQ. 0 ) .AND. & ( ( field(1:7) .EQ. 'CANOPYM' ) .OR. ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) ) THEN dum2d(j,i) = 0. END IF ! Test #3.1: All flag values, this is over land and this is a land field, specifically, ! a test for soil moisture. ELSE IF ( ( ( array(ilo ,jla ) .LE. big_negative ) .OR. ( array(ilo ,jla ) .GE. 0.9 ) ) .AND. & ( ( array(ilo+1,jla ) .LE. big_negative ) .OR. ( array(ilo+1,jla ) .GE. 0.9 ) ) .AND. & ( ( array(ilo ,jla+1) .LE. big_negative ) .OR. ( array(ilo ,jla+1) .GE. 0.9 ) ) .AND. & ( ( array(ilo+1,jla+1) .LE. big_negative ) .OR. ( array(ilo+1,jla+1) .GE. 0.9 ) ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILM' ) .OR. ( field(1:5) .EQ. 'SOILW' ) ) ) THEN ! Check to see how many values of real data are within the surrounding ! 16 point square. close_count = 0 close_data = 0 DO j_close = MAX (1,jla-1) , MIN(jdim,jla+2) DO i_close = MAX (1,ilo-1) , MIN(idim,ilo+2) IF ( ( array(i_close,j_close) .GT. big_negative ) .AND. & ( array(i_close,j_close) .LT. 0.5 ) ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) ELSE IF ( ( array(i_close,j_close) .GE. 0.5 ) .AND. & ( array(i_close,j_close) .LT. 0.9 ) ) THEN close_count = close_count + 1 close_data = close_data + 0.3 END IF END DO END DO ! Test #3.1a: There is someone close to use as data. Use the average of all ! of the available (yet distant) data. IF ( close_count .GT. 0 ) THEN dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Island-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Test #3.1b: There is nobody close to this grid point, and it is soil moisture. ! There is a predetermined value of 0.3 that is used. ELSE IF ( close_count .EQ. 0 ) THEN dum2d(j,i) = 0.3 IF ( rec4%print_mask ) THEN PRINT '(A,I4,A,I4,A)', & 'Warning: Had to put in a bogus moisture value at (j,i) = (',i,',',j,').' END IF END IF ! Test #4: All flag values, this is over water and this is a water field, a test for ! the seaice field. ELSE IF ( ( array(ilo ,jla ) .LT. big_negative ) .AND. & ( array(ilo+1,jla ) .LT. big_negative ) .AND. & ( array(ilo ,jla+1) .LT. big_negative ) .AND. & ( array(ilo+1,jla+1) .LT. big_negative ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) .AND. & ( field(1:6) .EQ. 'SEAICE' ) ) THEN ! Is anyone within an additional grid point? close_count = 0 close_data = 0 DO j_close = MAX (1,jla-1) , MIN(jdim,jla+2) DO i_close = MAX (1,ilo-1) , MIN(idim,ilo+2) IF ( array(i_close,j_close) .GT. big_negative ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) END IF END DO END DO ! Test #4a: There is someone close to use as data. IF ( close_count .GT. 0 ) THEN dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Lake-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Test #4b: There are no close points to use as data for the sea ice. ELSE IF ( ( close_count .EQ. 0 ) .AND. ( field(1:6) .EQ. 'SEAICE' ) ) THEN dum2d(j,i) = 0 IF ( rec4%print_mask ) THEN PRINT '(A,I4,A,I4,A)', & 'Warning: Had to put in a bogus sea ice value at (j,i) = (',i,',',j,').' END IF END IF ! Test #5a: Some flag values, this is over land and this is a land field, for either ! soil temperature or soil height. ELSE IF ( ( ( array(ilo ,jla ) .LT. big_negative ) .OR. & ( array(ilo+1,jla ) .LT. big_negative ) .OR. & ( array(ilo ,jla+1) .LT. big_negative ) .OR. & ( array(ilo+1,jla+1) .LT. big_negative ) ) .AND. & ! ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILT' ) .OR. ( field(1:7) .EQ. 'SOILHGT' ) .OR. & ( field(1:6) .EQ. 'MONALB' ) .OR. ( field(1:8) .EQ. 'ALBSNOMX' ) .OR. & ( field(1:7) .EQ. 'CANOPYM' ) .OR. ( field(1:7) .EQ. 'LANDSEA' ) .OR. & ( field(1:5) .EQ. 'SNOWH' ) .OR. ( field(1:5) .EQ. 'WEASD' ) ) ) THEN ! Count how many real data points exist in the surrounding nearest ! box of 4 points. It has to be at least 1 since there was already a test ! for zero real-data points. close_count = 0 close_data = 0 DO j_close = MAX (1,jla) , MIN(jdim,jla+1) DO i_close = MAX (1,ilo) , MIN(idim,ilo+1) IF ( array(i_close,j_close) .Ge. big_negative ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) END IF END DO END DO dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Coastal-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Test #5b: Some flag values, this is over land and this is a land field, specifically, this is the ! coastal areas for soil moisture. ELSE IF ( ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:5) .EQ. 'SOILM' ) .OR. ( field(1:5) .EQ. 'SOILW' ) ) ) THEN ! Count how many real data points exist in the surrounding nearest ! box of 4 points. It has to be at least 1 since there was already a test ! for zero real-data points. close_count = 0 close_data = 0 DO j_close = MAX (1,jla) , MIN(jdim,jla+1) DO i_close = MAX (1,ilo) , MIN(idim,ilo+1) IF ( ( array(i_close,j_close) .GE. big_negative ) .AND. & ( array(i_close,j_close) .LE. 0.5 ) ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) ELSE IF ( ( array(i_close,j_close) .GT. 0.5 ) .AND. & ( array(i_close,j_close) .LE. 0.9 ) ) THEN close_count = close_count + 1 close_data = close_data + 0.3 END IF END DO END DO dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Coastal-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Test #6: Some flag values, this is over water and this is a water field, specifically seaice. ELSE IF ( ( ( array(ilo ,jla ) .LT. big_negative ) .OR. & ( array(ilo+1,jla ) .LT. big_negative ) .OR. & ( array(ilo ,jla+1) .LT. big_negative ) .OR. & ( array(ilo+1,jla+1) .LT. big_negative ) ) .AND. & ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) .AND. & ( field(1:3) .EQ. 'SEA' ) ) THEN ! Count how many real data points exist in the surrounding nearest ! box of 4 points. It has to be at least 1 since there was already a test ! for zero real-data points for sea ice above. close_count = 0 close_data = 0 DO j_close = MAX (1,jla) , MIN(jdim,jla+1) DO i_close = MAX (1,ilo) , MIN(idim,ilo+1) IF ( array(i_close,j_close) .GT. big_negative ) THEN close_count = close_count + 1 close_data = close_data + array(i_close,j_close) END IF END DO END DO dum2d(j,i) = close_data / REAL ( close_count) IF ( rec4%print_mask ) THEN PRINT '(A,A,A,I4,A,I4,A)', & 'Coastal-type value for ',field(1:8),' found for (j,i)=(',i,',',j,').' END IF ! Mask test version #2: This is shut off to allow soil temps under the sea ice. Who'd'a thunk? ! ! ! Test #7: This is a land field and this is over water. ! ! ELSE IF ( ( ter_data(land_mask_index)%array(j,i) .EQ. 0 ) .AND. & ! ( ( field(1:4) .EQ. 'SOIL' ) ) ) THEN ! ! dum2d(j,i) = 0 ! Test #8: This is a water field and this is over land. ELSE IF ( ( ter_data(land_mask_index)%array(j,i) .EQ. 1 ) .AND. & ( ( field(1:6) .EQ. 'SEAICE' ) ) ) THEN dum2d(j,i) = 0 ELSE PRINT '(A,A,A,2I4,A)','You dropped through the IF tests for the masked interpolation: ',field(1:8),' at (j,i)=',j,i,'.' PRINT '(A)','This is a situation to which attention needs to be given.' PRINT '(A)','Drop_through_mask_interpolation' call exit(1) END IF END DO mask_j_loop END DO mask_i_loop ! If the data do not assume any sort of masking, we can do a 4 point interpolation. ! This guy is considerably easier to follow. ELSE ! Store the latitude and longitude from the TERRAIN file into some ! more attractively named local variables. We can split off the dot and cross ! point variables for the interpolation locations. IF ( ( field(1:8) .EQ. 'U ' ) .OR. & ( field(1:8) .EQ. 'V ' ) ) THEN latitude = ter_data(latitude_dot_index)%array longitude = ter_data(longitude_dot_index)%array icd_old = icd icd = 0 ELSE latitude = ter_data(latitude_cross_index)%array longitude = ter_data(longitude_cross_index)%array icd_old = icd icd = 1 END IF ! For the regular lat/lon grids, the longitude field needs to be ! massaged to allow correct horizontal interpolation across ! Greenwich and the date line. IF ( llflag .EQ. 0 )THEN IF ( deltalon .LT. 0 ) THEN WHERE ( longitude .GT. startlon ) longitude = longitude - 360. END WHERE ELSE IF ( deltalon .GT. 0 ) THEN WHERE ( longitude .LT. startlon ) longitude = longitude + 360. END WHERE END IF END IF ! Loop through each of the model grid points for the unmasked fields. ! pjj/cray !$OMP parallel do default(shared) & !$OMP private(i,j,yla,xlo,ilo,jla,dlon,dlat,xlonc_mod,deltalat,deltalon,wrap,deltalon_save) interpolate_i_loop : DO i = 1, ixf interpolate_j_loop : DO j = 1, jxf ! Calculate values of first guess at the model grid point locations: ! note that i is in west-east (or x) direction, and j in south-north ! (or y) direction. This is a weighted mean in both the x and y ! directions. The interpolation is done in latitude/longitude space ! for a regular latitude longitude projection, and the interpolation ! uses grid point space for a Lambert or polar stereographic projection. ! For the regular latitude longitude projection: ! The latitude/longitude of each of the grid points in the regular ! lat/lon grid is known: (startlat/deltalat, startlon/deltalon). The ! lat/lon of each model grid point has been previously computed. ! For the Lambert and polar stereographic projections: ! The individual latitude and longitude of the eventual output domain ! are used to compute the (i,j) index into the input array for interpolation. ! ! (ilo,jla+1) (ilo+1,jla+1) ! - ------------- ! deltalat-dlat | | | ! | | | ! - | * | ! dlat | | (i,j) | ! | | | ! - ------------- ! (ilo,jla) (ilo+1,jla) ! ! |--|--------| ! dlon deltalon-dlon yla = latitude(j,i) xlo = longitude(j,i) ! Make the interpolation specific to the first guess input projection. For ! a regular latitude longitude projection, the set of values (deltalat, ! deltalon, dlat, dlon, startlat, startlon, xlo and yla) are latitudes and ! longitudes. The size of the above grid box is deltalat by deltalon. For ! the Lambert and polar projections, the set of values (dlat, dlon, deltalat ! and deltalon) are grid point values. The above grid box is then assumed to ! be 1 grid unit on each side. wrap = .FALSE. IF ( llflag .EQ. 0 ) THEN xloc(j,i) = (xlo-startlon)/deltalon + 1. yloc(j,i) = (yla-startlat)/deltalat + 1. ilo = MIN(FLOOR((xlo-startlon)/deltalon) + 1, idim-1) jla = MIN(FLOOR((yla-startlat)/deltalat) + 1, jdim-1) dlon = xlo - ( startlon + (ilo-1) * deltalon ) IF ( ( deltalon .GT. 0 ) .AND. & ( FLOOR((xlo-startlon)/deltalon) + 1 .GE. idim ) .AND. & ( xlo + deltalon .GE. startlon + 360 ) ) THEN IF ( rec4%print_interp) THEN PRINT '(A,I3,A,I3,A,A)','Wrapping at point (',i,',',j,'), variable = ',field(1:8) END IF dlon = xlo - ( startlon + (idim-1) * deltalon ) wrap = .TRUE. ilo = idim ELSE IF ( ( deltalon .LT. 0 ) .AND. & ( FLOOR((xlo-startlon)/deltalon) + 1 .GE. idim ) .AND. & ( xlo + deltalon .LE. startlon - 360 ) ) THEN IF ( rec4%print_interp) THEN PRINT '(A,I3,A,I3,A,A)','Wrapping at point (',i,',',j,'), variable = ',field(1:8) END IF dlon = xlo - ( startlon + (idim-1) * deltalon ) wrap = .TRUE. ilo = idim ELSE IF ( ABS(dlon) .GT. ABS(deltalon) ) THEN IF ( rec4%print_interp) THEN PRINT '(A,I3,A,I3,A,A)','Regular lat/lon horizontal extrapolation at point (',i,',',j,'), variable = ',field(1:8) END IF dlon = 1. END IF dlat = yla - ( startlat + (jla-1) * deltalat ) ELSE IF ( llflag .EQ. 4 ) THEN CALL llxy_gauss( yla , xlo , xloc(j,i) , yloc(j,i) , wrap , & idim , startlat , startlon , deltalon , nlat ) IF ( .NOT. wrap ) THEN ilo = MIN(FLOOR(xloc(j,i)) , idim-1) ELSE ilo = idim END IF jla = MIN(FLOOR(yloc(j,i)) , jdim-1) dlon = xloc(j,i) - ilo dlat = yloc(j,i) - jla deltalat = 1. deltalon_save = deltalon deltalon = 1. ELSE IF ( ( llflag .EQ. 1 ) .OR. ( llflag .EQ. 3 ) .OR. ( llflag .EQ. 5 ) ) THEN IF (xlonc .GT. 180 ) THEN xlonc_mod = xlonc - 360. ELSE IF (xlonc .LT. -180 ) THEN xlonc_mod = xlonc + 360. ELSE xlonc_mod = xlonc END IF IF ( ( .NOT. similar ) .OR. ( icd .NE. icd_old ) ) THEN IF ( ( llflag .EQ. 1 ) .OR. & ( llflag .EQ. 3 ) .OR. & ( llflag .EQ. 5 ) ) THEN CALL latlon_to_ij(proj,yla,xlo,xloc(j,i) , yloc(j,i)) END IF IF ( ( NINT(xloc(j,i)) .GT. idim ) .OR. & ( NINT(yloc(j,i)) .GT. jdim ) .OR. & ( NINT(xloc(j,i)) .LT. 1 ) .OR. & ( NINT(yloc(j,i)) .LT. 1 ) ) THEN PRINT '(A)','The chosen domain from the terrain data is not completely contained & &within the input first guess domain.' PRINT '(A,I8,A,A,A)','This is the ',NINT(level),' Pa level of ',field(1:8),'.' PRINT '(A)','Model_domain_outside_first_guess' call exit(1) END IF END IF ilo = MIN(FLOOR(xloc(j,i)) , idim-1) jla = MIN(FLOOR(yloc(j,i)) , jdim-1) dlon = xloc(j,i) - ilo dlat = yloc(j,i) - jla deltalat = 1. deltalon = 1. END IF IF ( ( rec2%linear_interpolation ) .AND. ( .NOT. wrap ) ) THEN IF ( ( ilo .GE. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= ( ( deltalon-dlon ) * ( ( deltalat-dlat ) * array ( ilo , jla ) & + dlat * array ( ilo , jla+1 ) ) & + dlon * ( ( deltalat-dlat ) * array ( ilo+1 , jla ) & + dlat * array ( ilo+1 , jla+1 ) ) ) & / ( deltalat * deltalon ) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= array ( 1 , jla ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .GE. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( ilo , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( 1 , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) END IF ELSE IF ( ( .NOT. rec2%linear_interpolation ) .AND. ( .NOT. wrap ) ) THEN IF ( ( ilo .GE. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i) = bint ( xloc(j,i) , yloc(j,i) , array , idim , jdim , 0 ) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .GE. 1 ) ) THEN dum2d(j,i)= array ( 1 , jla ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .GE. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( ilo , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) ELSE IF ( ( ilo .LT. 1 ) .AND. ( jla .LT. 1 ) ) THEN dum2d(j,i)= array ( 1 , 1 ) PRINT '(A,I3,A,I3,A,A)','Extrapolating at point (',i,',',j,'), variable = ',field(1:8) END IF ELSE IF ( wrap ) THEN IF ( ilo .EQ. 0 ) THEN PRINT '(A)','Wrapping for ilo=0.' ilo = 1 ELSE IF ( ilo .LT. 0 ) THEN PRINT '(A)','Wrapping troubles for ilo.' PRINT '(A)','Wrapping_woes_ilo' call exit(1) END IF IF ( jla .EQ. 0 ) THEN PRINT '(A)','Wrapping for jla=0.' jla = 1 ELSE IF ( jla .LT. 0 ) THEN PRINT '(A)','Wrapping troubles for jla.' PRINT '(A)','Wrapping_woes_jla' call exit(1) END IF dum2d(j,i)= ( ( deltalon-dlon ) * ( ( deltalat-dlat ) * array ( ilo , jla ) & + dlat * array ( ilo , jla+1 ) ) & + dlon * ( ( deltalat-dlat ) * array ( 1 , jla ) & + dlat * array ( 1 , jla+1 ) ) ) & / ( deltalat * deltalon ) END IF IF ( llflag .EQ. 4 ) THEN deltalon = deltalon_save END IF END DO interpolate_j_loop END DO interpolate_i_loop END IF IF ( rec4%print_interp ) THEN PRINT '(//,A,1X,I8)',field,NINT(level) DO j=jxf,1,-10 PRINT '(10f9.1)',dum2d(j,1:ixf:10) END DO END IF ! Store this 2d data on the model grid. CALL proc_grid_store ( field // ' ' // units // ' ' // description , level , dum2d , ixf , jxf ) ! These are the old projection values. Save them for a comparison. If they ! are repeated, then we do not have to recompute the center latitude ! again. idim_old = idim jdim_old = jdim llflag_old = llflag IF ( ( llflag .EQ. 1 ) .OR. ( llflag .EQ. 3 ) .OR. ( llflag .EQ. 5 ) ) THEN lat1_old = lat1 lon1_old = lon1 deltax_old = deltax deltay_old = deltay xlonc_old = xlonc truelat1_old = truelat1 truelat2_old = truelat2 END IF END DO this_fg_file_only ! CLOSE this file. We can reuse this UNIT again if required. CLOSE (UNIT=first_guess_file) IF ( rec4%print_file ) THEN PRINT '(A,A/)' , 'CLOSE file: ',TRIM(file_name) END IF END DO diff_fg_files ! Is this array still allocated? If so, deallocate it at your earliest convenience. IF ( ALLOCATED ( array ) ) THEN DEALLOCATE (array) END IF ! 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.' PRINT '(A)','No_files_to_open' call exit(1) END IF ! 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 ( rec4%print_debug ) 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 ( rec4%print_debug ) 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(landuse) DEALLOCATE(tempgrd) DEALLOCATE(latitude) DEALLOCATE(longitude) DEALLOCATE(xloc) DEALLOCATE(yloc) END SUBROUTINE proc_ingest_first_guess