! This is a function to be called by the WEPS main subroutine mcip_weps ! Description: ! ! This program extracts from MCIP's METCRO2D and METCRO3D files ! variables needed for WEPS for the first 49 hours of output. ! The output is comma-delimited file. ! ! ! implicit none !...include file from the netCDF library !...(www.unidata.ucar.edu/software/netcdf) include 'netcdf.inc' !...include files from the Models-3 I/O API library !...(http://www.baronams.com/products/ioapi) include 'PARMS3.EXT' ! i/o api constants include 'FDESC3.EXT' ! i/o api file description data structure include 'IODECL3.EXT' ! i/o api function declarations include 'p1werm.inc' include 'w1wind.inc' include 'main/w1win.inc' include 'w1clig.inc' include 'airpact/spatialGIS.inc' ! external functions integer, external :: GETEFILE ! an i/o api function ! parameters/datapact/mcip4weps. character(len=namlen3), parameter :: progname='mcip4weps' integer, parameter :: iunit = 12 ! variables to extract from METCRO2D file integer, parameter :: nvar_metcro2d = 6 character(len=namlen3), parameter, dimension(nvar_metcro2d) :: & vname_metcro2d=(/ & & 'TEMP2 ' & & ,'WSPD10' & & ,'WDIR10' & & ,'RGRND ' & & ,'RN ' & & ,'RC ' & /) integer, parameter :: iin_metcro2d_temp2=1, iin_metcro2d_wspd10=2 & ,iin_metcro2d_wdir10=3, iin_metcro2d_rgrnd=4, iin_metcro2d_rn=5 & , iin_metcro2d_rc=6 ! variables to extract from METCRO3D file integer, parameter :: nvar_metcro3d = 3 character(len=namlen3), parameter, dimension(nvar_metcro3d) :: & vname_metcro3d=(/ & & 'TA ' & & ,'QV ' & & ,'PRES' & /) integer, parameter :: iin_metcro3d_ta=1, iin_metcro3d_qv=2 & ,iin_metcro3d_pres=3 ! variables to output integer, parameter :: nvar_out = 6 character(len=namlen3), parameter, dimension(nvar_out) :: & vname_out=(/ & & 'TEMP2m ' & & ,'TDEW ' & & ,'WSPD10m' & & ,'WDIR10m' & & ,'SWDOWN ' & & ,'PRECIP ' & /) character(len=namlen3), parameter, dimension(nvar_out) :: units_out=(/ & & 'oC ' & & ,'oC ' & & ,'m/s ' & & ,'deg ' & & ,'langley/hr' & & ,'mm ' & /) integer, parameter :: iout_temp2m=1, iout_tdew=2,iout_wspd10m=3 & ,iout_wdir10m=4, iout_swdown=5, iout_precip=6 !maximum number of timestep to output ! (WEPS' EROSION uses daily timesteps, so must be multiples of 24 plus 1) integer :: nhour_max ! integer, parameter :: nday=2 ! integer, parameter :: n_grid = 95 ! local scalar variables integer :: i,j,status,ivar,icount,itime integer :: logdev integer :: rdev ! unit number for output file integer :: readTime, readDate integer :: YYYY,DOY,HH,MM,DD,SS,MN real :: temperature, psat, pvap, gamma, RH integer :: k integer :: d !counting dates real :: temp ! CHARACTER (len=256) :: INFILE CHARACTER (len=256) :: MESG ! mesage statement CHARACTER (len=256) :: OUTFILE ! output file ! Local arrays real, allocatable :: inbuff(:,:) , data_metcro2d(:,:,:), data_out(:,:,:,:) real, allocatable :: data_metcro3d(:,:,:) ! real ::temp_max(n_grid,n_grid,nday) ! real ::temp_min(n_grid,n_grid,nday) ! real ::temp_ave(n_grid,n_grid,nday) ! real ::temp_sum(n_grid,n_grid,nday) ! real ::awe_rad(n_grid,n_grid,nday) ! real ::awe_prec(n_grid,n_grid,nday) ! real ::wdsp_hr(n_grid,n_grid,nday,24) real :: temp_hr(n_grid,n_grid,nday,24) ! End header --------------------------------------------------------- ! -------------------------------------------------------------------- ! Get environment variables ! -------------------------------------------------------------------- ! Initialiation of the Models-3 I/O API logdev = init3() ! initialization returns unit # for log ! Set input and output file names ! Set unit number of output file and open the output file ! (GETFILE is a Models-3 I/O API library function) RDEV = GETEFILE( 'OUTFILE', .FALSE., .TRUE., trim(progname) ) ! -------------------------------------------------------------------- ! Read data from file and output to csv file ! -------------------------------------------------------------------- ! open METCRO2D file if ( .not. open3( 'METCRO2D', fsread3, progname ) ) then mesg = "could not open file METCRO2D" call m3exit( progname, 0, 0, mesg, 2 ) endif ! get input file description if ( .not. desc3('METCRO2D') ) then mesg = "could not get desc of METCRO2D" call m3exit( progname, 0 , 0, mesg, 2 ) endif ! Note: After desc3 call above variables such as ncols3d, nrows3d, ! nvars3d, vname3d are populated with values relating to the ! the file content in 'METCRO2D'. ! ! ncols3d = number of columns ! nrows3d = number of rows ! nlays3d = number of vertical layers ! nvars3d = number of variables (excluding TFLAG) ! vnamed3d = name of variables (dimension nvars3d) ! ! See the include file FDESC3D.EXT for a list of all parameters ! associated with file content of an Models-3 I/O API file ! open METCRO3D file if ( .not. open3( 'METCRO3D', fsread3, progname ) ) then mesg = "could not open file METCRO2D" call m3exit( progname, 0, 0, mesg, 2 ) endif nhour_max = mxrec3d ! added in May 12 ! Write "header" to output file " mcip_weps write(rdev,fmt='("GRIDID_SA,I,J,",100(A,","))',advance='no') & (trim(vname_out(ivar)),ivar=1,nvar_out-1) write(rdev,fmt='(A)')trim(vname_out(nvar_out)) write(rdev,fmt='(" , , , ,",100(A,","))',advance='no') (trim(units_out(ivar)),ivar=1,nvar_out-1) write(rdev,fmt='(A)')trim(units_out(nvar_out)) ! initialize local variables do i=1,n_grid do j=1,n_grid do d=1,nday temp_max(i,j,d) = 0.0 temp_min(i,j,d) = 0.0 temp_sum(i,j,d) = 0.0 awe_prec(i,j,d) = 0.0 awe_rad(i,j,d) = 0.0 awe_dew(i,j,d) = 0.0 wu_max(i,j,d) = 0.0 wu_min(i,j,d) = 0.0 do ivar=1,24 temp_hr(i,j,d,ivar) = 0.0 end do end do end do end do ! read in data allocate (inbuff(ncols3d,nrows3d), data_metcro2d(ncols3d,nrows3d,nvar_metcro2d)) allocate (data_metcro3d(ncols3d,nrows3d,nvar_metcro3d), data_out(ncols3d,nrows3d,nvar_out,nhour_max) ) readDate = sdate3d ; readTime = stime3d loop_time : do itime = 1, nhour_max d = itime/24.1 + 1 ! find the number of dates k = mod(itime,24) if (k .eq. 0) k = 24 ! hourly setting ! read in values, one variable at time do ivar = 1, nvar_metcro2d ! read in only the first layer and date/time of interest if ( .not. read3( 'METCRO2D', vname_metcro2d(ivar), 1, readDate, readTime, inbuff ) ) then mesg = "could not read from METCRO2D" call m3exit(progname,0,0,mesg,2) else data_metcro2d(:,:,ivar) = inbuff(:,:) endif end do ! ivar = 1, nvar_metcro2d do ivar = 1, nvar_metcro3d ! read in only the first layer and date/time of interest if ( .not. read3( 'METCRO3D', vname_metcro3d(ivar), 1, readDate, readTime, inbuff ) ) then mesg = "could not read from METCRO3D" call m3exit(progname,0,0,mesg,2) else data_metcro3d(:,:,ivar) = inbuff(:,:) endif end do ! ivar = 1, nvar_metcro3d ! calculate output variables ! 2-meter temperature, convert from Kelvin to Celsius data_out(:,:,iout_temp2m,itime) = data_metcro2d(:,:,iin_metcro2d_temp2) - 273.15 ! 10-meter windspeed [=] m/s data_out(:,:,iout_wspd10m,itime) = data_metcro2d(:,:,iin_metcro2d_wspd10) ! 10-meter wind direction [=] deg data_out(:,:,iout_wdir10m,itime) = data_metcro2d(:,:,iin_metcro2d_wdir10) ! precipitation [=] mm data_out(:,:,iout_precip,itime) = 10.*(data_metcro2d(:,:,iin_metcro2d_rc) + data_metcro2d(:,:,iin_metcro2d_rn)) ! downward radiation at the surface, convert from W/m2 to langley/hr data_out(:,:,iout_swdown,itime) = data_metcro2d(:,:,iin_metcro2d_rgrnd)/ 11.622 ! dew point temperature do j = 1, nrows3d do i = 1, ncols3d if (d .LE. nday) then wdsp_hr(i,j,d,k) = data_metcro2d(i,j,iin_metcro2d_wspd10) ! hourly wind speed data ! find the daily wind direction which represents the max wind spead if (wdsp_hr(i,j,d,k) .GT. wu_max(i,j,d)) then wu_max(i,j,d) = wdsp_hr(i,j,d,k) awe_dir(i,j,d) = data_metcro2d(i,j,iin_metcro2d_wdir10) end if if (wdsp_hr(i,j,d,k) .LT. wu_min(i,j,d)) then wu_min(i,j,d) = wdsp_hr(i,j,d,k) end if ! 1st-layer temperature [=] o temperature = data_metcro3d(i,j,iin_metcro3d_ta) - 273.15 temp_hr(i,j,d,k) = temperature if (k .eq. 1) then temp_max(i,j,d) = temp_hr(i,j,d,1) temp_min(i,j,d) = temp_hr(i,j,d,1) end if if (temp_max(i,j,d) .LT. temperature) temp_max(i,j,d)=temperature if (temp_min(i,j,d) .GT. temperature) temp_min(i,j,d)=temperature temp_sum(i,j,d)= temp_sum(i,j,d)+temperature ! adding daily radiation awe_rad(i,j,d) = awe_rad(i,j,d) +data_out(i,j,iout_swdown,itime) ! adding daily precipitation awe_prec(i,j,d) = awe_prec(i,j,d)+ data_out(i,j,iout_precip,itime) ! 1st-layer saturation vapor pressure (following method used by CMAQ's aerosol_driver.F) [=] Pa psat = 610.94 * EXP( 17.625 * temperature / ( temperature + 243.04 ) ) ! 1st-layer water vapor pressure (following method used by CMAQ's aerosol_driver.F) [=] Pa pvap = data_metcro3d(i,j,iin_metcro3d_pres) * data_metcro3d(i,j,iin_metcro3d_qv) & / ( 0.622015 + data_metcro3d(i,j,iin_metcro3d_qv) ) ! 1st-layer relative humidity; limit to range [0.01,1.] RH = max (0.01, min(1.,pvap/psat)) gamma = 17.271*temperature/(237.7+temperature) + log(RH) data_out(i,j,iout_tdew,itime) = 237.7*gamma/(17.271-gamma) awe_dew(i,j,d) = awe_dew(i,j,d)+data_out(i,j,iout_tdew,itime)/24.0 ! average dew point end if end do end do call nextime( readDate, readTime, tstep3d ) end do loop_time !Insert a loop to print file ! write(*,*) 'Max&Min',icount,i,j,d,temp_max(46,64,1),temp_min(46,64,1),awe_prec(46,64,1) ! Output data in comma delimted file icount = 0 do i = 1, nrows3d do j = 1, ncols3d icount = icount + 1 do d= 1,nday write(rdev,*)icount,i,j,d,temp_max(i,j,d),temp_min(i,j,d), & & awe_prec(i,j,d),(wdsp_hr(i,j,d,ivar),ivar=1,24) end do end do end do ! end insert ! call nextime( readDate, readTime, tstep3d ) ! end do loop_time ! output data ! -------------------------------------------------------------------- ! Exit ! -------------------------------------------------------------------- ! call m3exit(progname,0,0,'',0) end