!-------------------------------------------------------------------- !Name: ESI_composite_mod ! !Type: F90 module ! !Function: ! This module contains major subroutines to compute ! the composite ESI for 2, 4, 8, 12 weeks. ! !Description: ! The computation is to read in all the smoothed ESI from ! previous days (84 or 364 for 1 year total) and generate ! the average ESI values for the time period. The average ! process will exclude snow days and days with fill ESI ! values. ! Compute the 1 week composite ! US_FPPM_SMN_1WK_2015175.dat ! US_FPPM_RAW_1WK_2015175.dat ! Compute the 2/4/8/12 week anomalies ! Inputs are ! AVG_SMN_2000175.dat ! SD_SMN_2000175.dat ! US_DFPPM_2WK_2015175.dat ! US_DFPPM_4WK_2015175.dat ! US_DFPPM_8WK_2015175.dat ! US_DFPPM_12WK_2015175.dat ! Computes a mask that will flag any pixels in the ! 4-week composite that is older than 14 days ! LAST_DAY_2015175.dat ! NCLEAR_2015175.dat ! MASK_2015175.dat ! !Files Needed: ! composite_ESI8km.cfg - configuration file ! ESI[YYYY][DDD].bin - smoothed ESI inputs ! !Modules Needed: ! - type_kinds ! - GETD_constant_mod ! - GETD_log_mod ! - GETD_constant_mod ! !Subroutines Contained: ! - N/A ! !Calling Sequence: ! N/A ! !Inputs: ! N/A ! !Outputs: ! N/A ! !System Calls: ! N/A ! !History: ! 2014-09-25 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS ! 2015-08-25 Added 1 week composite and anomalies computation ! 2016-02-15 Added seasonal snow mask for the ESI ! results and the QC flags. Used ET QC flag ! from ET results and added the seasonal mask in. !-------------------------------------------------------------------- MODULE ESI_composite_mod USE type_kinds USE GETD_constants_mod USE GETD_log_mod IMPLICIT NONE !By default, all the functions and variables are only to be used within this !moduel PRIVATE !Declare following subroutines to be accessed by outside programs PUBLIC :: ESI_composite PRIVATE :: read_composite_config PRIVATE :: compute_weekly_anom ! Input for composite ESI data should contain ! the data directory and the data name conventions ! The input data will follow the format ! Note the name has day of the year information. ! Added the seasonal mask info. TYPE composite_ESI_type CHARACTER(LONG_STR) :: in_esi_ave_sd_dir CHARACTER(LONG_STR) :: in_esi_ave_header CHARACTER(LONG_STR) :: in_esi_sd_header CHARACTER(LONG_STR) :: in_esi_dir CHARACTER(LONG_STR) :: in_esi_file_name CHARACTER(LONG_STR) :: in_et_qc_file_name CHARACTER(LONG_STR) :: in_mask_dir CHARACTER(LONG_STR) :: in_mask_file_name ! Output reference ET and ESI files CHARACTER(LONG_STR) :: out_esi_dir CHARACTER(LONG_STR) :: out_esi_comp1_file_name CHARACTER(LONG_STR) :: out_esi_comp2_file_name CHARACTER(LONG_STR) :: out_esi_comp4_file_name CHARACTER(LONG_STR) :: out_esi_comp8_file_name CHARACTER(LONG_STR) :: out_esi_comp12_file_name CHARACTER(LONG_STR) :: out_esi_qc_file_name ! Spatial information for the output domain INTEGER(Long) :: in_row_num, in_col_num REAL(Single) :: in_east_lon, in_west_lon REAL(Single) :: in_north_lat, in_south_lat ! REAL(Single) :: in_lat_resolution REAL(Single) :: in_lon_resolution REAL(Single) :: in_fill_value END TYPE composite_ESI_type CONTAINS !=============================================================== ! Name: ESI_composite ! ! Type: Subroutine ! ! Function: ! Compute the composite ESI from daily smoothed ESI. ! ! Descriptions: ! Read in the daily ESI file, ! 1. Create 1 week composite fPET using smoothed inputs. ! 2. Compute the anamoly using MEAN and SD from restrospective data ! for the 2, 4, 8, 12 weekly ESI. ! ! Files Needed: ! - configure file ! - Input ESI data file ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL ESI_composite( in_config_file, & ! in_year, in_doy, return_status ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_config_file I Name of the input config ! in_year I Year ! in_doy I Day of the year ! ! Output: ! Name Type Description ! --------------------------------------------------- ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! Note: ! 2015-09-16 ! Add the 1 week composite and anaomoly computation ! separate the anamoly computation into a new function. !=============================================================== SUBROUTINE ESI_composite( in_config_file, & in_year, in_doy, return_status) USE file_io_mod USE matrix_mod USE sunrise_mod USE GETD_log_mod IMPLICIT NONE ! --------- ! Arguments ! --------- CHARACTER(LONG_STR), INTENT( in) :: in_config_file INTEGER(Long), INTENT( in) :: in_year, in_doy INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j, k, iday INTEGER(Long) :: col, row INTEGER(Long) :: start_year, start_doy INTEGER(Long) :: cur_year INTEGER(Long) :: cur_doy, edoy INTEGER(Long) :: max_days INTEGER(Long) :: istat, icount TYPE(composite_ESI_type) :: ESI CHARACTER(LONG_STR) :: file_name, path_name ! Input matrix holding daily fPET for 1 week composite ESI REAL(Single), ALLOCATABLE, DIMENSION(:,:,:) :: in_ESI_daily ! Temp matrix for ESI computing REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ESI_sum INTEGER(Long), ALLOCATABLE, DIMENSION(:,:) :: ESI_count ! Weekly composited ESI REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_comp1 ! Read in the configuration return_status = SUCCESS CALL read_composite_config(in_config_file, ESI, istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Read input config Error' return_status = FAILURE RETURN END IF ! Init the matrix with the information from the input ! We only need 12 weeks maximum, so set this to 84 ! Changed this to 1 week composite only max_days = 7 ALLOCATE(in_ESI_daily(ESI%in_col_num, ESI%in_row_num, max_days),& stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI inputs error' return_status = FAILURE RETURN END IF ALLOCATE(ESI_sum(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI SUM error' return_status = FAILURE RETURN END IF ALLOCATE(ESI_count(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI count error' return_status = FAILURE RETURN END IF ALLOCATE(out_ESI_comp1(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI outputs error' return_status = FAILURE RETURN END IF ! Init the input ESI matrix will fill data ! in_ESI_daily = ESI%in_fill_value init_if: IF ( return_status /= SUCCESS ) THEN WRITE(GETD_logunit, *)'Error in allocating the ESI space,', & 'Program exit error in', ' ESI_composite' ELSE ! Loop to read in all the ESI data, using the information ! from the configuration file to set the input data path WRITE(GETD_logunit, *) 'Composite on ', in_year, in_doy ! Decide the start composite day information for 12 weeks call get_start_yyyydoy_by_yyyydoy(in_year, in_doy, max_days, & start_year, start_doy, istat) write(GETD_logunit,*) 'Start from', start_year, start_doy, & ' to ', in_year, in_doy ! if the start day is valid, continue to open the files ! from the start day and read in all the data. valid_start_if: if ( istat == SUCCESS ) then ! Count how many days has valid input iday = 0 ! Loop variable to track the index i = 1 year_loop: DO cur_year = start_year, in_year ! Decide the start doy in the current year if ( cur_year > start_year ) then cur_doy = 1 else cur_doy = start_doy endif ! Decide the end doy in the current year if ( cur_year < in_year ) then if ( is_leap_year(cur_year) ) then edoy = 366 else edoy = 365 endif else edoy = in_doy endif ! Open the files during the period to loadinthe ! fPET daily data. doy_loop: DO cur_doy = cur_doy, edoy WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%in_esi_dir), '/', cur_year, cur_doy, & '/', TRIM(ESI%in_esi_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ESI_daily(:,:,i), istat) if ( istat == SUCCESS ) then WRITE(GETD_logunit,*)'Read daily fPET ',& trim(file_name), ' successfully.' iday = iday + 1 endif i = i + 1 end do doy_loop END DO year_loop ! If input ESI has more than one day's data ! continue to do the compositing process. write(GETD_logunit,*)'Total ', iday, ' read in ',& edoy - cur_doy + 1 load_ESI_if: IF ( iday > 0 ) THEN !reset the values ESI_sum = 0 ! Not ESI%in_fill_value, we will fill in later ESI_count = 0 out_ESI_comp1 = ESI%in_fill_value ! The composite process only counts the ESI values that are ! between 0 and 5, no snow cover, and no filled value. ! Note, since we already used the snow mask in the ET ! computation, we don't need to check snow again. t_loop: DO k = max_days, 1, -1 out_col_loop: DO i = 1, ESI%in_col_num out_row_loop: DO j = 1, ESI%in_row_num IF ( in_ESI_daily(i,j, k) > 0 .and. & in_ESI_daily(i,j, k) < 5. ) THEN ESI_sum(i,j) = ESI_sum(i,j) + in_ESI_daily(i,j, k) ESI_count(i,j) = ESI_count(i,j) + 1 ENDIF END DO out_row_loop END DO out_col_loop END DO t_loop ! 1 weeks composite DO i = 1, ESI%in_col_num DO j = 1, ESI%in_row_num IF ( ESI_count(i,j) > 0 ) THEN out_ESI_comp1(i,j) = ESI_sum(i,j) / ESI_count(i,j) ENDIF END DO END DO ! Write the 1 week composite to the output direcotry ! Note this step create an additional level WRITE(path_name,'(a, a1, i4.4, i3.3)') & TRIM(ESI%out_esi_dir), '/',in_year, in_doy CALL SYSTEM('mkdir -p ' // TRIM(path_name)) WRITE(file_name,'(a,a1, a)') & TRIM(path_name), & '/', TRIM(ESI%out_esi_comp1_file_name) ! Write(GETD_logunit,*)'Write to ', trim(file_name) CALL write_out_2d_real(out_ESI_comp1, ESI%in_col_num, & ESI%in_row_num, file_name, istat) IF ( ALLOCATED(in_ESI_daily) ) DEALLOCATE(in_ESI_daily) IF ( ALLOCATED(ESI_sum) ) DEALLOCATE(ESI_sum) IF ( ALLOCATED(ESI_count) ) DEALLOCATE(ESI_count) IF ( ALLOCATED(out_ESI_comp1)) DEALLOCATE(out_ESI_comp1) ENDIF load_ESI_if ENDIF valid_start_if ENDIF init_if !2. Check if the input ESI AVE and SD parameters exist ! If both exists, continue to compute the anamoly ! otherwise, return with warning ! The input ESI AVE and SD file names should follow ! the format as ! SD_SMN_2001_2014241.dat ! AVG_SMN_2001_2014241.dat ! !3. Compute the 2, 4, 8, 12 weekly ! Note this anamoly only applies to ! the pixels without seasonal snow mask. ! WRITE(GETD_logunit,*)'Compute anamoly with seasonal mask on', & in_year, in_doy CALL compute_weekly_anom( ESI, & in_year, in_doy, istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Compute anamoly for ESI error!' return_status = FAILURE END IF END SUBROUTINE ESI_composite !=============================================================== ! Name: read_composite_config ! ! Type: Subroutine ! ! Function: ! Read in the spatial information of composite data sets. ! ! Description: ! N/A ! Files Needed: ! - in_file_name ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! call read_composite_config(in_file_name, ESI, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_file_name I Input configure file name ! ESI I/O Input ESI data ! ! Output: ! Name Type Description ! --------------------------------------------------- ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! None ! ! History: ! 2014-09-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS ! 2015-08-25 Updated with new configuration file format ! !=============================================================== SUBROUTINE read_composite_config(in_file_name, ESI, return_status) ! !USES:: USE GETD_constants_mod USE file_io_mod USE GETD_log_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- CHARACTER(LONG_STR), INTENT(in) :: in_file_name TYPE(composite_ESI_type), INTENT(out) :: ESI INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- INTEGER(Long) :: file_status, i, input_unit INTEGER(Long) :: istat CHARACTER(LONG_STR) one_line CHARACTER(LONG_STR) str_temp CHARACTER(len=256) :: env_var_value INTEGER :: env_var_len, env_var_status input_unit = GETD_get_next_unit_number() env_var_len = len("HOME_DIR") CALL GET_ENVIRONMENT_VARIABLE("HOME_DIR", env_var_value, env_var_len, env_var_status) OPEN(input_unit, & IOSTAT = file_status, & file = in_file_name, & form = 'formatted', & status = 'old') IF (file_status /= SUCCESS) THEN WRITE(GETD_logunit,*)'Error in opening the ESI configuration file' return_status = FAILURE RETURN ELSE ! Read in the information in the configuration file DO i = 1, 5 READ(input_unit,10) one_line !Skip the commented line END DO ! Input ESI file domain information READ(input_unit,25) ESI%in_east_lon READ(input_unit,25) ESI%in_west_lon READ(input_unit,25) ESI%in_lon_resolution READ(input_unit,25) ESI%in_north_lat READ(input_unit,25) ESI%in_south_lat READ(input_unit,25) ESI%in_lat_resolution READ(input_unit,15) ESI%in_col_num READ(input_unit,15) ESI%in_row_num READ(input_unit,25) ESI%in_fill_value ! READ(input_unit,10) ESI%in_esi_ave_sd_dir READ(input_unit,10) str_temp ESI%in_esi_ave_sd_dir = trim(env_var_value) // str_temp READ(input_unit,10) ESI%in_esi_ave_header READ(input_unit,10) ESI%in_esi_sd_header ! Input daily ET computed in ESI direcotry ! This ET has no RC correction. !READ(input_unit,10) ESI%in_esi_dir READ(input_unit,10) str_temp ESI%in_esi_dir = trim(env_var_value) // str_temp READ(input_unit,10) ESI%in_esi_file_name READ(input_unit,10) ESI%in_et_qc_file_name !READ(input_unit,10) ESI%in_mask_dir READ(input_unit,10) str_temp ESI%in_mask_dir = trim(env_var_value) // str_temp READ(input_unit,10) ESI%in_mask_file_name ! Output reference ET and ESI files READ(input_unit,10) one_line !comment line !READ(input_unit,10) ESI%out_esi_dir READ(input_unit,10) str_temp ESI%out_esi_dir = trim(env_var_value) // str_temp READ(input_unit,10) ESI%out_esi_comp1_file_name READ(input_unit,10) ESI%out_esi_comp2_file_name READ(input_unit,10) ESI%out_esi_comp4_file_name READ(input_unit,10) ESI%out_esi_comp8_file_name READ(input_unit,10) ESI%out_esi_comp12_file_name READ(input_unit,10) ESI%out_esi_qc_file_name END IF return_status = SUCCESS CLOSE(input_unit) CALL GETD_release_unit_number(input_unit) !formats to read in the inputs 10 FORMAT (30x, a300) !format to read in the string after 30 spaces 15 FORMAT (30x, i10) !format to read in the integer number, one per line 20 FORMAT (30x, f7.2) !format to read in the real value, one per line 25 FORMAT (30x, e12.3) 30 FORMAT (30x, a4) !format to read in the string after 30 spaces END SUBROUTINE read_composite_config !=============================================================== ! Name: compute_weekly_anom ! ! Type: Subroutine ! ! Function: ! Compute the composite ESI from daily smoothed ESI. ! ! Descriptions: ! Used the 1 week composite ESI value as inputs and ! retrospective data (Mean and SD) to compute the anamoly ! for multiple weeks. ! It will read in the input data files in multiple weeks ! to compute the 2, 4, 8, 12 weeks ESI. ! ! Files Needed: ! - configure file ! - Input 1 week composite fPET data files for the ! past 12 weeks. ! - Input directory and file headers for mean and sd ! to compute the anomoly in the ! static data directory. ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL compute_weekly_anom( ESI , & ! in_year, in_doy, return_status ) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! ESI I Input data structure holding ! the static climatology AVE and SD. ! in_year I Year ! in_doy I Day of the year ! ! Output: ! Name Type Description ! --------------------------------------------------- ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! Note: ! For leap year, we use 365 instead of 366 in loading ! the MEAN and SD from the restrospective data. This is due to ! the small number of 366 in the restrospective data. ! !=============================================================== SUBROUTINE compute_weekly_anom( ESI, & in_year, in_doy, & return_status) USE file_io_mod USE matrix_mod USE GETD_log_mod IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(composite_ESI_type), INTENT(in) :: ESI INTEGER(Long), INTENT( in) :: in_year, in_doy INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j, k, iday INTEGER(Long) :: col, row INTEGER(Long) :: cur_year, cur_doy INTEGER(Long) :: istat, icount INTEGER(bByte) :: qc_byte CHARACTER(LONG_STR) :: file_name ! Input Weekly composited ESI and ! mean and sd from retrospective data REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_ESI_comp1 REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_ESI_ave REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_ESI_sd REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_ESI_mask ! Temp matrix for ESI computing REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ESI_comp REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ESI_sd REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ESI_ave INTEGER(Long), ALLOCATABLE, DIMENSION(:,:) :: ESI_count ! Output anomoly matrix REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_comp2 REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_comp4 REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_comp8 REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_comp12 REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_ESI_qc ! Allocate the memory for input matrix ALLOCATE(in_ESI_comp1(ESI%in_col_num, ESI%in_row_num),& stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI weekly inputs error' return_status = FAILURE RETURN END IF ALLOCATE(in_ESI_ave(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI ave error' return_status = FAILURE RETURN END IF ALLOCATE(in_ESI_sd(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI sd error' return_status = FAILURE RETURN END IF ALLOCATE(in_ESI_mask(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI mask error' return_status = FAILURE RETURN END IF ! Allocate the memory for temp matrix ALLOCATE(ESI_comp(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI temp comp error' return_status = FAILURE RETURN END IF ALLOCATE(ESI_ave(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI temp ave error' return_status = FAILURE RETURN END IF ALLOCATE(ESI_sd(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI temp sd error' return_status = FAILURE RETURN END IF ALLOCATE(ESI_count(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI temp count error' return_status = FAILURE RETURN END IF ! Allocate the memory for output matrix ALLOCATE(out_ESI_comp2(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI output2 error' return_status = FAILURE RETURN END IF ALLOCATE(out_ESI_comp4(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI output4 error' return_status = FAILURE RETURN END IF ALLOCATE(out_ESI_comp8(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI output8 error' return_status = FAILURE RETURN END IF ALLOCATE(out_ESI_comp12(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI output 12 error' return_status = FAILURE RETURN END IF ALLOCATE(out_ESI_qc(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI QC error' return_status = FAILURE RETURN END IF out_ESI_comp2 = ESI%in_fill_value out_ESI_comp4 = ESI%in_fill_value out_ESI_comp8 = ESI%in_fill_value out_ESI_comp12 = ESI%in_fill_value out_ESI_qc = 0 ESI_comp = 0 ESI_sd = 0 ESI_ave = 0 ESI_count = 0 return_status = SUCCESS ! Multi weekly loop weekly_loop: DO k = 1, 12 ! Find the right file for the doy ! for all 12 weeks composition cur_year = in_year ! Get the day of the 1 weekly composite cur_doy = in_doy - (7 * (k - 1)) if ( cur_doy < 1 ) then cur_doy = 365 + cur_doy cur_year = in_year - 1 else if ( cur_doy > 365 ) then ! For leap year 366, use 365 instead cur_doy = 365 endif ! Read in the 1 week composite for each input week ! Note this is the output file from previous composite ! step. WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', cur_year, cur_doy, & '/', TRIM(ESI%out_esi_comp1_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ESI_comp1(:,:), istat) if ( istat /= SUCCESS) then WRITE(GETD_logunit, *) 'Warning, missing', & trim(ESI%out_esi_comp1_file_name) endif ! Read in the MEAN and SD files, They should ! be in the position and follow the name convention. WRITE(file_name,'(a,a1,a,i3.3,a4)') & TRIM(ESI%in_esi_ave_sd_dir), '/', & TRIM(ESI%in_esi_ave_header), cur_doy, '.dat' CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ESI_ave(:,:), istat) if ( istat /= SUCCESS ) then WRITE(GETD_logunit,*)'Error: Read ESI AVE for day ',& cur_doy, ' failed.' return_status = FAILURE return endif WRITE(file_name,'(a,a1,a,i3.3,a4)') & TRIM(ESI%in_esi_ave_sd_dir), '/', & TRIM(ESI%in_esi_sd_header), cur_doy, '.dat' CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ESI_sd(:,:), istat) if ( istat /= SUCCESS ) then WRITE(GETD_logunit,*)'Read ESI SD for day ',& cur_doy, ' failed.' return_status = FAILURE return ENDIF ! Loop through the domain and only compute the ! valid input 1 week composite fPET do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num if ( & abs(in_ESI_comp1(j,i) - ESI%in_fill_value) > 1e-6 .and. & abs(in_ESI_ave(j,i) - ESI%in_fill_value) > 1e-6 .and. & abs(in_ESI_sd(j,i) - ESI%in_fill_value) > 1e-6 & ) then ESI_comp(j,i) = ESI_comp(j,i) + in_ESI_comp1(j,i) ESI_ave(j,i) = ESI_ave(j,i) + in_ESI_ave(j,i) ESI_sd(j,i) = ESI_sd(j,i) + in_ESI_sd(j,i) ESI_count(j,i) = ESI_count(j,i) + 1 end if enddo enddo ! Compute the composition when certain weeks met ! 2 weeks if ( k == 2 ) then do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num if ( ESI_count(j,i) > 0 .and. & ESI_comp(j,i) > 0 ) then out_ESI_comp2(j,i) = ( & ESI_comp(j,i)/ESI_count(j,i) & - ESI_ave(j,i)/ESI_count(j,i) )& /(ESI_sd(j,i) / ESI_count(j,i)) if( (out_ESI_comp2(j,i) >= -3.5) .and. (out_ESI_comp2(j,i) <= 3.5)) then out_ESI_qc(j,i) = 0 else out_ESI_comp2(j,i) = ESI%in_fill_value out_ESI_qc(j,i) = 1 endif else ! No valid data in the past 2 weeks out_ESI_qc(j,i) = 1 endif enddo enddo endif ! 4 weeks if ( k == 4 ) then do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num if ( ESI_count(j,i) > 0 .and. & ESI_comp(j,i) > 0 ) then out_ESI_comp4(j,i) = ((ESI_comp(j,i))/ESI_count(j,i) & - ESI_ave(j,i) / ESI_count(j,i) ) / & (ESI_sd(j,i) / ESI_count(j,i)) if( (out_ESI_comp4(j,i) >= -3.5) .and. (out_ESI_comp4(j,i) <= 3.5)) then else out_ESI_comp4(j,i) = ESI%in_fill_value endif endif enddo enddo endif ! 8 weeks if ( k == 8 ) then do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num if ( ESI_count(j,i) > 0 .and. & ESI_comp(j,i) > 0 ) then out_ESI_comp8(j,i) = ((ESI_comp(j,i))/ESI_count(j,i) & - ESI_ave(j,i) / ESI_count(j,i) ) / & (ESI_sd(j,i) / ESI_count(j,i)) if( (out_ESI_comp8(j,i) >= -3.5) .and. (out_ESI_comp8(j,i) <= 3.5)) then else out_ESI_comp8(j,i) = ESI%in_fill_value endif endif enddo enddo endif ! 12 weeks if ( k == 12 ) then do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num if ( ESI_count(j,i) > 0 .and. & ESI_comp(j,i) > 0 ) then out_ESI_comp12(j,i) = ((ESI_comp(j,i))/ESI_count(j,i) & - ESI_ave(j,i) / ESI_count(j,i) ) / & (ESI_sd(j,i) / ESI_count(j,i)) if( (out_ESI_comp12(j,i) >= -3.5) .and. (out_ESI_comp12(j,i) <= 3.5)) then else out_ESI_comp12(j,i) = ESI%in_fill_value endif endif enddo enddo endif ENDDO weekly_loop 100 format(1x, i4,i3, 4F10.4, i10) ! Added the seasonal mask for the final composite data ! ! Read in the ET QC flags to use to produce the output QC ! Read in the seasonal mask and put in the bit flag. WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', cur_year, cur_doy, & '/', TRIM(ESI%in_et_qc_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, out_ESI_qc, istat) if ( istat /= SUCCESS) then WRITE(GETD_logunit, *) 'Warning, missing', & trim(file_name) endif ! Read ESI seasonal MASK, notice that the seasonal ! mask is upside down comparing with the input ET ! data WRITE(file_name,'(a,a1,a,i3.3,a4)') & TRIM(ESI%in_mask_dir), '/', & TRIM(ESI%in_mask_file_name), in_doy, '.dat' !print*, "finalmask: ", trim(file_name) CALL read_2d_real_by_col_row_reverse(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ESI_mask(:,:), istat) if ( istat /= SUCCESS ) then WRITE(GETD_logunit,*)'Error! Read ESI seasonal mask for day ',& cur_doy, ' failed.' return_status = FAILURE return ENDIF ! The ESI seasonal mask is defined as ! 1.0 is the valid pixel ! 0.0 is the pixel need to mask out ! -9999. is the fill value do j = 1, ESI%in_col_num do i = 1, ESI%in_row_num ! clear the pixel when not the filled value ! and not a valid pixe (1.0) if ( abs(in_ESI_mask(j,i) - 1.0) > 1e-3 & .and. (in_ESI_mask(j,i) - ESI%in_fill_value) > 1e-3 ) then out_ESI_comp2(j,i) = ESI%in_fill_value out_ESI_comp4(j,i) = ESI%in_fill_value out_ESI_comp8(j,i) = ESI%in_fill_value out_ESI_comp12(j,i) = ESI%in_fill_value ! Set the seasonal flag qc_byte = INT(out_ESI_qc(j,i)) qc_byte = IBSET(qc_byte, 7) out_ESI_qc(j,i) = REAL(qc_byte) endif enddo enddo ! Write out the results write_out_if: IF ( return_status == SUCCESS ) THEN ! Write ESI ! WRITE(GETD_logunit,*)'Write output composite ESI to', & ! TRIM(ESI%out_esi_dir) ! 2 weekly composite WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', in_year, in_doy, & '/', TRIM(ESI%out_esi_comp2_file_name) CALL write_out_2d_real(out_ESI_comp2, ESI%in_col_num, & ESI%in_row_num, file_name, istat) ! 4 weekly composite WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', in_year, in_doy, & '/', TRIM(ESI%out_esi_comp4_file_name) CALL write_out_2d_real(out_ESI_comp4, ESI%in_col_num, & ESI%in_row_num, file_name, istat) ! 8 weekly composite WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', in_year, in_doy, & '/', TRIM(ESI%out_esi_comp8_file_name) CALL write_out_2d_real(out_ESI_comp8, ESI%in_col_num, & ESI%in_row_num, file_name, istat) ! 12 weekly composite WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', in_year, in_doy, & '/', TRIM(ESI%out_esi_comp12_file_name) CALL write_out_2d_real(out_ESI_comp12, ESI%in_col_num, & ESI%in_row_num, file_name, istat) ! ESI QC WRITE(file_name,'(a,a1,i4.4,i3.3,a1, a)') & TRIM(ESI%out_esi_dir), '/', in_year, in_doy, & '/', TRIM(ESI%out_esi_qc_file_name) CALL write_out_2d_real(out_ESI_qc, ESI%in_col_num, & ESI%in_row_num, file_name, istat) END IF write_out_if IF ( ALLOCATED(in_ESI_comp1)) DEALLOCATE(in_ESI_comp1) IF ( ALLOCATED(in_ESI_ave)) DEALLOCATE(in_ESI_ave) IF ( ALLOCATED(in_ESI_sd)) DEALLOCATE(in_ESI_sd) IF ( ALLOCATED(in_ESI_mask)) DEALLOCATE(in_ESI_mask) IF ( ALLOCATED(ESI_comp)) DEALLOCATE(ESI_comp) IF ( ALLOCATED(ESI_sd)) DEALLOCATE(ESI_sd) IF ( ALLOCATED(ESI_ave)) DEALLOCATE(ESI_ave) IF ( ALLOCATED(ESI_count)) DEALLOCATE(ESI_count) IF ( ALLOCATED(out_ESI_comp2)) DEALLOCATE(out_ESI_comp2) IF ( ALLOCATED(out_ESI_comp4)) DEALLOCATE(out_ESI_comp4) IF ( ALLOCATED(out_ESI_comp8)) DEALLOCATE(out_ESI_comp8) IF ( ALLOCATED(out_ESI_comp12)) DEALLOCATE(out_ESI_comp12) IF ( ALLOCATED(out_ESI_qc)) DEALLOCATE(out_ESI_qc) END SUBROUTINE compute_weekly_anom END MODULE ESI_composite_mod !EOP