!-------------------------------------------------------------------- !Name: ESI_smooth_mod ! !Type: F90 module ! !Function: ! This module smooth the ESI values from the input files ! and then write the results to multiple output files. ! !Description: ! There are four components to compute the smooth from ET data ! 1. Compute the daily reference ET using climatology input solar radiation. ! 2. Compute the daily smooth by dividing the ET with reference ET. ! 3. Smooth the ! 4. Compute the composite smooth for 2, 4, 8 and 12 weeks. ! Note, this process doesn't handle any resolution changes, the input ! and output data will have the same resolution. ! !Files Needed: ! smooth.cfg - configuration file ! snow_mask_yyyyddd.bin - Input snow mask time series ! ESI_yyyyddd.bin - Input ESI time series ! !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-03-25 Changed the input file names and directory. ! Instead of all the inputs/outputs in a single directory, each day ! has its own directory. The program automatically create the path ! by using the input time period information. !-------------------------------------------------------------------- MODULE ESI_smooth_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 :: smooth_ESI PUBLIC :: smooth_point PRIVATE :: read_smooth_config PRIVATE :: read_in_files PRIVATE :: write_out_files PRIVATE :: smooth_region PRIVATE :: filter_series1 PRIVATE :: filter_series2 PRIVATE :: smooth_series PRIVATE :: fill_series PRIVATE :: smooth_init PRIVATE :: smooth_finalize INTEGER(Long), PARAMETER :: ibuff = 50 ! parameter(ibuff=50) !utl subroutines to run the smooth algorithm ! The INPUTS configs TYPE smooth_ESI_type CHARACTER(LONG_STR) :: in_esi_path CHARACTER(LONG_STR) :: in_esi_file_name CHARACTER(LONG_STR) :: in_snowmask_path CHARACTER(LONG_STR) :: in_snowmask_file_name ! Output reference ET and smooth files CHARACTER(LONG_STR) :: out_esi_path CHARACTER(LONG_STR) :: out_esi_file_name INTEGER(Long) :: start_year, start_doy INTEGER(Long) :: end_year, end_doy INTEGER(Long) :: in_day_num INTEGER(Long) :: ifill ! 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 ! corrected temperature at time 1 and 2 REAL(Single), ALLOCATABLE, DIMENSION(:,:,:) :: in_esi_array REAL(Single), ALLOCATABLE, DIMENSION(:,:,:) :: in_snowmask_array REAL(Single), ALLOCATABLE, DIMENSION(:,:,:) :: out_esi_array END TYPE smooth_ESI_type CONTAINS !=============================================================== ! Name: smooth_ESI ! ! Type: Subroutine ! ! Function: ! Main procedure for smoothing multiple input ESI data files and ! output the smoothed files. ! ! Descriptions: ! Reads in the ESI input files and ! applies Savitsky-Golay smoothing function to time series ! on ALEXI grid and writes out time-smoothed daily files. ! ! Files Needed: ! - configure file ! - Input ESI data file ! - Input Snow Mask files ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL smooth_ESI( in_config_file, & ! start_year, start_doy, & ! end_year, end_doy, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_config_file I Name of the input config ! start_year I Start Year ! start_doy I Start day of the year ! end_year I End Year ! end_doy I End day of the year ! ! ! Output: ! Name Type Description ! --------------------------------------------------- ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! Add check of the input day length to prevent !=============================================================== SUBROUTINE smooth_ESI( in_config_file, & start_year, start_doy, & end_year, end_doy, return_status) USE file_io_mod USE matrix_mod IMPLICIT NONE ! --------- ! Arguments ! --------- CHARACTER(LONG_STR), INTENT(in) :: in_config_file INTEGER(Long), INTENT( in) :: start_year INTEGER(Long), INTENT( in) :: start_doy INTEGER(Long), INTENT( in) :: end_year INTEGER(Long), INTENT( in) :: end_doy INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j INTEGER(Long) :: istrt,iend INTEGER(Long) :: istat, icount INTEGER :: cur_col, cur_row CHARACTER(LONG_STR) :: file_name TYPE(smooth_ESI_type) :: smooth ! Set the configuration using the inputs smooth%start_year = start_year smooth%end_year = end_year smooth%start_doy = start_doy smooth%end_doy = end_doy ! Compute the number of days from the input data length IF ( smooth%start_year == smooth%end_year ) THEN smooth%in_day_num = smooth%end_doy - smooth%start_doy + 1 ELSE !Only need to consider the start year IF ( MOD(smooth%start_year, 4) == 0 .AND. & MOD(smooth%start_year, 400)/=0 ) THEN smooth%in_day_num = smooth%end_doy - smooth%start_doy + 1 & + (smooth%end_year - smooth%start_year) * 365 + 1 ELSE smooth%in_day_num = smooth%end_doy - smooth%start_doy + 1 & + (smooth%end_year - smooth%start_year) * 365 ENDIF ENDIF !Note, the smoothing process can only handle data less than !4 years long, not more than 4 years num_days_if: if ( smooth%in_day_num > 1400 ) then return_status = FAILURE else return_status = SUCCESS CALL read_smooth_config(in_config_file, smooth, istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Read input config Error' return_status = FAILURE RETURN END IF WRITE(GETD_logunit,*)'From ', smooth%start_year, smooth%start_doy, & smooth%end_year, smooth%end_doy WRITE(GETD_logunit,*)'total days', smooth%in_day_num ! Init the matrix using input configuration data CALL smooth_init(smooth, istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory error!' return_status = FAILURE RETURN ELSE WRITE(GETD_logunit,*)'Init done, reading files' END IF ! Read in the ESI data CALL read_in_files(smooth, istat) read_in_if: IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit, *)'Error in Reading in the ESI/Snow Mask,', & 'Program exit error in', 'smooth_ESI' return_status = FAILURE ELSE ! Smoothing the input ESI data ! Smooth time series WRITE (GETD_logunit,*) "Smoothing time series..." istrt = smooth%start_doy iend = smooth%end_doy CALL smooth_region(smooth, istat) ! ! write out the results to the output files ! if the computation is successful write_out_if: IF ( return_status == SUCCESS ) THEN WRITE(GETD_logunit,*)'Write output smooth files' call write_out_files(smooth, istat) return_status = istat END IF write_out_if END IF read_in_if CALL smooth_finalize(smooth, istat) endif num_days_if END SUBROUTINE smooth_ESI !=============================================================== ! Name: smooth_region ! ! Type: Subroutine ! ! Function: ! Smooth the ESI in the region. ! ! Description: ! Go through all the pixels and perform the smoothing for valid ! pixels. ! ! Files Needed: ! N/A ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! N/A ! ! Input: ! Name Type Description ! --------------------------------------------------- ! smooth I/O Data structure for SMOOTH ! ! 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 ! !=============================================================== SUBROUTINE smooth_region(smooth, return_status) !USES:: USE GETD_constants_mod USE GETD_log_mod USE file_io_mod !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- CHARACTER(LONG_STR) :: file_name INTEGER(Long) :: i, j, idoy,istrt,iend INTEGER(Long) :: input_unit, output_unit INTEGER(Long) :: record_length, option_flag INTEGER(Long) :: istat INTEGER(Long) :: fill_count, valid_count REAL(Single) :: cur_lon, cur_lat, cur_dist REAL(Single), ALLOCATABLE, DIMENSION(:) :: var_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varfilt1_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varfilt2_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varsm_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: var_fill ALLOCATE(var_array(smooth%in_day_num), stat = istat) ALLOCATE(varfilt1_array(smooth%in_day_num), stat = istat) ALLOCATE(varfilt2_array(smooth%in_day_num), stat = istat) ALLOCATE(varsm_array(smooth%in_day_num), stat = istat) ALLOCATE(var_fill(smooth%in_day_num), stat = istat) ! Main loop for the regional smoothing process ! xmask is not used allocate_if: IF (istat == SUCCESS ) THEN valid_count = 0 col_loop: DO i = 1, smooth%in_col_num row_loop: DO j = 1, smooth%in_row_num fill_count = 0 var_array = smooth%in_fill_value varfilt1_array = smooth%in_fill_value varfilt2_array = smooth%in_fill_value varsm_array = smooth%in_fill_value var_fill = smooth%in_fill_value ! Get the input time series variable for regional ! data matrix, count the number of valid values ! reset the time series array before performing ! the smooth process. in_doy_loop: DO idoy=1,smooth%in_day_num IF ( ABS(smooth%in_esi_array(i,j,idoy) - & smooth%in_fill_value) > 1e-6 ) THEN var_array(idoy)=smooth%in_esi_array(i,j,idoy) ELSE fill_count = fill_count + 1 ENDIF ENDDO in_doy_loop !call prefillseries(val,smooth%in_day_num) ! Start fill in the series if there are valid values ! Filter out outliers according to the data statistics fill_if: IF ( fill_count < smooth%in_day_num ) THEN CALL filter_series1(var_array,smooth%in_day_num, & smooth%in_fill_value, varfilt1_array, istat) CALL filter_series2(varfilt1_array,smooth%in_day_num, & smooth%in_fill_value,varfilt2_array,istat) CALL smooth_series(varfilt2_array,smooth%in_day_num, & smooth%in_fill_value, varsm_array,istat) ! if (ifill == 1) then ! call fill_series(valsm_array,valfill,istrt,iend,smooth%in_day_num) ! endif ! Put the smoothed results to the output matrix DO idoy = 1, smooth%in_day_num!istrt,iend IF (smooth%ifill == 1) THEN smooth%out_esi_array(i,j,idoy)=var_fill(idoy) ELSE smooth%out_esi_array(i,j,idoy)=varsm_array(idoy) ENDIF ENDDO valid_count = valid_count + 1 ENDIF fill_if END DO row_loop END DO col_loop ELSE WRITE(GETD_logunit,*) 'Failed to allocate the memory ! Program exit!' return_status = FAILURE ENDIF allocate_if ! Release the temp memory DEALLOCATE(var_fill, stat = istat) DEALLOCATE(var_array, stat = istat) DEALLOCATE(varfilt1_array, stat = istat) DEALLOCATE(varfilt2_array, stat = istat) DEALLOCATE(varsm_array, stat = istat) END SUBROUTINE smooth_region !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Debugging tool, load in a time series of data ! in ASCII file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE smooth_point(smooth, return_status) USE GETD_constants_mod USE GETD_log_mod USE file_io_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- CHARACTER(LONG_STR) :: file_name CHARACTER(LONG_STR) :: str_temp INTEGER(Long) :: i, j, idoy,istrt,iend INTEGER(Long) :: input_unit, output_unit INTEGER(Long) :: record_length, option_flag INTEGER(Long) :: istat INTEGER(Long) :: fill_count, file_status REAL(Single) :: cur_lon, cur_lat, cur_dist REAL(Single), ALLOCATABLE, DIMENSION(:) :: var_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varfilt1_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varfilt2_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: varsm_array REAL(Single), ALLOCATABLE, DIMENSION(:) :: var_fill REAL(Single), ALLOCATABLE, DIMENSION(:) :: alexi77_results ALLOCATE(var_array(smooth%in_day_num), stat = istat) ALLOCATE(varfilt1_array(smooth%in_day_num), stat = istat) ALLOCATE(varfilt2_array(smooth%in_day_num), stat = istat) ALLOCATE(varsm_array(smooth%in_day_num), stat = istat) ALLOCATE(var_fill(smooth%in_day_num), stat = istat) ALLOCATE(alexi77_results(smooth%in_day_num), stat = istat) var_array = smooth%in_fill_value var_fill = smooth%in_fill_value varfilt1_array = smooth%in_fill_value varfilt2_array = smooth%in_fill_value varsm_array = smooth%in_fill_value ! Main loop for the regional smoothing process ! xmask is not used allocate_if: IF (istat == SUCCESS ) THEN fill_count = 0 ! Get the input time series variable for regional ! data matrix, count the number of valid values input_unit = GETD_get_next_unit_number() OPEN(input_unit, & iostat = file_status, & file ='ESI_point.txt', & form = 'formatted', & status = 'old') IF (file_status /= SUCCESS) THEN WRITE(GETD_logunit,*)'Error in opening test point file' return_status = FAILURE RETURN ELSE ! Read in the point level data DO i = 1, 5 READ(input_unit,*) str_temp !Skip the commented line END DO DO idoy=1,smooth%in_day_num READ(input_unit,100) str_temp, i, var_array(idoy), & alexi77_results(idoy) ENDDO 100 FORMAT(a6, i4, 2f10.3) ENDIF ! Use read in from the file to replace ! the retrieved command line arguments ! Read in the information in the configuration file !c call prefillseries(val,smooth%in_day_num) ! Start fill in the series if there are valid values ! Filter out outliers according to the data statistics fill_if: IF ( fill_count < smooth%in_day_num ) THEN CALL filter_series1(var_array,smooth%in_day_num, & smooth%in_fill_value, varfilt1_array, istat) CALL filter_series2(varfilt1_array,smooth%in_day_num, & smooth%in_fill_value,varfilt2_array,istat) CALL smooth_series(varfilt2_array,smooth%in_day_num, & smooth%in_fill_value, varsm_array,istat) ! if (ifill == 1) then ! call fill_series(valsm_array,valfill,istrt,iend,smooth%in_day_num) ! endif ! Put the smoothed results to the output matrix WRITE(GETD_logunit,*)'Doy, IN ESI, ALEXI77, ALEXI90' DO idoy = 1, smooth%in_day_num!istrt,iend IF ( ABS(var_array(idoy) - smooth%in_fill_value) > 1e-6) THEN WRITE(GETD_logunit,'(i4,5(a1, f10.3))') idoy, ',',& var_array(idoy), ',', alexi77_results(idoy),',', & varfilt1_array(idoy), ',',& varfilt2_array(idoy), ',', varsm_array(idoy) ENDIF ENDDO ENDIF fill_if ELSE WRITE(GETD_logunit,*) 'Failed to allocate the memory ! Program exit!' return_status = FAILURE ENDIF allocate_if ! Release the temp memory DEALLOCATE(var_fill, stat = istat) DEALLOCATE(var_array, stat = istat) DEALLOCATE(varfilt1_array, stat = istat) DEALLOCATE(varfilt2_array, stat = istat) DEALLOCATE(varsm_array, stat = istat) CLOSE (input_unit) CALL GETD_release_unit_number(input_unit) END SUBROUTINE smooth_point !=============================================================== ! Name: read_smooth_config ! ! Type: Subroutine ! ! Function: ! Read in the spatial information of smooth EAST/WEST and ! meteorological domain. ! Description: ! N/A ! Files Needed: ! - in_file_name ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL read_smooth_config(in_file_name, smooth, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_file_name I Input configure file name ! ! Output: ! Name Type Description ! --------------------------------------------------- ! smooth I/O Data structure for SMOOTH ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! None ! ! History: ! 2014-09-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS ! !=============================================================== SUBROUTINE read_smooth_config(in_file_name, smooth, return_status) ! !USES:: USE GETD_constants_mod USE file_io_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- CHARACTER(LONG_STR), INTENT(in) :: in_file_name TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- INTEGER(Long) :: file_status, i, input_unit 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 ! Use read in from the file to replace ! the retrieved command line arguments ! Read in the information in the configuration file DO i = 1, 5 READ(input_unit,10) str_temp !Skip the commented line END DO ! Read in the input satellite information ! and the data spatial domain. ! Input smooth file ! READ(input_unit,10) smooth%in_esi_path READ(input_unit,10) str_temp smooth%in_esi_path = trim(env_var_value) // str_temp READ(input_unit,10) smooth%in_esi_file_name !READ(input_unit,10) smooth%in_snowmask_path READ(input_unit,10) str_temp smooth%in_snowmask_path = trim(env_var_value) // str_temp READ(input_unit,10) smooth%in_snowmask_file_name ! these info may be replaced by future domain file READ(input_unit,25) smooth%in_east_lon READ(input_unit,25) smooth%in_west_lon READ(input_unit,25) smooth%in_lon_resolution READ(input_unit,25) smooth%in_north_lat READ(input_unit,25) smooth%in_south_lat READ(input_unit,25) smooth%in_lat_resolution READ(input_unit,15) smooth%in_col_num READ(input_unit,15) smooth%in_row_num READ(input_unit,25) smooth%in_fill_value ! Output domain READ(input_unit,10) str_temp !comment line !READ(input_unit,10) smooth%out_esi_path READ(input_unit,10) str_temp smooth%out_esi_path = trim(env_var_value) // str_temp READ(input_unit,10) smooth%out_esi_file_name ! WRITE(GETD_logunit,*) 'smooth%in_col_num', smooth%in_col_num ! WRITE(GETD_logunit,*) 'smooth%in_row_num', smooth%in_row_num 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_smooth_config !=============================================================== ! Name: filter_series1 ! ! Type: Subroutine ! ! Function: ! Filter the input time series data (from filterseries:smooth_alexi.f). ! Iteratively screen all compressed data (BAD data removed) within a +/-IWIN window ! and flag pixels more that 2sigma above mean and 1.8sigma below mean (below avg ! values tend to be cloud-related, so more stringent test applied there.) ! ! Descriptions: ! First compress the input data by find the total ! number of valid data, then screen all the valid data. ! Files Needed: ! - N/A ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL filter_series1(in_vars, in_num, in_fillvalue, & ! out_vars, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_vars I Input time series data. ! in_num I Number of data in the in_vars ! in_fillvalue I Filled value, indicate the data is BAD ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_vars O Output filtered time series ! input time series data. ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! !=============================================================== SUBROUTINE filter_series1(in_vars, in_num, in_fillvalue, & out_vars, return_status) ! !USES:: USE GETD_constants_mod USE GETD_log_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- REAL(Single), DIMENSION(in_num), & INTENT(in) :: in_vars INTEGER(Long), INTENT( in) :: in_num REAL(Single), INTENT( in) :: in_fillvalue REAL(Single), DIMENSION(in_num), & INTENT(out) :: out_vars INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j, iter, j1, j2, jj INTEGER(Long) :: win_size, pre_valid_num, valid_num INTEGER(Long) :: x_num REAL :: x_sum, x_sum2, x_avg, x_sdev REAL(Single), DIMENSION(in_num) :: & temp_data, temp_dataf INTEGER, DIMENSION(in_num) :: data_index, dataf_index !Define the smooth window size win_size = 7 ! Compress data by filter out the filled value out_vars = in_fillvalue valid_num = 0 compress_loop: DO i = 1, in_num IF ( ABS(in_vars(i) - in_fillvalue) > 1e-6 ) THEN valid_num = valid_num + 1 temp_data(valid_num)=in_vars(i) data_index(valid_num)=i ENDIF END DO compress_loop ! Screen data for 6 times iter_loop: DO iter=1,6 pre_valid_num=valid_num filter_series_loop: DO j=1,valid_num ! find the valid data to compute the ! mean and variance to filter ! the outliers. ! For the values in the middle of the series, ! use the one window before ! and one window after the current value. ! But for the values in the beginning time ! series. ! Find the data in the two windows later. ! For the values in the end of series. IF (j < win_size) THEN j1=1 j2=win_size*2+1 j2=MIN(j2,valid_num) ELSE IF (j > valid_num-win_size) THEN j1=valid_num-(win_size*2+1) j1=MAX(j1,1) j2=valid_num ELSE j1=j-win_size j2=j+win_size ENDIF ! Now compute the statistics ! for the variables in the window x_num=0. x_sum=0. x_sum2=0 DO jj=j1,j2 x_num=x_num+1 x_sum=x_sum + temp_data(jj) x_sum2=x_sum2 + temp_data(jj)*temp_data(jj) ENDDO ! filter out the data if the ! value is out of 2sigma above mean and 1.8sigma below mean x_num_if: IF (x_num >= 3) THEN x_avg = x_sum/x_num x_sdev = SQRT((x_sum2/x_num)-x_avg*x_avg) IF ( (ABS(temp_data(j) - x_avg) >= 3.0*x_sdev) .OR. & (x_avg - temp_data(j) >= 2.8*x_sdev)) THEN temp_dataf(j)= in_fillvalue ELSE temp_dataf(j)=temp_data(j) ! Note, this is added by zp to track the index of the smoothed value dataf_index(j) = data_index(j) ENDIF ELSE temp_dataf(j)=temp_data(j) dataf_index(j) = data_index(j) ENDIF x_num_if ENDDO filter_series_loop !finished one series loop, ! Recompress, removing flagged points ! by recount the valid pixel numbers j=0 recompress_loop: DO i=1,valid_num IF ( ABS(temp_dataf(i) - in_fillvalue) > 1e-6) THEN j=j+1 temp_data(j)=temp_dataf(i) data_index(j)=dataf_index(i) ENDIF ENDDO recompress_loop valid_num=j ! Exit when we don't have any points ! out of range IF (valid_num == pre_valid_num) EXIT ENDDO iter_loop !Fill in the filtered data to the output matrix DO i=1,valid_num out_vars(data_index(i))=temp_data(i) ENDDO return_status = SUCCESS END SUBROUTINE filter_series1 !=============================================================== ! Name: filter_series2 ! ! Type: Subroutine ! ! Function: ! If there is a big gap at the beginning of the year, ! don't use data from last year. ! This helps to reduce smoothing across snow-induced gaps. ! ! Descriptions: ! First compress the input data by find the total ! number of valid data, then screen all the valid data. ! Files Needed: ! - N/A ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL filter_series2(in_vars, in_num, in_fillvalue, & ! out_vars, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_vars I Input time series data. ! in_num I Number of data in the in_vars ! in_fillvalue I Filled value, indicate the data is BAD ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_vars O Output filtered time series ! input time series data. ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! !=============================================================== SUBROUTINE filter_series2(in_vars, in_num, in_fillvalue, & out_vars, return_status) ! !USES:: USE GETD_constants_mod USE GETD_log_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- REAL(Single), DIMENSION(in_num), & INTENT(in) :: in_vars INTEGER(Long), INTENT( in) :: in_num REAL(Single), INTENT( in) :: in_fillvalue REAL(Single), DIMENSION(in_num), & INTENT(out) :: out_vars INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j, iter, j1, j2, jj INTEGER(Long) :: win_size, pre_valid_num, valid_num INTEGER(Long) :: x_num, ibuff, max_gap, gap_days ibuff = 50 max_gap=30 ! Copy the intput to the output list DO i=1,in_num out_vars(i)=in_vars(i) ENDDO ! go through the output series, filter ! out the previous year's data if the ! gap is too big in the input time series. i_loop: DO i = ibuff + 1, in_num IF ( ABS(in_vars(i) - in_fillvalue ) > 1e-6) THEN gap_days = i - ibuff IF (gap_days > max_gap)THEN ! The gap is too big DO j=1,ibuff out_vars(j)=in_fillvalue ENDDO ENDIF ! goto 100 EXIT i_loop ENDIF ENDDO i_loop !100 continue return_status = SUCCESS END SUBROUTINE filter_series2 !=============================================================== ! Name: fill_series ! ! Type: Subroutine ! ! Function: ! Fill in the time series data. ! ! Descriptions: ! Fill in big gaps using linear interpolation, otherwise spline can give big overshoots ! Spline interpolation of cloudy/missing pixels if grid filling is requested ! Flag times where big gaps exist at beginning or ending of timeseries ! Spline curve will head towards crazy values ! Files Needed: ! - N/A ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL fill_series(in_vars, in_start, in_end, in_num, & ! in_fillvalue, out_vars, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_vars I Input time series data. ! in_start I Start positions in_vars ! in_end I End position in_vars ! in_num I Number of data in the in_vars ! in_fillvalue I Filled value, indicate the data is BAD ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_vars O Output filtered time series ! input time series data. ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! Note !subroutine fillseries(val,valfill,istrt,iend,id) !=============================================================== SUBROUTINE fill_series(in_vars, in_start, in_end, & in_num, in_fillvalue, & out_vars, return_status) ! !USES:: USE GETD_constants_mod USE GETD_log_mod USE fft_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- REAL(Single), DIMENSION(in_num), & INTENT(in) :: in_vars INTEGER(Long), INTENT( in) :: in_start INTEGER(Long), INTENT( in) :: in_end INTEGER(Long), INTENT( in) :: in_num REAL(Single), INTENT( in) :: in_fillvalue REAL(Single), DIMENSION(in_num), & INTENT(out) :: out_vars INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, ii, n, iBAD INTEGER(Long) :: ni, nf, nlast, ngap INTEGER :: x_num, ibuff, max_gap REAL(Single), DIMENSION(in_num) :: temp_vars, x, y, y2 ! these are control the spline functions, ! since we only use natural sprline fit, ! these values are set to 1e30 REAL(Single) :: xi, yi, yp1, ypn !1. Move the input values into a temporary matrix ! and fill in big gaps (interval > 5 ) with linear method. ! Fill in big gaps using linear interpolation, otherwise spline can give big overshoots iBAD = -9999 ! This integer value is used to indicate the index value is empty nlast = iBAD fill_loop1: DO i = 1, in_num temp_vars(i) = in_vars(i) IF ( ABS(temp_vars(i) - in_fillvalue ) > 1e-6 ) THEN ! When there is a good value, check the gap ! between the current good value and the last good value n = i IF (nlast /= iBAD .AND. (n - nlast) > 5) THEN ! If the gap is large than 5 ! Linearly interpolate the values beween the gas ! for the 5 missing values. DO ii = nlast, n, 5 temp_vars(ii) = temp_vars(nlast)+ & float(ii-nlast)*(temp_vars(n)-temp_vars(nlast))/float(n-nlast) ENDDO ENDIF ! Record the position of the last good value as current. nlast = n ENDIF ENDDO fill_loop1 ! Spline interpolation of cloudy/missing pixels if grid filling is requested yp1=1e30 ! Natural spline fit ypn=1e30 n=0 ! fill_loop2: DO i=1,in_num IF (ABS(temp_vars(i) - in_fillvalue) > 1e-6) THEN n = n + 1 x(n)=i y(n)=temp_vars(i) ENDIF out_vars(i)= in_fillvalue ENDDO fill_loop2 ! When there is no valid values in the ! input time series, skip the steps !if (n.eq.0) goto 100 valid_value_exists_if: IF ( n > 0 ) THEN n_one_if: IF ( n == 1 ) THEN !When we only have one single value ! use this value to fill in all the DO i = 1, in_num out_vars(i)=temp_vars(INT(x(n))) END DO ! if (n.eq.1) then ! i=x(n) ! out_vars(i)=temp_vars(i) ! goto 100 ! endif ELSE ! More than one valid values in the ! time series and no gap is larger ! than 5 (or already been filled). ! Use Spline fit to fill in the values CALL spline(x,y,in_num,yp1,ypn,y2) DO i=in_start,in_end xi=i CALL splint(x,y,y2,in_num,xi,yi) out_vars(i)=yi ENDDO !c Flag times where big gaps exist at beginning or ending of timeseries !c Spline curve will head towards crazy values ! zp comments: This may not happen since we have already fill in these ! large gaps ni=iBAD nf=iBAD ngap=5 ! go through the input time series and save the position ! of the 1st valid value in ni, ! last valid value in nf DO i=in_start,in_end IF ((ABS(temp_vars(i) - in_fillvalue) > 1e-6 ) & .AND. ni == iBAD) ni=i IF (ABS(temp_vars(i) - in_fillvalue) > 1e-6 ) nf=i ENDDO ! If either no start or no end value, IF (ni == iBAD .OR. nf == iBAD) THEN DO i=in_start,in_end out_vars(i)=in_fillvalue ENDDO ELSE IF (ni-in_start.GT.ngap) THEN DO i=in_start,ni-1 out_vars(i) = in_fillvalue ENDDO ENDIF IF ( (in_end - nf) > ngap) THEN DO i=nf-1,in_end out_vars(i)= in_fillvalue ENDDO ENDIF ENDIF ENDIF n_one_if END IF valid_value_exists_if ! fill process done !100 continue return_status = SUCCESS END SUBROUTINE fill_series !=============================================================== ! Name: smooth_series ! ! Type: Subroutine ! ! Function: ! This is the main subroutine to smooth the time series ! variables. ! ! Descriptions: ! First send compressed time-series through Savitsky-Golay ! Filter, then apply derived convolution factors to get smoothed ! series. ! ! Files Needed: ! - N/A ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL smooth_series(in_vars, in_num, in_fillvalue, & ! out_vars, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_vars I Input time series data. ! in_num I Number of data in the in_vars ! in_fillvalue I Filled value, indicate the data is BAD ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_vars O Output smoothed time series ! data. ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! !=============================================================== SUBROUTINE smooth_series(in_vars, in_num, in_fillvalue, & out_vars, return_status) ! !USES:: USE GETD_constants_mod USE GETD_log_mod USE fft_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- REAL(Single), DIMENSION(in_num), & INTENT(in) :: in_vars INTEGER(Long), INTENT( in) :: in_num REAL(Single), INTENT( in) :: in_fillvalue REAL(Single), DIMENSION(in_num), & INTENT(out) :: out_vars INTEGER(Long), INTENT(out) :: return_status ! Local variables ! These are the parameters used to call the convolution ! function ! NSIZE is the maximum number for the smoothing matrix !The actual size is npt, this shoudl be less than NSIZE INTEGER(Long), PARAMETER :: NSIZE = 512 INTEGER(Long), PARAMETER :: NP = 11 INTEGER(Long), PARAMETER :: NL = 5 INTEGER(Long), PARAMETER :: NR = 5 INTEGER(Long), PARAMETER :: MPARM = 2 INTEGER(Long) :: ld, npt, in_size INTEGER(Long) :: i, j, k REAL(DOUBLE), DIMENSION(NSIZE) :: cur_data REAL(DOUBLE), DIMENSION(NSIZE) :: cur_ans REAL(DOUBLE), DIMENSION(NSIZE) :: cur_c REAL(DOUBLE), DIMENSION(NSIZE) :: cur_respns REAL(DOUBLE) :: cur_sum, pad_val1, pad_val2, xin ! Send compressed time-series through Savitsky-Golay Filter, then apply derived ! convolution factors to get smoothed series. ! ld is the flag used to ld = 0 DO k=1,NSIZE cur_data(k)=0. cur_ans(k)=0. ENDDO ! go through the input time series, Compress data and pad out front ! and back end to accommodate window +nr/-nl pixels wide ! Pad values represent averages of 5 first/last values. ! If there are less than 5 good values, average over all values j = 0 i_loop: DO i = 1, in_num IF ( ABS(in_vars(i) - in_fillvalue ) > 1e-6) THEN j = j + 1 cur_data(j+NL) = in_vars(i) ENDIF ENDDO i_loop npt = j in_size = MIN(npt, 5) cur_sum = 0 !Right data +nr DO k = 1, in_size cur_sum = cur_sum + cur_data(k + NL) END DO pad_val1 = cur_sum / REAL(in_size) !Left data /-nl cur_sum = 0 DO k = 1, in_size cur_sum = cur_sum + cur_data(npt - in_size + k + NL) END DO pad_val2 = cur_sum / REAL(in_size) DO k = 1, NL cur_data(k) = pad_val1 END DO DO k = npt + 1, npt + NR cur_data(k+NL) = pad_val2 END DO ! Call SG filter, then apply convolution, then decompress CALL savgol(cur_c, NP, NL, NR, ld, MPARM) DO i = 1,NP cur_respns(i)=cur_c(i) ENDDO CALL convlv(cur_data,NSIZE,cur_respns,NP,1,cur_ans) j=0 out_vars = in_fillvalue out_loop: DO i = 1, in_num IF ( ABS(in_vars(i) - in_fillvalue ) > 1e-6) THEN j = j + 1 out_vars(i) = cur_ans(j + NL) ENDIF END DO out_loop return_status = SUCCESS END SUBROUTINE smooth_series !=============================================================== ! Name: read_in_files ! ! Type: Subroutine ! ! Function: ! This subroutine is to read in the ESI and snow mask data ! in multiple days for the smoothing process. ! ! Description: ! Instead of check for bad values, only put in the valid ! value when there is no snow, and fraction is valid ! Files Needed: ! - esi_series.bin ! - snow_mask_yyyyddd.bin ! ! Modules Needed: ! - TYPE_kinds ! - GETD_constants_mod ! - GETD_log_mod ! - file_io_mod ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL read_in_files(smooth, return_status) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! smooth I/O Number of channels ! ! Output: ! None ! ! History: ! 2014-09-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS !note: the input has to be in the correct format ! 2014-03-25 Major update with change in the file name ! and the location. ! !=============================================================== SUBROUTINE read_in_files(smooth, return_status) ! ---------- ! Module use ! ---------- USE type_kinds USE GETD_constants_mod USE GETD_log_mod USE file_io_mod ! --------------------------- ! Disable all implicit typing ! --------------------------- IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! Local variables CHARACTER(LONG_STR) :: file_name !matrix holding daily information REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: snow_mask REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: daily_esi INTEGER(Long) :: istat, i, j, m, n INTEGER(Long) :: icount, day_num INTEGER(Long) :: cur_year, cur_doy INTEGER(Long) :: start_doy, end_doy return_status = SUCCESS ! Allocate memory for two temporary matrix ALLOCATE(snow_mask(smooth%in_col_num, & smooth%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate 2D snow mask error' return_status = FAILURE RETURN END IF ALLOCATE(daily_esi(smooth%in_col_num, & smooth%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate 2D input ESI error' return_status = FAILURE RETURN END IF !Main loop to read in all the input time series !Note, this requires the correct location for the input data icount = 0 day_num = 0 ! Load in the ESI and snow mask data information year_loop: DO i = smooth%start_year, smooth%end_year ! ONly for the start year, start from input ! otherwise it should start from the 1st day of the year if ( i > smooth%start_year ) then start_doy = 1 else start_doy = smooth%start_doy endif ! Only the last year used the end of day if ( i < smooth%end_year ) then ! could be a leap year IF ( ((MOD(i, 4) == 0) .AND. (MOD(i, 100) /= 0)) .OR. & (MOD(i, 400) == 0)) THEN end_doy = 366 ELSE end_doy = 365 ENDIF else end_doy = smooth%end_doy endif day_loop: DO j = start_doy, end_doy ! Read in daily snow mask into the time series ! Note this need to be the same as input file name and path. WRITE(file_name,'(a, a1, i4.4, i3.3,a1, a)') & TRIM(smooth%in_esi_path), '/', i, j, & '/', TRIM(smooth%in_esi_file_name) ! e.g. ! /data/data123/zpli/data/VLAI_ESI/2014224/ CALL read_2d_real_by_col_row(file_name, & smooth%in_col_num, smooth%in_row_num, smooth%in_fill_value,& daily_esi, istat) ! This value just record the number of days, not ! judge if the read is successful or not day_num = day_num + 1 ! Save the data into the input matrix if read successful ! Otherwise skip it IF ( istat == SUCCESS ) THEN DO m = 1, smooth%in_col_num DO n = 1, smooth%in_row_num smooth%in_esi_array(m,n,day_num) = daily_esi(m,n) END DO END DO icount = icount + 1 ELSE write(GETD_logunit,*)'Error reading ', file_name ENDIF !Snow Mask file ! WRITE(file_name,'(a, a1, i4.4, a1, a, i4.4, i3.3, a4)') & ! TRIM(smooth%in_snowmask_path), '/', i, & ! '/', TRIM(smooth%in_snowmask_file_name), i, j, '.dat' WRITE(file_name,'(a, a1, i4.4, i3.3,a1, a)') & TRIM(smooth%in_snowmask_path), '/', i, j, & '/', TRIM(smooth%in_snowmask_file_name) CALL read_2d_real_by_col_row(file_name, & smooth%in_col_num, smooth%in_row_num, smooth%in_fill_value,& snow_mask, istat) IF ( istat == SUCCESS ) THEN DO m = 1, smooth%in_col_num DO n = 1, smooth%in_row_num smooth%in_snowmask_array(m,n,day_num) = snow_mask(m,n) END DO END DO ELSE write(GETD_logunit,*)'Error reading ', file_name ENDIF END DO day_loop END DO year_loop WRITE(GETD_logunit,*)'Total ', day_num, ' days with ', icount, & ' ESI read successfully.' DEALLOCATE(snow_mask, stat = istat) DEALLOCATE(daily_esi, stat = istat) END SUBROUTINE read_in_files !=============================================================== ! Name: smooth_init ! ! Type: Subroutine ! ! Function: ! Allocate the memory for SNOW module. ! Description: ! ! Files Needed: ! - None ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! N/A ! ! Input: ! Name Type Description ! --------------------------------------------------- ! -smooth Derived ! ! Output: ! Name Type Description ! --------------------------------------------------- ! smooth D Allocatable memory ! return_status I 0 - Success, 1 - Failure ! ! System Calls: ! None ! ! History: ! 2014-09-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS !=============================================================== SUBROUTINE smooth_init(smooth, return_status) USE GETD_constants_mod IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! --------------- ! Local variables ! --------------- INTEGER(Long) :: istat return_status = SUCCESS ! ------------------------------------------------------------------- ! Allocate the memory for the 2-D smooth data boundary_if: IF (smooth%in_col_num < 1 .OR. & smooth%in_row_num < 1 .OR. smooth%in_day_num < 1) THEN WRITE(GETD_logunit,*)'Input boundary error ! ESI init process exit.' return_status = FAILURE ELSE ! Input data, smooth input mask. ALLOCATE(smooth%in_esi_array(smooth%in_col_num, & smooth%in_row_num, smooth%in_day_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate input ESI error' return_status = FAILURE ELSE smooth%in_esi_array = smooth%in_fill_value END IF !snow mask ALLOCATE(smooth%in_snowmask_array(smooth%in_col_num, & smooth%in_row_num, smooth%in_day_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate input snow mask error' return_status = FAILURE ELSE smooth%in_snowmask_array = smooth%in_fill_value END IF ! Output data ALLOCATE(smooth%out_esi_array(smooth%in_col_num, & smooth%in_row_num, smooth%in_day_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*) 'col_num:', smooth%in_col_num WRITE(GETD_logunit,*) 'row_num:', smooth%in_row_num WRITE(GETD_logunit,*) 'day_num:', smooth%in_day_num WRITE(GETD_logunit,*)'Allocate out ESI array error' return_status = FAILURE ELSE smooth%out_esi_array = smooth%in_fill_value END IF END IF boundary_if END SUBROUTINE smooth_init !=============================================================== ! Name: smooth_finalize ! ! Type: Subroutine ! ! Function: ! Release the memory ! Description: ! ! Files Needed: ! - None ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! N/A ! ! Input: ! Name Type Description ! --------------------------------------------------- ! -smooth I/O Data type holding the data ! ! Output: ! Name Type Description ! --------------------------------------------------- ! return_status I 0 - Success, 1 - Failure ! ! System Calls: ! None ! ! History: ! 2014-09-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS !=============================================================== SUBROUTINE smooth_finalize(smooth, return_status) USE GETD_constants_mod IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! --------------- ! Local variables ! --------------- INTEGER(Long) :: istat IF ( ALLOCATED(smooth%in_esi_array) ) THEN DEALLOCATE(smooth%in_esi_array, stat = istat) ENDIF IF ( ALLOCATED(smooth%in_snowmask_array) ) THEN DEALLOCATE(smooth%in_snowmask_array, stat = istat) ENDIF IF ( ALLOCATED(smooth%out_esi_array) ) THEN DEALLOCATE(smooth%out_esi_array, stat = istat) ENDIF return_status = SUCCESS END SUBROUTINE smooth_finalize !=============================================================== ! Name: write_out_files ! ! Type: Subroutine ! ! Function: ! This subroutine is to write the smoothed ESI data ! into multiple daily files. ! ! Description: ! Write out the smoothed ESI into destination files. ! ! Files Needed: ! ! Modules Needed: ! - file_io_mod ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL write_out_files(smooth, return_status) ! ! Input: ! Name Type Description ! -------------------------------------------------------- ! smooth I ! ! Output: ! Name Type Description ! -------------------------------------------------------- ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! History: ! 2015-01-15 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS ! !=============================================================== SUBROUTINE write_out_files(smooth, return_status) ! !USES:: USE GETD_constants_mod USE GETD_log_mod USE file_io_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- TYPE(smooth_ESI_type) :: smooth INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- CHARACTER(LONG_STR) :: path_name CHARACTER(LONG_STR) :: file_name INTEGER(Long) :: i, j, icount, istat INTEGER(Long) :: start_doy, end_doy ! Write smoothed results to multiple files icount = 1 year_loop: DO i = smooth%start_year, smooth%end_year ! ONly for the start year, start from input ! otherwise it should start from the 1st day of the year if ( i > smooth%start_year ) then start_doy = 1 else start_doy = smooth%start_doy endif ! Only the last year used the end of day if ( i < smooth%end_year ) then ! could be a leap year IF ( ((MOD(i, 4) == 0) .AND. (MOD(i, 100) /= 0)) .OR. & (MOD(i, 400) == 0)) THEN end_doy = 366 ELSE end_doy = 365 ENDIF else end_doy = smooth%end_doy endif day_loop: DO j = start_doy, end_doy ! Create the directory if not exist WRITE(path_name,'(a, a1, i4.4, i3.3, a1)') & TRIM(smooth%out_esi_path), '/', i, j,'/' CALL SYSTEM('mkdir -p ' // TRIM(path_name)) WRITE(file_name,'(a, a1, i4.4, i3.3, a1, a)') & TRIM(smooth%out_esi_path), '/', i, j, & '/', TRIM(smooth%out_esi_file_name) if ( icount <= smooth%in_day_num ) then write(*,*)'Write to ', trim(file_name) CALL write_out_2d_real( smooth%out_esi_array(:,:,icount),& smooth%in_col_num, smooth%in_row_num, file_name, istat) icount = icount + 1 else write(GETD_logunit,*)'Output days is larger than data content' endif END DO day_loop END DO year_loop END SUBROUTINE write_out_files END MODULE ESI_smooth_mod !EOP