!-------------------------------------------------------------------- !Name: ESI_data_mod ! !Type: F90 module ! !Function: ! This module compute the ESI from ALEXI ET and other input files. ! !Description: ! There are four components to compute the ESI from ET data ! 1. Compute the daily reference ET using climatology input solar radiation. ! 2. Compute the daily ESI by dividing the ET with reference ET. ! 3. Smoothe the ESI for multiple days. ! 4. Compute the composite ESI 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: ! ESI.cfg - configuration file ! insolar.bin - climatology input radiation ! temp.bin - Daily surface air temperature ! wind.bin - Daily surfce wind ! q2.bin - climatology input radiation! ! !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-04-28 Zhengpeng Li (UMD ESSIC) ! ! !-------------------------------------------------------------------- MODULE ESI_data_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 :: Compute_REFET_ESI PRIVATE :: read_ESI_config PRIVATE :: compute_REFET_data PRIVATE :: ESI_fao_PM TYPE ESI_data_type CHARACTER(LONG_STR) :: in_path CHARACTER(LONG_STR) :: in_et_file_name CHARACTER(LONG_STR) :: in_rnet2_file_name CHARACTER(LONG_STR) :: in_sdn_file_name CHARACTER(LONG_STR) :: in_ta_file_name CHARACTER(LONG_STR) :: in_ea_file_name CHARACTER(LONG_STR) :: in_pres_file_name CHARACTER(LONG_STR) :: in_wind_file_name ! Output reference ET and ESI files CHARACTER(LONG_STR) :: out_path CHARACTER(LONG_STR) :: out_esi_file_name CHARACTER(LONG_STR) :: out_ref_et_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 ! corrected temperature at time 1 and 2 ! REAL(Single), allocatable, dimension(:,:) :: in_ET_data ! ! REAL(Single), allocatable, dimension(:,:) :: out_fPET_data ! REAL(Single), allocatable, dimension(:,:) :: out_col_array ! REAL(Single), allocatable, dimension(:,:) :: out_row_array END TYPE ESI_data_type CONTAINS !=============================================================== ! Name: Compute_REFET_ESI ! ! Type: Subroutine ! ! Function: ! Compute the daily ESI value from input ET and calculated ! reference ET in the domain. ! ! Descriptions: ! Read in the Meteorological input files to compute the ! reference ET. Then read in the daily ET value, to compute the ! daily ESI. ! ! Files Needed: ! - configure file ! - Input ET data file ! - Input Meterorlogical data files ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! CALL Compute_REFET_ESI( 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 ! The conversion between ALEXI outputs ! with ET data is: ! Equivalent radiation in megajoules per square ! metre per day (MJ m-2 day-1) ! 1 MJ m-2 day-1 = 0.408 mm day-1 !=============================================================== SUBROUTINE Compute_REFET_ESI( in_config_file, & in_year, in_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) :: in_year, in_doy INTEGER(Long), INTENT(out) :: return_status ! Local variables INTEGER(Long) :: i, j INTEGER(Long) :: col, row INTEGER(Long) :: istat, icount INTEGER :: cur_col, cur_row CHARACTER(LONG_STR) :: file_name TYPE(ESI_data_type) :: ESI REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_ET_data REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: in_RNET2_data REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ref_ET_data REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: out_fPET_data ! Read in the configuration return_status = SUCCESS CALL read_ESI_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 ALLOCATE(ref_ET_data(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for reference ET error' return_status = FAILURE RETURN END IF ! Compute the Reference ET DATA CALL compute_REFET_data(ESI, ref_ET_data, istat) refet_if: IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit, *)'Error in computing the Reference ET,', & 'Program exit error in', 'Compute_REFET_ESI' return_status = FAILURE ELSE ! Load the input ET data ALLOCATE(in_ET_data(ESI%in_col_num, ESI%in_row_num), stat = istat) WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_et_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_ET_data, istat) WRITE(GETD_logunit, *)'Input ET file,', TRIM(file_name) load_ET_if: IF ( istat == SUCCESS ) THEN ! Load the input RNET2 data ALLOCATE(in_RNET2_data(ESI%in_col_num, ESI%in_row_num), stat = istat) WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_rnet2_file_name) WRITE(GETD_logunit, *)'Input RNET2 file,', TRIM(file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, in_RNET2_data, istat) ! Compute the ESI data ALLOCATE(out_fPET_data(ESI%in_col_num, ESI%in_row_num), stat = istat) IF ( istat /= SUCCESS ) THEN WRITE(GETD_logunit,*)'Allocate memory for ESI data error' return_status = FAILURE RETURN END IF WRITE(GETD_logunit,*) 'Compute ESI using ET and reference ET:' out_fPET_data = ESI%in_fill_value out_col_loop: DO i = 1, ESI%in_col_num out_row_loop: DO j = 1, ESI%in_row_num if(in_RNET2_data(i,j) > 0 .AND. & ref_ET_data(i,j)/in_RNET2_data(i,j) < 0.5) then in_ET_data(i,j) = ESI%in_fill_value endif IF ( in_ET_data(i,j) > 0 .AND. & ref_ET_data(i,j) > 0 ) THEN out_fPET_data(i,j) = in_ET_data(i,j) / ref_ET_data(i,j) if ( out_fPET_data(i,j) > 5.0 ) then if ( GETD_debug > 0 ) then write(GETD_logunit,*) 'ESI outlier:', i, j, & in_ET_data(i,j), ref_ET_data(i,j) ENDIF endif ENDIF END DO out_row_loop END DO out_col_loop ! write out the results to the output files ! if the computation is successful write_out_if: IF ( return_status == SUCCESS ) THEN ! Write reference ET WRITE(file_name, '(a,a1,a)') TRIM(ESI%out_path), & '/', TRIM(ESI%out_ref_et_file_name) CALL write_out_2d_real(ref_ET_data, ESI%in_col_num, & ESI%in_row_num, file_name, istat) ! Write ESI WRITE(file_name, '(a,a1,a)') TRIM(ESI%out_path), & '/', TRIM(ESI%out_esi_file_name) CALL write_out_2d_real(out_fPET_data, ESI%in_col_num, & ESI%in_row_num, file_name, istat) return_status = SUCCESS END IF write_out_if ELSE return_status = FAILURE ENDIF load_ET_if END IF refet_if IF ( ALLOCATED(in_ET_data) ) DEALLOCATE(in_ET_data) IF ( ALLOCATED(ref_ET_data)) DEALLOCATE(ref_ET_data) IF ( ALLOCATED(out_fPET_data)) DEALLOCATE(out_fPET_data) END SUBROUTINE Compute_REFET_ESI !=============================================================== ! Name: read_ESI_config ! ! Type: Subroutine ! ! Function: ! Read in the spatial information of ESI EAST/WEST and ! meteorological domain. ! Description: ! N/A ! Files Needed: ! - in_file_name ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! N/A ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_file_name I configure file name ! ! Output: ! Name Type Description ! --------------------------------------------------- ! ESI O Data type holding configure ! information and dynamic memory. ! 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_ESI_config(in_file_name, ESI, return_status) ! ---------- ! Module use ! ---------- USE type_kinds USE GETD_constants_mod USE GETD_log_mod USE file_io_mod ! !PUBLIC TYPES: IMPLICIT NONE ! --------- ! Arguments ! --------- CHARACTER(LONG_STR), INTENT(in) :: in_file_name TYPE(ESI_data_type) :: ESI INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- INTEGER(Long) :: file_status, i, input_unit CHARACTER(LONG_STR) str_temp input_unit = GETD_get_next_unit_number() 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) str_temp !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_path READ(input_unit,10) ESI%in_et_file_name READ(input_unit,10) ESI%in_rnet2_file_name !Climatology insolation READ(input_unit,10) ESI%in_sdn_file_name !Daily meteorological data READ(input_unit,10) ESI%in_ta_file_name READ(input_unit,10) ESI%in_ea_file_name READ(input_unit,10) ESI%in_pres_file_name READ(input_unit,10) ESI%in_wind_file_name ! Output reference ET and ESI files READ(input_unit,10) str_temp !comment line READ(input_unit,10) ESI%out_path READ(input_unit,10) ESI%out_ref_et_file_name READ(input_unit,10) ESI%out_esi_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_ESI_config !=============================================================== ! Name: compute_REFET_data ! ! Type: Subroutine ! ! Function: ! Compute the reference ET from input meteorological data. ! Description: ! Use P-M formula to compute the ! Files Needed: ! - in_data_file ! ! Modules Needed: ! - None ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! call compute_REFET_data(file_name, ref_ET_data, istat) ! ! Input: ! Name Type Description ! --------------------------------------------------- ! ESI O Data type holding configure ! information and dynamic memory. ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_et_array I Data matrix holding computed ET ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! None ! ! History: ! 2015-01-06 Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS !=============================================================== SUBROUTINE compute_REFET_data(ESI, out_et_array, return_status) ! !USES:: ! ---------- ! 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(ESI_data_type), INTENT(in) :: ESI REAL(Single), DIMENSION(:,:), INTENT(out) :: out_et_array INTEGER(Long), INTENT(out) :: return_status ! ---------------- ! Local variables ! --------------- INTEGER(Long) :: file_status, i, j INTEGER(Long) :: istat INTEGER(Long) :: icount, ierror! statistics CHARACTER(LONG_STR) file_name REAL :: temp_et REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: sdn_array REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ta_array REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: ea_array REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: pres_array REAL(Single), ALLOCATABLE, DIMENSION(:,:) :: wind_array ALLOCATE (sdn_array(ESI%in_col_num, ESI%in_row_num), stat = istat) ALLOCATE (ta_array(ESI%in_col_num, ESI%in_row_num), stat = istat) ALLOCATE (ea_array(ESI%in_col_num, ESI%in_row_num), stat = istat) ALLOCATE (pres_array(ESI%in_col_num, ESI%in_row_num), stat = istat) ALLOCATE (wind_array(ESI%in_col_num, ESI%in_row_num), stat = istat) ! Fill in the output with filled values out_et_array = ESI%in_fill_value return_status = SUCCESS ! Read in all the input meteorological variables ! Insolation input (from climatology WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_sdn_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, sdn_array, istat) ! Air temperature WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_ta_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, ta_array, istat) if ( istat /= SUCCESS ) then return_status = FAILURE endif ! SPFH WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_ea_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, ea_array, istat) ! Surface pressure WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_pres_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, pres_array, istat) ! Surface wind WRITE(file_name,'(a,a1,a)')TRIM(ESI%in_path), '/', & TRIM(ESI%in_wind_file_name) CALL read_2d_real_by_col_row(file_name, & ESI%in_col_num, ESI%in_row_num, & ESI%in_fill_value, wind_array, istat) IF (istat /= SUCCESS) THEN WRITE(GETD_logunit,*)'Error in opening the input data files' return_status = FAILURE ELSE ! Read in the ESI mask data file ! Skip the comment lines icount = 0 ierror = 0 out_col_loop: DO i = 1, ESI%in_col_num out_row_loop: DO j = 1, ESI%in_row_num IF ( ABS(pres_array(i,j) - ESI%in_fill_value) > 1e-6 & .AND. ABS(ta_array(i,j) - ESI%in_fill_value) > 1e-6 & .AND. ABS(ea_array(i,j) - ESI%in_fill_value) > 1e-6 & .AND. ABS(wind_array(i,j) - ESI%in_fill_value) > 1e-6 & .AND. ABS(sdn_array(i,j) - ESI%in_fill_value) > 1e-6 & ) THEN ! Note, the input should be Air temperature in Kelvin ! Pa and q, the conversion will be done inside the ! FAO function. CALL ESI_fao_PM(sdn_array(i,j), ta_array(i,j), & pres_array(i,j),ea_array(i,j),wind_array(i,j), & out_et_array(i,j), istat) IF ( istat == SUCCESS ) THEN icount = icount + 1 ELSE ierror = ierror + 1 ENDIF ENDIF END DO out_row_loop END DO out_col_loop ENDIF WRITE(GETD_logunit,*) 'Valid', icount, 'Error', ierror DEALLOCATE (sdn_array, stat = istat) DEALLOCATE (ta_array, stat = istat) DEALLOCATE (ea_array, stat = istat) DEALLOCATE (pres_array, stat = istat) DEALLOCATE (wind_array, stat = istat) return_status = SUCCESS END SUBROUTINE compute_REFET_data !=============================================================== ! Name: ESI_fao_PM ! ! Type: Subroutine ! ! Function: ! Compute the reference ET with climatology insolation, used ! for ESI computation only. No time of day and zenith inputs, ! Comparing with ALEXI_fao_PM, this one used constant partition ! for the input radiation. ! ! Description: ! Computes grass reference ET based primarily on FAO Penman-Monteith ! (Equation 53 for hourly ETo) with the assumption that ! albedo is 0.23 constant. ! ! Files Needed: ! - None ! ! Modules Needed: ! - GETD_constant_mod ! ! Subroutines Contained: ! N/A ! ! Calling Sequence: ! ! Input: ! Name Type Description ! --------------------------------------------------- ! in_sdn I Solar Insolation ! in_ta I Air temperature in Celcius ! in_pres I Air pressure ( 0.01 bar) ! in_ea I Specific Humidity (1000) ! in_wind I Wind speed ! ! Output: ! Name Type Description ! --------------------------------------------------- ! out_eref O Reference ET computed ! return_status O Flag, 0 - SUCCESS, 1 - FAILURE ! ! System Calls: ! None ! ! History: ! 08/31/2011 Created by Martha Anderson ! 04/04/2014 Modified by Zhengpeng Li (UMD ESSIC) NOAA/NESDIS/STAR/CICS ! Add a check for the input data, if not valid return with filled value. ! !=============================================================== SUBROUTINE ESI_fao_PM(in_sdn, & in_ta,in_pres,in_ea,in_wind, & out_eref, return_status) ! ---------- ! Module use ! ---------- USE type_kinds USE GETD_constants_mod USE GETD_log_mod ! --------------------------- ! Disable all implicit typing ! --------------------------- IMPLICIT NONE ! ARGUMENTS: REAL, INTENT(in) :: in_sdn REAL, INTENT(in) :: in_ta ! air temperature in Celcius REAL, INTENT(in) :: in_pres REAL, INTENT(in) :: in_ea REAL, INTENT(in) :: in_wind REAL, INTENT(out) :: out_eref INTEGER, INTENT(out):: return_status ! Local variables: REAL :: albedo, xlam, psych, esat, s, d REAL :: rnet, temp_sdn REAL :: rnetmj, gmj, xnum1, xnum2, xnum, xden REAL :: erefmmhr REAL :: tak, esky_clear, esky, esfc, fclear REAL :: rnet_canopy, rnet_soil ! variables used to convert the units REAL :: cur_ta, cur_ea, cur_pres INTEGER :: istate ! ALEXI_SPECIFIC_AIR_HEAT REAL, PARAMETER :: cp = 1010.e-6 ! [MJ/kg-K], REAL, PARAMETER :: f_epsilon = 0.622 REAL, PARAMETER :: sigma = 5.67e-8 ! [W/(m2-K4)] ! Convert the input variable units cur_ea = in_ea * 1000 cur_ta = in_ta - 273.15 cur_pres = in_pres * 0.01 albedo = 0.23 ! [ Hypothetical grass reference] fclear = 1.0 ! clear sky fraction, xlam=(2.501-0.00237*cur_ta) ! latent heat o'vaporization [MJ kg-1] ! psych=0.1*pres*pa/(xlam*f_epsilon) ! psychometric constant [kPa C-1] ! Note that ALEXI_SPECIFIC_AIR_HEAT is in [J/kg-K] and this formula needs ! [MJ/kg-K] ! Note, the input Air temperature should be ! These are different than ALEXI_fao_PM tak=cur_ta + 273.15 psych=0.1*cur_pres*cp /(xlam * f_epsilon) esat=0.6108*EXP(17.2694*cur_ta/(237.3+cur_ta)) ! [kPa] s=17.2694*237.3*esat/(237.3+cur_ta)**2 ! [kPa C-1] d = esat - 0.1 * cur_ea ! [kPa] ! Clear: Swinbank '63 (C&N pg 163) esky_clear=9.2e-6*tak*tak ! Cloudy: Monteith & Unsworth '90 (C&N pg 164) esky=(1-0.84*fclear)*esky_clear + 0.84*fclear esfc=0.98 rnet_canopy=esfc*(esky-1.0)*sigma*(tak**4) rnet_canopy=esfc*esky*sigma*(tak**4)-esfc*sigma*((tak+4)**4) rnet_soil=in_sdn*(1-albedo) rnet = rnet_canopy + rnet_soil ! Convert the unit, note that rnet cannot be negative IF ( rnet > 1e-6 ) THEN rnetmj=rnet*3600./1.e6 ! [MJ/m2-hr] gmj=0.1*rnetmj xnum1=0.408*s*(rnetmj-gmj) ELSE ! Added this to prevent the error when input net radiation ! is 0 rnetmj = 0 gmj = 0 xnum1 = 0 ENDIF xnum2=psych*37.*in_wind*d/tak xnum=xnum1+xnum2 xden=s+psych*(1.+0.24*in_wind) erefmmhr=xnum/xden ! [mm/hr] out_eref = erefmmhr * xlam * 1.e6/3600. ! [W/m2] IF ( ABS(out_eref) > 1000 ) THEN WRITE(GETD_logunit,*) 'erefmmhr', erefmmhr, 'xlam', xlam, 'xden', xden WRITE(GETD_logunit,*) 'Error in eref', out_eref WRITE(GETD_logunit,*) 'SDN, TA inputs', in_sdn,in_ta WRITE(GETD_logunit,*) 'PRES, EA, WIND inputs', in_pres, in_ea, in_wind ENDIF return_status = SUCCESS END SUBROUTINE ESI_fao_PM !--------------------------------------------------------------------- END MODULE ESI_data_mod !EOP