module ncoda_climo2_module ! - This module contains subroutines for calculating interplated data from the ! updated 2005 - 2017 NCODA and DSST climatology for a given lat, lon, and time ! - written by G. Chirokova (CIRA/CSU) Dec 28, 2017, ! parts of the code adapted from the original ncodaclimo.f ! and John Knaff's (STAR/NESDIS) climatology code ! - v1.1 - Dec 29, 2017 - added option to select needed climo variables ! from the calling subroutine (G. Chirokova) ! - v1.2 - Jan 1, 2018 - added more logging messages ! and fixed minor bugs with error handling. ALso, added ensuring that on ! any error all returned climo values are set to missing (G. Chirokova) ! - v1.3 - Jan 2, 2018 - fixed bugs in leap_year and length of cvar. Also ! changed lu_climo to 2l. Tested th ecode with pgf90 compiler (G Chirokova) ! - v2.1 - Aug 10, 2018 - modified to process both 0.25 deg and 0.5 deg ! climatology. (G Chirokova) ! - this verion works with 0.25 deg climatology ver 2017_v2 ! and with 0.5 deg climatology ver 2018_v2.1 ! - v2.2 - Aug 13, 2018 - modified by G Chirokova and J. Knaff ! - cleaned up the code to remove unnesessary duplicate arrays ! - v2.2.1 - Sep 16, 2018 - updated README.md ( G Chirokova ) ! - v2.3.0 - Nov, 2018 (G. Chirokova) ! - renamed to add *module* to a filename ! - added option to use differnt resolution, variable res ! - added subroutine get_ncodaclimo_profile that ! can be called from eohcadd.f90 ! - v2.3.1 - Nov, 2018 (G. Chirokova) ! - fixed issues with renaming variables in filter2 sub ! - this version works again with less than 26 climo variables ! - v2.3.2 - (GC) - redirected all output to log file, added 2nd log file ! - v2.3.3 - (GC) - fixed bug with interplation climo values between December and ! January ! - v2.3.4 - (GC) - fixed bug with error handling. Tested that if no climo is ! avaialble: ! - climo var names written to output file ! - all climo values set to missing ! - code will stop for devlop verion and keep running ! for real-time version ! - v2.3.5 - (GC) - removed filter2 subroutine ! - v2.3.6 - (GC) - replaced global variable n_data_arr by ncvars and made ! it local. That fixes bug with eohcadd and small (==3) ! number of climo variables ! - v2.3.6 - (GC) - renamed n_climos_data_arr_max to ncvars_max ! - v2.4.0 - (SS) - removed ioper flag ! - v2.4.1 - (MD) - commas removed from legacy write statements (Jan 2024) implicit none private public :: ncodaclimo2 public :: filter_out_bad_climo public :: get_ncodaclimo_profile public :: get_final_ncodaclimo_profile !data grid parameters integer,parameter :: n_lons_reg_max = 1440 integer,parameter :: n_lats_reg_max = 522 !use with 0.25 deg resolution !integer,parameter :: n_lons_reg = 1440 !integer,parameter :: n_lats_reg = 522 !use with 0.50 deg resolution !integer,parameter :: n_lons_reg = 720 !integer,parameter :: n_lats_reg = 261 integer,parameter,public :: ncvars_max = 26 integer,parameter :: n_months = 12 ! current number of used climo variables requested by user !integer :: ncvars character(len=4),public :: var_names_arr_nclimo2(ncvars_max) character(len=4) :: all_var_names_arr(ncvars_max) real, allocatable :: data_2dr_arr(:,:,:,:) character(len=27) :: climo_fnames(ncvars_max) integer,parameter :: climo_n_months = 12 integer :: climo_nlons,climo_nlats real :: climo_dlon,climo_dlat, climo_slon, climo_slat integer :: climo_syear, climo_eyear real :: climo_rmiss character(len=25) :: climo_data_source logical :: climo_open = .false. real,parameter,private :: rmiss = -999.9 integer,parameter,private :: imiss = 9999 contains subroutine get_final_ncodaclimo_profile(mft,ilon,ilat, adjust_climo,& climo_var_names_arr,iyear4,imonth,iday,resolution,& ncvars,ncvars_lsdiag,& ierr_ncoda_climo2,& climo_profile_2d_arr,& lu_log_info,lu_log_debug,lu_climo,develop) ! gets interpolated and filtered ncodaclimo profile for ! all forecast times and for all ncoda_climo2 variables at a single ! point implicit none ! input variables -------------------------------------------- ! query year integer :: iyear4 ! query month integer :: imonth ! query day integer :: iday ! note: lonv can be 0 to 360 or ! 0 to 180 for east longitude, and -180 to 0 for west longitude !number of forecast times integer :: mft ! integer lat array for lsdiag integer:: ilat(-2:mft) ! integer lon array for lsdiag integer:: ilon(-2:mft) ! set to ! true : filter out bad values from climo ! false : do not filter out bad values from climo ! SHOULD be set to true in most cases ! setting this to 'false' can be only useful for ! code debugging logical :: adjust_climo !integer :: max number of climo variables ! do NOT change, must be alwasy = 26 integer :: ncvars ! ! number of climo variables to write to lsdiag file ! use only variables you need. 1st ncvars_lsdiag ! variables will be read and written to lsdiag ! refer to climo_var_names_arr for the names of the 1st 19 variables ! and rules for changing them integer :: ncvars_lsdiag ! climo var names character(len=4) :: climo_var_names_arr(ncvars) ! resolution: ! 'r025' - use .25 deg resolution for climo ! 'r025' - use .50 deg resolution for climo character(len=4) :: resolution ! sum of all error returned by ncoda_climo2 ! 0 - normal output ! >0 - look inside subroutine get_ncodaclimo_profile.f90 for ! specific error !TODO: make sure these errors are output to the log file integer :: ierr_ncoda_climo2 ! current forecast time integer:: iftime ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! logical unit for reading climo data integer :: lu_climo ! ! output variables ------------------------------------------------ ! ncoda climo profile for a given point for all forecast times real :: climo_profile_2d_arr(ncvars_lsdiag,-2:mft) ! !local variables ------------------------------------------------ real :: tlat_inp ! query latitude (deg n) real :: tlon_inp ! query longitude (deg e) ! ! ------------------------------------------------------------------------------------ ! - populate array climo_profile_2d_arr(ncvars,-2:mft) ! that array contains all 26 (25 +dsst) ncoda climo variables ! for all forecast times ! for a single lsdiag case ! ------------------------------------------------------------------------------------ ! initialize climo_profile to missing climo_profile_2d_arr(:,:) = rmiss !get interpolated and filtered ncoda_climo2 profile do iftime = -2,mft if (( ilat(iftime) .lt. 9000 ) .and. (ilon(iftime) .lt. 9000 )) then tlat_inp = 0.1*float(ilat(iftime)) tlon_inp = -0.1*float(ilon(iftime)) !if ( loncon .eq. 'W' ) then ! tlon_inp = -0.1*float(ilon(iftime)) !endif !----------------------------------------------- ! get interpolated and filtered ncodaclimo profile !init climo error to zero ! if non-zaero after returning from subroutine, then check ! inside subroutine for individual errors ierr_ncoda_climo2 = 0 call get_ncodaclimo_profile(tlat_inp, tlon_inp, adjust_climo,& climo_var_names_arr,iyear4,imonth,iday,resolution,ncvars_lsdiag,& ierr_ncoda_climo2,climo_profile_2d_arr(:,iftime),& lu_log_info, lu_log_debug,lu_climo,develop) !----------------------------------------------- else climo_profile_2d_arr(:,iftime) = rmiss endif enddo return end subroutine get_final_ncodaclimo_profile subroutine get_ncodaclimo_profile(tlat, tlon, adjust_climo,& climo_var_names_arr_inp,ityr4,itmon,itday,res,ncvars,& ierr_ncoda_climo2,climo_profile_arr,& lu_log_info,lu_log_debug,lu_climo,develop) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Calculates interpolated ncoda climo profile for first ncvars variables in climo_var_names_arr_inp !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit none ! !input variables --------------------------------------------- ! query latitude (deg n) real :: tlat ! query longitude (deg e) !note: lonv can be 0 to 360 or ! 0 to 180 for east longitude, and -180 to 0 for west longitude real :: tlon ! query year integer :: ityr4 ! query month integer :: itmon ! query day integer :: itday ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! logical unit for reading climo data integer :: lu_climo ! number of user-requested climo variables ! - number of variables to use for climatology calculations ! - change this to only use climatology data for the 1st ncvars ! variables ! - max available number of var = 26; that includes ! DSST climo integer :: ncvars ! user_requested names for climo variables ! could be a subset of all available climo variables character(len=4) climo_var_names_arr_inp(ncvars) ! array of climo variable names character(len=4) climo_var_names_arr(ncvars_max) !must be set to true to use d16-d32 and dml logical :: adjust_climo ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! resolution of climo data to use character(len=4) :: res ! 'r025' - quarter degree ! 'r050' - half degree !output variables ----------------------------------------------- ! output array with climo data for all variables ! - ncoda climatology values for query lat/lon/date real :: climo_profile_arr(ncvars) ! local variables ----------------------------------------------- ! iterator integer :: icarr ! error flags ! climo interpolation error ! 0 = normal ! 1 = invalid lon/lat ! 2 = invalid date integer :: ierr_date_space ! IO climo error ! 0 = normal ! > 0 = error open/read climo file integer :: ierr_climo_io ! climo variable names error ! 0 = normal ! 1 = fatal error climo var names not set integer :: ierr_climo_names ! error in climo data: interpolation or ! missing values found ! 0 = normal ! 1 = interpolation error, WARNING ! 2 = missing values found in climo error; FATAL integer :: ierr_climo_data ! error filytering out bad climo data ! 0 = climo depths were adjusted for Tmax and Dfloor ! 1 = either Dfloor var missing; ! unable to adjust climo data ! 2 = either Tmax var missing; ! unable to adjust climo data ! 3 = user requested to not set bad values to missing integer :: ierr_adj ! master error for ncodaclimo ! - = 0 - evrything works as expected ! /= 0 - error wereencountered integer :: ierr_ncoda_climo2 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !!!!!! define climo variables in the calling subroutine !!!!!!! !rules ! 1) - always define all 26 variables; put in the beginning the ! one you actually need ! 2) - only the 1st ncvars variables will be used, and only ! data frm ncvar climo files will be read ! 3) if use any of dml, d32-d16, d_tmax, then MUST ! use d_floor ! 4) if use any of d32 - d16 then MUST use Tmax ! ! List of all available climo variable names ! do NOT uncomment. Climo names must be defined in teh calling subroutine ! climo_var_names_arr(1) = 'XLON' ! climo_var_names_arr(2) = 'XLAT' ! climo_var_names_arr(3) = 'XDST' ! climo_var_names_arr(4) = 'XNST' ! climo_var_names_arr(5) = 'XDFR' ! climo_var_names_arr(6) = 'XTMX' ! climo_var_names_arr(7) = 'XDTX' ! climo_var_names_arr(8) = 'XDML' ! climo_var_names_arr(9) = 'XD30' ! climo_var_names_arr(10) = 'XD28' ! climo_var_names_arr(11) = 'XD26' ! climo_var_names_arr(12) = 'XD24' ! climo_var_names_arr(13) = 'XD22' ! climo_var_names_arr(14) = 'XD20' ! climo_var_names_arr(15) = 'XD18' ! climo_var_names_arr(16) = 'XD16' ! climo_var_names_arr(17) = 'XTFR' ! climo_var_names_arr(18) = 'XOHC' ! climo_var_names_arr(19) = 'XO20' ! ! climo_var_names_arr(20) = 'XD32' ! climo_var_names_arr(21) = 'XT05' ! climo_var_names_arr(22) = 'XT10' ! climo_var_names_arr(23) = 'XT16' ! climo_var_names_arr(24) = 'XT20' ! climo_var_names_arr(25) = 'XT25' ! climo_var_names_arr(26) = 'XT40' !populate array of varnames with the names provided by user do icarr = 1, ncvars_max climo_var_names_arr(icarr) = climo_var_names_arr_inp(icarr) enddo !call ncodaclimo and filter out bad data call ncodaclimo2(climo_var_names_arr,ncvars,res, & tlat,tlon, & ityr4,itmon,itday,climo_profile_arr,adjust_climo, & ierr_climo_data,ierr_date_space,ierr_climo_io, & ierr_adj,ierr_climo_names,& lu_log_info, lu_log_debug,lu_climo, develop) ierr_ncoda_climo2 = & ierr_date_space+ierr_climo_io+ierr_climo_names+ & ierr_climo_data+ierr_adj !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !write(lu_log_debug,*) "INFO: get_ncodaclimo_profile: & ! &ncodaclimo2 errors ", ierr_climo_data,ierr_date_space, & ! ierr_climo_io !do icarr = 1, ncvars ! write(lu_log_debug,*) "DEBUG: get_ncodaclimo_profile: & ! &", climo_var_names_arr(icarr), & ! ' ',climo_profile_arr(icarr) !enddo return end subroutine get_ncodaclimo_profile subroutine ncodaclimo2(users_var_names_arr,ncvars,& res,tlat,tlon,& ityr,itmon,itday,climo_values,adjust_climo, & ierr_climo_data,ierr_date_space,ierr_climo_io,& ierr_adj,ierr_climo_names,& lu_log_info,lu_log_debug,lu_climo,develop) implicit none ! input -------------------------------------- ! number of climo variables integer :: ncvars ! latitude of the current point real :: tlat ! longitude of the current point real :: tlon ! year of the current point integer :: ityr ! month of the current point integer :: itmon ! day of the current point integer :: itday ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! logical unit for reading climo data integer :: lu_climo ! resolution of climo data character(len=4) :: res ! number of climo longitudes integer :: n_lons_reg ! number of climo latitudes integer :: n_lats_reg ! climo var names provided by user. THat might be a subset of all ! available climo var names character(len=4) :: users_var_names_arr(ncvars_max) ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! flag to replace bad missing climo values by missing logical :: adjust_climo ! output ------------------------------------------- ! interpolated ncoda climo profile with bad data filtered out real :: climo_values(ncvars) ! IO climo error integer :: ierr_climo_io ! IO climo interpolation integer :: ierr_climo_data ! IO climo interpolation integer :: ierr_date_space ! error getting climo names integer :: ierr_climo_names ! error filtering out bad climo data integer :: ierr_adj ! local ------------------------------------------------ ! local copy of climo variable names character(len=4) :: var_names_arr_local(ncvars_max) ! iterator integer :: i !set all output values to missing climo_values(:) = rmiss ierr_climo_data = imiss ierr_date_space = imiss ierr_adj = imiss ierr_climo_names = imiss ! define number of lats and lons based on resolution ! res = 'r050' ! res = 'r025' if ( res .eq. 'r050') then n_lons_reg = 720 n_lats_reg = 261 else if ( res .eq. 'r025') then n_lons_reg = 1440 n_lats_reg = 522 endif if ( .not. climo_open ) then if ( allocated (data_2dr_arr) ) then deallocate (data_2dr_arr) endif allocate (data_2dr_arr(n_lons_reg,n_lats_reg,0:n_months+1,ncvars_max)) data_2dr_arr(:,:,:,:) = rmiss ierr_climo_io = imiss call ncoda_climo_readall(ncvars,res,n_lons_reg,n_lats_reg,users_var_names_arr,ierr_climo_names,& ierr_climo_io,& lu_climo,lu_log_info, lu_log_debug,develop) if (ierr_climo_io .eq. 0 ) then climo_open = .true. ierr_climo_names = 0 else climo_open = .false. 800 format(/,'CRITICAL ERROR: ncodaclimo2: & &climo files are not open in ncodaclimo2 on return from ncoda_climo_readall') write(lu_log_info,800) ierr_climo_names = 1 ! stop for developmental version; ! keep running for real-time version if ( develop ) then stop else return endif endif endif ! if ( climo_open ) then ierr_climo_names = 0 call interp_ncodaclimo_all(ncvars,res,n_lons_reg,n_lats_reg,tlat,tlon,& ityr,itmon,itday,& climo_values,ierr_climo_data,ierr_date_space,& lu_log_info, lu_log_debug, develop) endif if ( adjust_climo ) then var_names_arr_local = var_names_arr_nclimo2 call filter_out_bad_climo(var_names_arr_local,climo_values,ncvars,ierr_adj,& lu_log_info, lu_log_debug,develop) else 804 format(/,'WARNING: ncodaclimo2: & & bad climo values for d16 -d31, dml, and dtmax were not set to missing & & per user request ') write(lu_log_info,804) ierr_adj = 0 endif return end subroutine ncodaclimo2 subroutine ncoda_climo_readall(ncvars,res,n_lons_reg,n_lats_reg,users_var_names_arr,& ierr_climo_names,ierr_all_climo,& lu_climo,lu_log_info,lu_log_debug,develop) !use vars_ncodaclimo implicit none !input ------------------------- ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! logical unit for reading climo data integer :: lu_climo ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! number of climo lons integer :: n_lons_reg ! number of climo lats integer :: n_lats_reg ! resolution of climo data character(len=4) :: res ! array of names for user-requested climo variables ! that can be a subset of all available climo variables ! see test_ncodaclimo for more details character(len=4) :: users_var_names_arr(ncvars_max) ! current number of used climo variables requested by user integer :: ncvars ! output --------------------------------- ! error flag for getting climo names integer :: ierr_climo_names ! master error for all error related to climo names and IO ! 0 = normal ! 1 = error open/read climo file integer :: ierr_all_climo ! local --------------------------------- ! filename for climo data character(len=200) :: fn_climo ! filename for climo data with path character(len = 200) :: fn_climo_path_local ! temporary array of climo data for eacg lon,lat, month real :: cdata(n_lons_reg,n_lats_reg,0:n_months+1) ! current climo variable name character(len=4) :: cvar ! master error for reading climo data ! will be non-zero if ! - cannot open climo, file ! - cannot read climo file header ! - cannot read climo file data integer :: ierr_climo_read_data ! iterator integer :: iarr ! initialize errors to missing ierr_climo_names = imiss ierr_climo_read_data = imiss ierr_all_climo = imiss !The full list of currently available variables for reference !The actual list is created in the calling subroutine outside !of the module ncoda_clim2 !all_var_names_arr(1) = 'XLON' !all_var_names_arr(2) = 'XLAT' !all_var_names_arr(3) = 'XDST' !all_var_names_arr(4) = 'XNST' !all_var_names_arr(5) = 'XTMX' !all_var_names_arr(6) = 'XDTX' !all_var_names_arr(7) = 'XDML' !all_var_names_arr(8) = 'XD32' !all_var_names_arr(9) = 'XD30' !all_var_names_arr(10) = 'XD28' !all_var_names_arr(11) = 'XD26' !all_var_names_arr(12) = 'XD24' !all_var_names_arr(13) = 'XD22' !all_var_names_arr(14) = 'XD20' !all_var_names_arr(15) = 'XD18' !all_var_names_arr(16) = 'XD16' !all_var_names_arr(17) = 'XDFR' !all_var_names_arr(18) = 'XTFR' !all_var_names_arr(19) = 'XOHC' !all_var_names_arr(20) = 'XO20' !all_var_names_arr(21) = 'XT05' !all_var_names_arr(22) = 'XT10' !all_var_names_arr(23) = 'XT16' !all_var_names_arr(24) = 'XT20' !all_var_names_arr(25) = 'XT25' !all_var_names_arr(26) = 'XT40' do iarr = 1,ncvars_max var_names_arr_nclimo2(iarr) = 'NONE' enddo do iarr = 1,ncvars_max var_names_arr_nclimo2(iarr) = users_var_names_arr(iarr) enddo !check variable names do iarr = 1,ncvars if (var_names_arr_nclimo2(iarr) .eq. 'NONE' ) then 940 format(/,'CRITICAL ERROR: ncoda_climo_readall: & &climo var names not set in ncoda_climo_readall.') write(lu_log_info,940) ierr_climo_names = 1 ierr_all_climo = ierr_climo_names ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif else ierr_climo_names = 0 endif enddo ierr_all_climo = ierr_climo_names !open and read into a single array all climo data do iarr = 1, ncvars cvar = var_names_arr_nclimo2(iarr) if (res .eq. 'r025' ) then fn_climo_path_local = './climo25/' if ( cvar .eq. 'XDST' ) then fn_climo = 'clim_'//cvar//'_dsst_025deg.dat' else fn_climo = 'clim_'//cvar//'_eohc_025deg.dat' endif else if (res .eq. 'r050' ) then fn_climo_path_local = './climo50/' if ( cvar .eq. 'XDST' ) then fn_climo = 'clim_'//cvar//'_dsst_050deg.dat' else fn_climo = 'clim_'//cvar//'_eohc_050deg.dat' endif endif call read_ncodaclimo(res,n_lons_reg,n_lats_reg,& fn_climo,fn_climo_path_local,lu_climo,& cdata,ierr_climo_read_data,& lu_log_info, lu_log_debug,develop) data_2dr_arr(:,:,:,iarr) = cdata ierr_all_climo = ierr_all_climo+ierr_climo_read_data enddo return end subroutine ncoda_climo_readall subroutine read_ncodaclimo(res,n_lons_reg,n_lats_reg,& fn_climo,fn_climo_path_local,lu_climo,& cdata_read,err_climo,& lu_log_info,lu_log_debug,develop) !input -------------------------------- !use vars_ncodaclimo character (len=80):: fn_climo_path ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop !number of climo lons integer :: n_lons_reg !number of climo lats integer :: n_lats_reg ! full filename of climo file with path character(len=256) :: fn_climo_wpath ! location of coefficients file (climofiles arein th esame place for ! operational runs) character(len=256) :: coef_location ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! resolution of climo data character(len=4) :: res ! location of climo files for non-operational runs character(len = 200) :: fn_climo_path_local ! climo filename without path character(len=200):: fn_climo ! logical unit for reading climo data integer :: lu_climo ! output -------------------------------- ! array of climo data for all lons, lats, months real :: cdata_read(n_lons_reg,n_lats_reg,0:n_months+1) ! IO climo master error integer :: err_climo ! local ---------------------------------- ! name of the current limo variable character(len=4) :: cvar ! array of month names character (len=3):: cmon(12) data cmon /'JAN','FEB','MAR','APR','MAY','JUN',& 'JUL','AUG','SEP','OCT','NOV','DEC'/ ! iterator for reading climo data integer :: ilats ! iterator for reading climo data integer :: ilons ! iterator for reading climo data integer :: imonth ! IO climo error opening file integer :: err_openclimo ! IO climo error reading header integer :: err_header_climo ! IO climo errors reading data integer :: err_data_climo ! initialize errors err_openclimo = imiss err_header_climo = imiss err_data_climo = imiss err_climo = imiss ! read climo data in ascii format !get filename and location call getenv( "SHIPS_COEF", coef_location ) fn_climo_wpath = trim( coef_location ) // trim( fn_climo ) write(lu_log_debug,*) "INFO: read_ncodaclimo: & &Reading ascii climo data from ", fn_climo_path_local write(lu_log_info,*) "INFO: read_ncodaclimo: & &Reading ascii climo data from ", fn_climo_wpath !open data open(unit=lu_climo,file=fn_climo_wpath,status='old',iostat=err_openclimo) 950 format(/,'CRITICAL ERROR: read_ncodaclimo: & &in read_ncodaclimo unable to open climo file: ',a200) write(lu_log_debug,*) "DEBUG: read_ncodaclimo: open ", err_openclimo if (err_openclimo /= 0) then write(lu_log_info,950) fn_climo ! stop code for developmental version ! but keep running for real-time verion err_climo = err_openclimo if ( develop ) then stop else return endif endif !read header read(lu_climo,'(a11,2x,i6,i7,2(i7,2f9.3),f9.3,a28)',iostat=err_header_climo) & cvar,climo_syear, climo_eyear, climo_nlons, climo_slon, & climo_dlon, climo_nlats, climo_slat, climo_dlat,& climo_rmiss,climo_data_source 1001 format(/,'ERROR: read_ncodaclimo: & &Error in read_ncodaclimo while reading header from ', i10,a60) write(lu_log_debug,*) "DEBUG: read_ncodaclimo: read header ", err_header_climo if (err_header_climo .ne. 0 ) then write(lu_log_info,1001) err_header_climo,fn_climo err_climo = err_header_climo if ( develop ) then stop else return endif endif write(lu_log_info,*) 'INFO: read_ncodaclimo: CVAR ', cvar write(lu_log_info,*) 'INFO: read_ncodaclimo: & &Header data for variable ', cvar,climo_syear, climo_eyear, climo_nlons,& climo_slon,climo_dlon, climo_nlats, climo_slat,& climo_dlat,climo_rmiss,climo_data_source !set all climo data to missing cdata_read(:,:,:) = climo_rmiss !read climo data write(lu_log_debug,*) "INFO: read_ncodaclimo: res ", res write(lu_log_debug,*) "INFO: read_ncodaclimo: res ", & &n_lons_reg,n_lats_reg,climo_nlons,climo_nlats do imonth = 1, climo_n_months read(lu_climo,'(a3)',iostat=err_data_climo) cmon(imonth) do ilats = climo_nlats, 1, -1 if (res .eq. 'r025' ) then read(lu_climo,'(1440(f7.1))',iostat=err_data_climo) & (cdata_read(ilons,ilats,imonth),ilons=1,climo_nlons) else if (res .eq. 'r050' ) then read(lu_climo,'(720(f7.1))',iostat=err_data_climo) & (cdata_read(ilons,ilats,imonth),ilons=1,climo_nlons) endif if (err_data_climo > 0 ) then write(lu_log_info,1003) err_data_climo,fn_climo 1003 format(/,'CRITICAL ERROR: read_ncodaclimo: & &Error in read_ncodaclimo while reading data from ', i10,a60) err_climo = err_data_climo ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif else if (err_data_climo < 0 ) then write(lu_log_info,1004) err_data_climo,fn_climo 1004 format(/,'INFO: read_ncodaclimo: & &EOF reached while reading data from ', i10,a60) ! set err_climo to 0 - this is normal code flow err_climo = 0 endif enddo if (imonth .eq. 1) cdata_read(:,:,climo_n_months+1)=cdata_read(:,:,imonth) if (imonth .eq. 12) cdata_read(:,:,0)=cdata_read(:,:,imonth) enddo close(lu_climo) ! at this point the sum of errors must be zero err_climo = err_openclimo+err_header_climo+err_data_climo if (err_climo .ne. 0 ) then write(lu_log_info,1002) err_climo,fn_climo 1002 format(/,'ERROR: read_ncodaclimo: & &IO errors encountered while executing read_ncodaclimo ', i10,a60) ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif endif return end subroutine read_ncodaclimo subroutine interp_ncodaclimo_all(ncvars,res,n_lons_reg,n_lats_reg,& tlat,tlon_inp,ityr,itmon,itday,& climo_values,ierr_climo,ierr_date_space,& lu_log_info,lu_log_debug,develop) !use vars_ncodaclimo implicit none ! input ------------------------------ ! query latitude - could be 0 - 360 or -180-180 real,intent(in) :: tlat ! query longitude real,intent(in) :: tlon_inp ! query year integer :: ityr ! query month integer :: itmon ! query day integer :: itday ! number of climo longitudes integer,intent(in) :: n_lons_reg ! number of climo latitudes integer,intent(in) :: n_lats_reg ! resolution of climo data character(len = 4) :: res ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! current number of used climo variables requested by user integer :: ncvars ! output ----------------------------------- ! TODO: only retyrn master error to calling subroutine ! error - incorrect date, lat, lon integer :: ierr_date_space ! master error for climo interpolation integer :: ierr_climo !array containing interpolated climo values for a given !lat,lon,date point real :: climo_values(ncvars) ! local ------------------------------------ ! local copy of query longitude - will be converted to 0 - 360 real :: tlon ! arrey of interpolation errors integer :: ierr_interp_arr(ncvars) ! arrey of missing values errors (ncoda climatology ! doe NOT have missing values) integer :: ierr_miss_arr(ncvars) !ndays in a month integer :: ndays_mon(12) ! sum of all error for missing values integer :: ierr_all_climo_miss ! sum of all error for interpolation integer :: ierr_all_climo_interp !iterator integer :: iarr ! climo filename without path character(len = 50) :: fn_climo ! maximum latitude available for climo data real :: max_lat ! maximum longitude available for climo data real :: max_lon ! variable name for current climo variable character(len=4) :: var_name ! missing values used by ncoda_climo (should not be encountered) real,parameter :: rmiss = -999.9 ! initialize all arrors to missing ierr_date_space = imiss ierr_climo = imiss ierr_all_climo_miss = imiss ierr_all_climo_interp = imiss !set output climo data to missing climo_values(:) = rmiss ! populate ndays_mon array ndays_mon(1) = 31 ndays_mon(2) = 28 ndays_mon(3) = 31 ndays_mon(4) = 30 ndays_mon(5) = 31 ndays_mon(6) = 30 ndays_mon(7) = 31 ndays_mon(8) = 31 ndays_mon(9) = 30 ndays_mon(10) = 31 ndays_mon(11) = 30 ndays_mon(12) = 31 ! fix febfuary ndays for leap year if ( leap_year(ityr) ) then ndays_mon(2) = 29 endif ! make sure lon is between 0 - 360 before calling interpolation ! 0-360 is the convention used by llintohc if (tlon_inp < 0 ) then tlon = 360.+tlon_inp else if (tlon_inp .ge. 360 ) then tlon = tlon_inp - 360. else tlon = tlon_inp endif !check valid date and lat lon max_lat = climo_slat+climo_dlat*(climo_nlats-1) max_lon = climo_slon+climo_dlon*(climo_nlons-1) if ( ((itmon .lt. 1) .or. (itmon .gt. 12) ) .or. & ((itday .lt. 1) .or. (itday .gt. ndays_mon(itmon)) ) )then ierr_date_space = 2 else if (((tlat .lt. climo_slat) .or. (tlat .gt. max_lat) ) .or. & ((tlon .lt. climo_slon) .or. (tlon .gt. max_lon) )) then ierr_date_space = 1 else ierr_date_space = 0 endif !return if bad date if (ierr_date_space .eq. 2) then write(lu_log_info,1002) itmon,itday 1002 format(/,'CRITICAL ERROR: interp_ncodaclimo_all: & & BAD time in interpolation for itmon,itday ', i10,i10) ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif endif !continue with WARNING if bad lon or lat ! llintsp should take care of that but it might ! be helpful to chack why bad lat/lon is encountered if (ierr_date_space .eq. 2) then write(lu_log_info,1003) tlon,tlat 1003 format(/,'WARNING: interp_ncodaclimo_all: & &Error space interpolation for tlon,tlat ',f9.3,f9.3) endif do iarr = 1, ncvars var_name = climo_fnames(iarr) if (res .eq. 'r025' ) then fn_climo = 'clim_'//var_name//'_eohc_025deg.dat' else if (res .eq. 'r050' ) then fn_climo = 'clim_'//var_name//'_eohc_050deg.dat' endif call interp_ncodaclimo(tlat,tlon,ityr,itmon,itday,& climo_nlons,climo_nlats,climo_n_months,& climo_dlon,climo_dlat,climo_slon,climo_slat,& var_names_arr_nclimo2(iarr),data_2dr_arr(:,:,:,iarr),& climo_rmiss,climo_values(iarr),& ierr_interp_arr(iarr),ierr_miss_arr(iarr),& lu_log_info, lu_log_debug,develop) enddo ! add all errors from interpolate_ncodaclimo ierr_all_climo_miss = sum(ierr_miss_arr) ierr_all_climo_interp = sum(ierr_interp_arr) ! continue or return depending on error value if ( ierr_all_climo_interp /= 0 ) then ierr_climo = 1 !error were encountered with interpolation and logged write(lu_log_info,801) 801 format(/,'WARNING: interp_ncodaclimo_all: & &interp_ncodaclimo exited with errors ') elseif ( ierr_all_climo_miss /= 0 ) then ierr_climo = 2 !fatal - the climatology does not have missing values !this condition indicates code did not run correctly write(lu_log_info,802) 802 format(/,'CRITICAL ERROR: interp_ncodaclimo_all: & &non-existing missing values fount in climo ') ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif else ierr_climo = 0 endif return end subroutine interp_ncodaclimo_all subroutine interp_ncodaclimo(tlat,tlon,ityr,itmon,itday,& nxmax,nymax,nmon,dlon,dlat,slon,slat,cvar,cdat,rmiss,cvalue,ierr_miss,ierr_int,& lu_log_info,lu_log_debug,develop) !returns climatological vales for a given lat, lon ! the climatological cvalue is interpolated in time between ! 2 monthly means ! ! implicit none ! input ------------------------------------------------- ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! query year integer :: ityr ! query month integer :: itmon ! query day integer :: itday ! query lat real :: tlat ! query lon real :: tlon ! logical unit numberfor regular log file integer :: lu_log_info ! logical unit numberfor debug log file integer :: lu_log_debug ! climo bottom (south) starting lat real:: slat ! climo left (west) starting lon real:: slon ! climo lon increment real:: dlon ! climo lat increment real:: dlat ! max possible number of lons integer :: nxmax ! max possible number of lats integer :: nymax ! number of months in a year integer :: nmon ! climatology for all lats, lons, months real:: cdat(nxmax,nymax,0:nmon+1) ! variable name for current climo array character (len=4) :: cvar ! real missing value for climo real :: rmiss !output --------------------------------------------------- ! interpolation error integer :: ierr_int ! missing values in climatology error: climatology does not have missing ! values. If missig valuyes are returned then something is wrong integer :: ierr_miss ! climatology value for single point and singleclimatology value for ! single point and single variable real :: cvalue ! local --------------------------------------------------- ! integer missing value for climo integer,parameter :: imiss = 9999 ! temporary month number integer :: ipoint2 ! 2nd temporary month number integer :: ipoint1 ! julian day integer :: jd ! array of midpoints (in days) for each month real :: midpt(0:nmon+1) ! temporary climatology value 1 real :: value1 ! temporary climatology value 2 real :: value2 ! weight1 for time difference from month mid-point real :: wt1 ! weight2 for time difference from month mid-point real :: wt2 ! temporary climatology array for a given time real:: dummy(nxmax,nymax) ! TODO: move that to calling sub !zonal periodicity for llinsp integer,parameter :: izp = 1 ! iterator integer::i ! iterator integer::ii ! iterator integer::jj ! iterator integer::j ! initialize errors to missing ierr_miss = imiss !will be set to 1 if missing values found ierr_int = imiss !will be reset is llinsp runs normally ! initialize temporary climo values to missing value1 = rmiss value2 = rmiss ! midpoint for each month midpt(0)=-14.5 midpt(1)=15.5 midpt(2)=45. midpt(3)=74.5 midpt(4)=105. midpt(5)=135.5 midpt(6)=166. midpt(7)=196.5 midpt(8)=227.5 midpt(9)=258. midpt(10)=288.5 midpt(11)=319. midpt(12)=349.5 midpt(13)=380.5 !correct for leap year if ( leap_year(ityr) ) then midpt(2) = 45.5 do i = 3,12 midpt(i) = midpt(i)+1. enddo endif ! calculate julian day corresponding to the input date call jday(itmon,itday,ityr,jd) !find month correspoding to the current date do i=1,nmon+1 if(jd .ge. midpt(i-1) .and. jd.lt.midpt(i)) ipoint2=i end do ! interpolate data for the current date from the ! climatology for this month and previous month !weights: time difference from the month mid-point wt1=0.0 wt2=0.0 ! index for month1 and month2 ipoint1 = ipoint2-1 ! calculate time weights wt2=(jd-midpt(ipoint1))/(midpt(ipoint2)-midpt(ipoint1)) wt1=1-wt2 ! save temporary climo array for single time for 1st month for interpolation dummy=cdat(:,:,ipoint1) !spatial interpolation for month1 CALL llintsp(dummy,slon,slat,dlon,dlat,nxmax,nymax,nxmax,nymax,value1,tlon,tlat,izp,ierr_int) if (ierr_int .ne. 0 ) then write(lu_log_info,1003) cvar, tlon,tlat,ityr,itmon,itday 1003 format(/,'WARNING: interp_ncodaclimo: & &in interp_ncodaclimo exited with errors: tlon,tlat,ityr,itmon,itday ',& a60, f9.3,f9.3,i5,i5,i5) endif ! save temporary climo array for single time for 2nd month for interpolation dummy=cdat(:,:,ipoint2) !spatial interpolation for month2 CALL llintsp(dummy,slon,slat,dlon,dlat,nxmax,nymax,nxmax,nymax,value2,tlon,tlat,izp,ierr_int) if (ierr_int .ne. 0 ) then write(lu_log_info,1003) cvar, tlon,tlat,ityr,itmon,itday endif ! find interpolated cvalue if (( value1 .gt. rmiss+1) .and. (value2 .gt. rmiss+1) ) then cvalue=value1*wt1+value2*wt2 ierr_miss = 0 else cvalue=rmiss write(lu_log_info,1002) cvar, tlon,tlat,ityr,itmon,itday ierr_miss = 1 1002 format(/,'CRITICAL ERROR: interp_ncodaclimo: while running interp_ncodaclimo & & missing values found for tlon,tlat,ityr,itmon,itday ',& &a60, f9.3,f9.3,i5,i5,i5) ! stop code for developmental version ! but keep running for real-time verion if ( develop ) then stop else return endif endif return end subroutine interp_ncodaclimo subroutine jday(imon,iday,iyear,julday) ! this routine calculates the julian day (julday) from ! the month (imon), day (iday), and year (iyear). the ! appropriate correction is made for leap year. implicit none ! input -------------------------------- ! query year integer :: iyear ! query month integer :: imon ! query day integer :: iday ! output -------------------------------- ! query julian day integer :: julday ! local ---------------------------------- ! array of number of days in each month integer :: ndmon(12) ! iterator integer :: i ! specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 ! correct for leap year !call check_leap_year(iyear,leap_year) if ( leap_year(iyear) ) then ndmon(2) = 29 endif ! check for illegal input if (imon .lt. 1 .or. imon .gt. 12) then julday=-1 return endif if (iday .lt. 1 .or. iday .gt. ndmon(imon)) then julday=-1 return endif ! calculate the julian day julday = iday if (imon .gt. 1) then do i=2,imon julday = julday + ndmon(i-1) end do endif return end subroutine jday function leap_year(iyear) implicit none integer :: iyear logical :: leap_year if ( (mod(iyear, 4) .eq. 0) .and. (mod(iyear, 100) .ne. 0) ) then leap_year = .true. else if ( (mod(iyear, 4) .eq. 0) .and. (mod(iyear, 100) .eq. 0) .and. (mod(iyear, 400).eq. 0) ) then leap_year = .true. else leap_year = .false. endif return end function leap_year subroutine filter_out_bad_climo(var_names_arr_local,climo_values,ncvars,ierr_adj,& lu_log_info,lu_log_debug,develop) ! this subroutine replaces non-phisical climo values by missing ! that is needed for d32 - d16, dml, dtmx, tmx to make sure ! that depths in the profile are never below the ocaen floor, ! and temperatures do not exceed maximum profile temperature ! !set to missing d16 - d32 if above t_max !or below d_floor implicit none ! input ---------------------------------------- ! number of climo variables requested by user. ! that can be less than the total available number of ! climo variables = 26 integer,intent(in) :: ncvars ! logical unit number for regular log file integer :: lu_log_info ! logical unit number for debug log file integer :: lu_log_debug ! climo var names character(len=4) :: var_names_arr_local(ncvars_max) ! flag for developmenal vs real-time data ! -- develop = .true. - code will stop on critical error ! = .false. - code will stop on critical error logical :: develop ! input and output ! input: non-filtered climo data profile for a single point ! output: filtered climo data profile for a single point real,intent(inout) :: climo_values(ncvars) ! output --------------------------------------- ! error flag for filtering bad data status ! = 0 if values were adjusted ! = 1 if d_floor and t_max needed for adjustment ! are not available integer :: ierr_adj ! local --------------------------------------- ! real missing values real,parameter :: rmiss = -999.9 ! integer missing values integer,parameter :: imiss = 9999 ! small delta to use with float comparison real,parameter :: delta_t = 0.000001 ! flag to determine if temperatures ara available for filtering bad data logical :: tdepths_used ! flag to determine if depths ara available for filtering bad data logical :: depths_used ! flag to determine if depth of the ocean floor is available logical :: dfloor_available ! flag to determine if tmax is available logical :: tmax_available !iterator integer :: ivar !iterator integer :: ifvar ! max number of variables used by filtering subroutine ! this is smaller than total number of EOHC variables ! since not all variables need to be filtered integer,parameter :: max_num_filter_vars = 13 ! variable names used by filtering subroutine - DO NOT CHANGE character(len=4) :: filter_var_names(max_num_filter_vars) ! array of indices for variables that should be filtered integer :: filter_var_idx(max_num_filter_vars) ! array of temperatures for variables that need to be filtered real :: filter_t_arr(max_num_filter_vars) ! integer :: tmax_idx ! integer :: n__d32_idx ! integer :: n__d30_idx ! integer :: n__d28_idx ! integer :: n__d26_idx ! integer :: n__d24_idx ! integer :: n__d22_idx ! integer :: n__d20_idx ! integer :: n__d18_idx ! integer :: n__d16_idx ! integer :: n_dflr_idx ! integer :: n__dml_idx ! integer :: ndtmax_idx ! initialize error to missing ierr_adj = imiss !DO NOT change the order of variables - ! changing order of variables will break the code filter_var_names(1) = 'XTMX' filter_var_names(2) = 'XDTX' filter_var_names(3) = 'XDML' filter_var_names(4) = 'XD32' filter_var_names(5) = 'XD30' filter_var_names(6) = 'XD28' filter_var_names(7) = 'XD26' filter_var_names(8) = 'XD24' filter_var_names(9) = 'XD22' filter_var_names(10) = 'XD20' filter_var_names(11) = 'XD18' filter_var_names(12) = 'XD16' filter_var_names(13) = 'XDFR' !order of this array MUST match order of ! filter_var_names filter_t_arr(1) = rmiss filter_t_arr(2) = rmiss filter_t_arr(3) = rmiss filter_t_arr(4) = 32. filter_t_arr(5) = 30. filter_t_arr(6) = 28. filter_t_arr(7) = 26. filter_t_arr(8) = 24. filter_t_arr(9) = 22. filter_t_arr(10) = 20. filter_t_arr(11) = 18. filter_t_arr(12) = 16. filter_t_arr(13) = rmiss !initialize all indices to missing do ifvar = 1,max_num_filter_vars filter_var_idx(ifvar) = imiss enddo depths_used = .false. tdepths_used = .false. dfloor_available = .false. tmax_available = .false. !find indices of tmax and all depth values do ivar = 1,ncvars do ifvar = 1,max_num_filter_vars if (var_names_arr_local(ivar) .eq. filter_var_names(ifvar) ) then filter_var_idx(ifvar) = ivar endif enddo enddo !set logical variables based on which data are available !this assumes the exsiting order of variables ! and will be broken if that order is changed if ( filter_var_idx(1) .lt. imiss ) tmax_available = .true. if ( filter_var_idx(13).lt. imiss ) dfloor_available = .true. do ifvar = 2,12 if ( filter_var_idx(ifvar) .lt. imiss ) then depths_used = .true. exit endif enddo do ifvar = 4,12 if ( filter_var_idx(ifvar) .lt. imiss ) then tdepths_used = .true. exit endif enddo !flter out d32 - d16 to exclude values below d_floor ierr_adj = 0 if ( ( depths_used ) .and. ( dfloor_available ) ) then do ifvar = 2,12 !loop through ndtmax, ndml, d32 - d16 !filter_var_idx(13) - MUST be d_floor if (filter_var_idx(ifvar) .lt. imiss) then if (climo_values(filter_var_idx(ifvar)) .gt. & climo_values(filter_var_idx(13))) then climo_values(filter_var_idx(ifvar)) = rmiss endif endif enddo ierr_adj = 0 else if ( ( depths_used ) .and. ( .not. dfloor_available ) ) then !set error flag and write log messages if d_floor !is needed for filtering but not available ierr_adj = 1 949 format(/,'ERROR: filter_out_bad_climo: & &Error adjusting climo depths in filter_out_bad_climo.f90 ',& "if using dml,d16-d32,dtmax then MUST also use d_floor") write(lu_log_info,949) endif !filter d32 - d16 based to exclude values above max temp if ( ( tdepths_used ) .and. ( tmax_available ) ) then do ifvar = 4,12 !loop through d32 - d16 !filter_var_idx(1) - MUST be t_max if (filter_var_idx(ifvar) .lt. imiss) then if (climo_values(filter_var_idx(1)) .lt. & (filter_t_arr(ifvar)-delta_t)) then climo_values(filter_var_idx(ifvar)) = rmiss endif endif enddo else if ( ( tdepths_used ) .and. ( .not. tmax_available ) ) then !set error flag and write log messages if t_max !is needed for filtering but not available ierr_adj = 2 950 format(/,'ERROR: filter_out_bad_climo: & &Error adjusting climo depths in filter_out_bad_climo.f90 ',& "if using d16-d32 then MUST also use t_max") write(lu_log_info,950) endif return end subroutine filter_out_bad_climo end module ncoda_climo2_module