!################################################################# !WP! !WP! Routine samples specific events in space and time, as read from a list of !WP! events. These events are created from the NOAA ESRL Flask database, and can !WP! thus be matched to observations using the eventnumbers. Co-sampling is based on !WP! a 4-hour window around the observations time and date, and slopes sampling is !WP! currently used for the interpolation. Contact Wouter Peters !WP! (Wouter.Peters@noaa.gov) to get the lists of NOAA ESRL event numbers in the !WP! proper format. Or copy and paste the following lines to a file (remove comment, no header): ! !ASK_01D0 23.18 5.42 2728.0 2000 1 1 13 30 0 9703 !CHR_01D1 1.70 -157.17 3.0 2000 1 1 3 40 0 42746 !SPO_01D0 -89.98 -24.80 2810.0 2000 1 1 5 55 0 129319 !STM_01D0 66.00 2.00 7.0 2000 1 1 11 15 0 134065 !KZD_01D0 44.45 75.57 412.0 2000 1 3 4 20 0 69908 !RPB_01D0 13.17 -59.43 45.0 2000 1 3 19 30 0 108606 !KEY_01D0 25.67 -80.20 3.0 2000 1 3 18 0 0 63211 !MLO_01D0 19.53 -155.58 3397.0 2000 1 3 20 37 0 84745 ! !WP! To turn this routine on, add the following RC items to your tm5.rc file: ! ! output.noaa : T ! input.noaa.filename : noaa.list # get from Wouter ! inputdir.noaa : /p73/co2/input/1x1/ # your noaa.list file location here ! output.noaa.filename: noaa.hdf !WP! !WP! Routine is off as default and will thus not bother other users !WP! !WP! Future improvements: !WP! - window size can be read from rc-file !WP! - tracer list can be specific !WP! - interpolation mode could be an rc-file item !WP! - or all 3 interpolations could be used !WP! - more extensive output on model sample (gph, grid, ...) !WP! !WP! Created June 9th, 2007 by Wouter Peters !WP! !### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module user_output_noaa use GO , only : gol, goErr, goPr, goLabel, goUpCase use dims, only : im, jm, lm, dx, xref, dy, yref, xbeg, xend, ybeg, yend use dims, only : nregions, region_name, itaui,itaue, xcyc use chem_param, only : fscale, ntracet use file_hdf, only : THdfFile, Init, Done, TSds, SetDim, WriteData, WriteAttribute use ParTools #ifdef MPI use mpi_const, only : my_real,mpi_comm_world,ierr,mpi_integer,mpi_character use mpi_comm, only : barrier #endif implicit none public :: init_noaa_events, get_noaa, write_noaa_events, noaa_data logical :: noaa_data=.false. private character(len=100) :: EventlistFilename = '' character(len=100) :: EventlistOutFilename = '' character(len=100) :: InputDir = './' integer :: n_event, n_temp integer,dimension(6) :: idate_event integer :: itau_event integer :: window = 4*3600 ! 2 hrs around obs hour type(THdfFile) :: io_hdf ! contains HDF output file character(len=*), parameter :: mname = 'user_output_noaa' type observation integer,dimension(6) :: idate ! date of observation character(len=24) :: ident ! station identifier real :: lat ! latitude of station real :: lon ! longitude of station real :: height ! height of station integer :: eventnumber ! unique reference to NOAA ESRL database real,dimension(ntracet) :: mmix=0.0 integer :: itau=0 integer :: n_accumulate=0 ! nr of times concentration was added in get_noaa integer :: region integer :: is integer :: js integer :: lsmin integer :: lsmax end type observation type(observation), dimension(:), allocatable :: concentration ! private integer,parameter :: nf_trace = 1 real :: rmf ! number of locations to be calculated for 1 model time integer,dimension(nf_trace) :: if_trace = (/ 1 /) contains subroutine write_noaa_events(status) !===================================== ! HDF output for noaa events ! (averaged model results) !===================================== use ParTools use Dims, only : itaui, itaue, ndyn_max, idate, idatei, idatee use file_hdf, only : Init, WriteAttribute, WriteData, TSds, Done use GO, only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile, outdir use chem_param, only : ntracet,names implicit none ! i/o integer :: status integer :: i_stat, tag,n type(TrcFile) :: rcF character(len=12) :: xname integer :: i, ii, j type(TSds), dimension(:), allocatable :: Sds type(TSds) :: ds real,dimension(:),allocatable :: data_array integer, dimension(:,:), allocatable :: iidate character(len=10) :: sidatei, sidatee ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/write_noaa_events' call evaluate_noaa() if(myid==root) then !WP! Create output file and write header call Init( rcF, rcfile,status) call ReadRc( rcF, 'output.noaa.filename', EventlistOutFilename,status) call Done( rcF ,status) IF_NOTOK_RETURN(status=1) !WP! Add the date to the filename i = index(EventlistOutFilename,'.hdf') write(sidatei,'(i4.4,3i2.2)') idatei(1:4) write(sidatee,'(i4.4,3i2.2)') idatee(1:4) if(i /= 0) then EventlistOutFilename = EventlistOutFilename(1:i-1)//'_'//sidatei//'_'//sidatee//'.hdf' else EventlistOutFilename = trim(EventlistOutFilename)//'_'//sidatei//'_'//sidatee//'.hdf' endif call Init( io_hdf, trim(OutDir)//'/'// trim(EventlistOutFilename),'create' , status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_lat' ,(/n_event/),'real(4)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteAttribute(ds, 'Unit', 'deg N',status) call WriteData( ds, concentration(1:n_event)%lat,status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_lon' ,(/n_event/),'real(4)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteAttribute(ds, 'Unit', 'deg E',status) call WriteData( ds, concentration(1:n_event)%lon,status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_height' ,(/n_event/),'real(4)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteData( ds, concentration(1:n_event)%height,status) call WriteAttribute(ds, 'Unit', 'meters a.s.l.',status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_region' ,(/n_event/),'integer(2)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteData( ds, concentration(1:n_event)%region,status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_is' ,(/n_event/),'integer(2)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteData( ds, concentration(1:n_event)%is,status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'event_js' ,(/n_event/),'integer(2)',status) call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteData( ds, concentration(1:n_event)%js,status) call Done(ds,status) IF_NOTOK_RETURN(status=1) call Init(ds, io_hdf, 'idate', (/6, n_event/), 'integer(2)',status) call SetDim( ds, 0, 'n_idate ', 'yymmddhhssmm', (/ (j, j=1,6) /) ,status) call SetDim( ds, 1, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteAttribute(ds, 'Unit', 'year month day hour minutes seconds',status) allocate(iidate(6,n_event)) do j=1,n_event iidate(:,j) = concentration(j)%idate enddo call WriteData( ds, iidate,status) IF_NOTOK_RETURN(status=1) deallocate(iidate) call Done(ds,status) IF_NOTOK_RETURN(status=1) call WriteAttribute( io_hdf, 'ntracet', ntracet,status) call WriteAttribute( io_hdf, 'n_event', n_event,status) call WriteAttribute( io_hdf, 'event_ident' , concentration(1:n_event)%ident,status) call WriteAttribute( io_hdf, 'tracer_names' , names(1:ntracet),status) call WriteAttribute( io_hdf, 'event_nsample', concentration(1:n_event)%n_accumulate,status) IF_NOTOK_RETURN(status=1) allocate(Sds(ntracet)) do i=1, ntracet xname = trim(names(i)) call init(Sds(i), io_hdf, xname, (/n_event/), 'real(4)',status) call SetDim( Sds(i), 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status) call WriteAttribute(Sds(i), 'Unit', 'Volume mixing ratio',status) enddo IF_NOTOK_RETURN(status=1) endif ! root !WP! Now write datasets from each PE to file allocate(data_array(n_event)) do i=1,ntracet data_array=concentration(1:n_event)%mmix(i) ! copy tracer to data_array, proper PE will send to root for I/O #ifdef MPI call MPI_BCAST(data_array(1),n_event,MY_REAL,tracer_id(i),mpi_comm_world,ierr) #endif if(myid==root) then call WriteData(Sds(i),data_array,status) IF_NOTOK_RETURN(status=1) call Done(Sds(i),status) IF_NOTOK_RETURN(status=1) endif enddo deallocate(data_array) if(myid==root)deallocate(Sds) deallocate(concentration) if(myid==root) then call Done(io_hdf,status) IF_NOTOK_RETURN(status=1) endif end subroutine write_noaa_events subroutine init_noaa_events(status) use Meteo, only : Set use Meteo, only : temper_dat, humid_dat, gph_dat use binas, only : ae, twopi, grav use dims, only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend use dims, only : nregions, region_name, itaur, xcyc, idate use chem_param, only : ntrace, ntracet, names, fscale use toolbox, only : escape_tm use datetime, only : date2tau use go_string, only : goUpCase use go_string, only : goNum2Str use GO, only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile implicit none ! i/o integer :: status ! local character(len=200) :: dummy integer :: fstatus integer :: i_stat, j_stat,n integer :: sunit0=123 integer :: region integer,dimension(6) :: idate_sim integer :: itau_sim, isim, ispec ! x,y resolution (in degrees) for current region real :: dxr, dyr character(len=8) :: name_spec !CMK changed from 6-->8 logical :: site_found type(TrcFile) :: rcF integer, dimension(6) :: temp_idate integer :: temp_itau,temp_eventnumber real :: temp_lon,temp_lat,temp_height character(len=30) :: temp_ident ! --- const ------------------------------ character(len=*), parameter :: rname = mname//'/init_noaa_events' ! start do region=1,nregions call Set( gph_dat(region), status, used=.true. ) IF_NOTOK_RETURN(status=1) call Set( temper_dat(region), status, used=.true. ) IF_NOTOK_RETURN(status=1) call Set( humid_dat(region), status, used=.true. ) IF_NOTOK_RETURN(status=1) enddo if(myid==root) then call Init( rcF, rcfile,status) call ReadRc( rcF, 'input.noaa.filename', EventlistFilename,status) call ReadRc( rcF, 'inputdir.noaa', InputDir,status) call Done( rcF,status) !======================================================= ! open StationFile (with list of stations) ! "station" means "observation" in this context (arj) !======================================================= open(UNIT=sunit0, FORM='FORMATTED', & FILE=trim(InputDir)//'/'//trim(EventlistFilename)//trim(goNum2Str(idate(1),'(i4.4)'))//trim(goNum2Str(idate(2),'(i2.2)'))//trim(goNum2Str(1,'(i2.2)')) , & ERR=1000) ! count number of stations n_temp=0 stat_loop: DO read(sunit0, '(a)', END=100) dummy n_temp = n_temp+1 end do stat_loop 100 continue print *, '_________________________________________________' & //'_____________________________' print *, ' preparing to read ', n_temp, ' NOAA events from file: ' //trim(InputDir)//'/'//trim(EventlistFilename)//trim(goNum2Str(idate(1),'(i4.4)'))//trim(goNum2Str(idate(2),'(i2.2)'))//trim(goNum2Str(idate(3),'(i2.2)')) print *, '_________________________________________________' & //'_____________________________' endif #ifdef MPI call MPI_BCAST(n_temp ,1, MPI_INTEGER, root ,MPI_COMM_WORLD,ierr) #endif if(n_temp.eq.0) call escape_tm( 'read_noaa_events: no obs in file ' // EventlistFilename) if(allocated(concentration)) deallocate(concentration) ! avoid double allocation allocate(concentration(n_temp)) ! very long array, only partly used dyr = dy/yref(1) dxr = dx/xref(1) !=============== ! read stations !=============== n_event=0 if(myid==root) then rewind(sunit0) do i_stat=1, n_temp read(sunit0, '(a24,1x,f8.2,f8.2,f8.1,1x,6i4,1x,i30)') & temp_ident, & temp_lat, & temp_lon, & temp_height, & temp_idate, & temp_eventnumber call date2tau(temp_idate,temp_itau) ! fill itau if(abs(temp_itau-itaui).lt.1) temp_itau=temp_itau+3600 if(abs(temp_itau-itaue).lt.1) temp_itau=temp_itau-3600 if(temp_itauitaue) then cycle ! exclude events not within this run endif if(temp_lat<-90.or.temp_lat>90) then cycle ! exclude events not within this run endif if(temp_lon<-180.or.temp_lon>180) then cycle ! exclude events not within this run endif !WP! Proceed to fill array concentration with obs within start and end of this run n_event=n_event+1 ! goUpCase is *probably* unnecessary concentration(n_event)%ident=goUpCase(trim(adjustl(temp_ident))) concentration(n_event)%lat=temp_lat concentration(n_event)%lon=temp_lon concentration(n_event)%itau=temp_itau concentration(n_event)%height=temp_height concentration(n_event)%idate=temp_idate concentration(n_event)%eventnumber=temp_eventnumber ! Note that the is and js components for the observation locations are ! set *before* the site moves are applied. Is this a bug? ! -Andy, 7 Nov 2006 concentration(n_event)%is = & int((concentration(n_event)%lon-float(xbeg(1)))/dxr + 0.99999) concentration(n_event)%js = & int((concentration(n_event)%lat-float(ybeg(1)))/dyr + 0.99999) concentration(n_event)%lsmin = 999 concentration(n_event)%lsmax = 0 print('(a24,1x,f8.2,f8.2,f8.1,1x,6i4,x,i30)'), & concentration(n_event)%ident, & concentration(n_event)%lat, & concentration(n_event)%lon, & concentration(n_event)%height, & concentration(n_event)%idate, & concentration(n_event)%eventnumber end do ! close stationfile close(sunit0) endif ! share the number of real events: n_event ! #ifdef MPI call MPI_BCAST(n_event ,1, MPI_INTEGER, root ,MPI_COMM_WORLD,ierr) #endif if(n_event==0) then ! no data in range write(gol,*) 'WARNING: No NOAA ESRL events within this run, skipping output ' ; call goPr deallocate(concentration) noaa_data=.false. status=0 return else write(gol,*) 'Found ',n_event,' NOAA ESRL events within this run' ; call goPr endif ! determine assignment of station to region ! (with highest refinement factor) #ifdef MPI do i_stat=1, n_event call MPI_BCAST(concentration(i_stat)%ident ,24, MPI_CHARACTER, & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%lat ,1, MY_REAL , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%lon ,1, MY_REAL , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%height,1, MY_REAL , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%idate,6, MPI_INTEGER , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%eventnumber,1, MPI_INTEGER , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%itau,1, MPI_INTEGER , & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%is ,1, MPI_INTEGER, & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%js ,1, MPI_INTEGER, & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%lsmin ,1, MPI_INTEGER, & root ,MPI_COMM_WORLD,ierr) call MPI_BCAST(concentration(i_stat)%lsmax ,1, MPI_INTEGER, & root ,MPI_COMM_WORLD,ierr) end do #endif do i_stat=1, n_event ! assume global region as default for average mixing ratios concentration(i_stat)%region = 1 concentration(i_stat)%is=1 concentration(i_stat)%js=1 concentration(i_stat)%n_accumulate = 0 ! no model values addded yet concentration(i_stat)%mmix(:) = 0.0 ! no model values addded yet do region=1, nregions dyr = dy/yref(region) dxr = dx/xref(region) if ( (concentration(i_stat)%lon .gt. xbeg(region) .and. & concentration(i_stat)%lon .lt. xend(region)) .and. & (concentration(i_stat)%lat .gt. ybeg(region) .and. & concentration(i_stat)%lat.lt.yend(region) ) ) then !===================== ! station is in region !===================== ! determine region with hightes refinement factor ! if (xref(region) > xref(concentration(i_stat)%region)) then concentration(i_stat)%region = region concentration(i_stat)%is = & int((concentration(i_stat)%lon-float(xbeg(region)))/dxr + 0.99999) concentration(i_stat)%js = & int((concentration(i_stat)%lat-float(ybeg(region)))/dyr + 0.99999) ! end if end if end do end do call Par_Barrier return ! error handling 1000 call escape_tm( 'read_noaa_events: could not open ' // trim(InputDir)//'/'//EventlistFilename) end subroutine init_noaa_events subroutine get_noaa(region,itau_f,force) ! ! This subroutine samples the model at given locations in array ! concentrations and puts the values in the tag: forecast ! This is done only in a certain time window, or when the keyword /force is ! used. Thus, a model field can be sampled at sites instantaneously (force), ! or at appropriate time ( omit force) ! ! Note: the values accumulate so care should be taken to clear the forecast ! and n_accumulate tags between functions ! ! use global_data, only: mass_dat, region_dat use Meteo, only: gph_dat, m_dat use ParTools, only: myid,ntracet_ar implicit none ! input/output integer,intent(in) :: region ! idate_f: date for which output required... integer,intent(in) :: itau_f logical,optional :: force ! local real,dimension(:,:,:), pointer :: m,gph real,dimension(:,:,:,:), pointer :: rm, rxm, rym, rzm real,dimension(0:lm(region)) :: height integer :: i,is,js,l,n,isn,jsn,ls,j,offsetj, lst, lstn, lmr, lsn real :: flon,flat,falt,ris,rjs,dxr,dyr,wcx,wcy,rls, wcz ! start m => m_dat(region)%data !pointers to global arrays... rm => mass_dat(region)%rm_t rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t gph => gph_dat(region)%data lmr=lm(region) fcast_loop: do n=1,n_event ! 0. Is idate equal to idate_site ! 1. Is the site in the area?--->no, then put -1 in c ! 2. Determine gridbox ! 3. Use slopes to determine concentration at the site. if ( (concentration(n)%itau.ge.itau_f-window/2.and.& concentration(n)%itau.lt.itau_f+window/2).or.force ) then dyr = dy/yref(region) dxr = dx/xref(region) flon=concentration(n)%lon flat=concentration(n)%lat falt=concentration(n)%height if (concentration(n)%region == region) then ris = (flon-float(xbeg(region)))/dxr + 0.99999 rjs = (flat-float(ybeg(region)))/dyr + 0.99999 is = int(ris) ! i-index of grid cell in which station is located js = int(rjs) ! j-index of grid cell in which station is located !fraction from the center of the is-box (-0.5---+0.5) ris = ris-is-0.5 !idem js rjs = rjs-js-0.5 !the neighbour for pressure interpolation if(ris .gt. 0) then isn = is+1 else isn = is-1 endif !the neighbour for y interpolation if(rjs .gt. 0) then jsn = js+1 else jsn = js-1 endif ! x- / y-weighting of grid cell in which station is located wcx = (1.0-abs(ris)) ! 1.0 ... 0.5 wcy = (1.0-abs(rjs)) ! 1.0 ... 0.5 !================================================================= ! if index of neighbour is exceeding range of region set ! neighbour = current cell (i.e. no interpolation) ! in case of cyclic x-boundaries take corresponding cyclic i index !================================================================= if ( jsn < 1) jsn=1 if ( jsn > jm(region) ) jsn=jm(region) ! isn-->jsn (wouter Peters) if ( xcyc(region) == 0 ) then ! non-cyclic boundaries if ( isn < 1) isn=1 if ( isn > im(region) ) isn=im(region) else ! cyclic x-boundaries if ( isn < 1 ) isn=im(region) if ( isn > im(region) ) isn=1 end if ! interpolate the altitude to site position... ls = 1 !layer do l=0,lm(region) height(l) = wcx * wcy* gph(is,js,l+1) + & (1.0-wcx)* wcy* gph(isn,js,l+1) + & wcx *(1.0-wcy)* gph(is,jsn,l+1) + & (1.0-wcx)*(1.0-wcy)* gph(isn,jsn,l+1) enddo do l=2,lm(region) ! selects layer , note that we start from second layer from surface if(height(l).gt.falt) exit enddo select case(l) case(0) if(myid==root) print*,'get_noaa: Warning..., forecast altitude ',& 'is below the surface height',height(0),' > ',falt,concentration(n)%ident ls = 1 rls = -0.5 !surface... case default ls = l !the site layer ! the off-set from the center of the layer (-0.5--->+0.5) ! (interpolation is in (m)) rls = (falt-height(l-1))/(height(l)-height(l-1)) - 0.5 end select !================================= !the neighbour for z interpolation !================================= if ( rls .gt. 0 ) then lsn = ls+1 else lsn = ls-1 end if ! z-weighting of grid cell in which station is located wcz = (1.0-abs(rls)) !.0 ... 0.5 !========================================================= ! if vertical neighbor is 0 (which does not exist) ! take vertical layer with l=2 for EXTRApolation to ground !========================================================= IF(lsn == 0) THEN lsn=2 wcz=1.0-rls ! 1.0 ... 1.5 ENDIF IF(lsn == lmr+1) THEN !========================================================= ! if vertical neighbor is lmr+1 (which does not exist) ! -> no interpolation !========================================================= lsn=lmr ! no interpolation wcz=1.0 ENDIF !from is,js,ls, ris,rjs,rls, determine the mixing ratio of each tracer ... do j=1,ntracetloc ! rm-value is obtained from rm + slopes. ! slope = rxm = (rm*dX/dx *deltaX/2) ! For EnKF forward, not sure that sum over a zero ! length subarray works properly. Let's be specific ! anyway. offsetj = sum(ntracet_ar(0:myid-1)) rmf = ( rm(is,js,ls,j) + & 2.0*(ris*rxm(is,js,ls,j) + & rjs*rym(is,js,ls,j) + & rls*rzm(is,js,ls,j) ) )/m(is,js,ls)*fscale(offsetj+j) concentration(n)%mmix(offsetj+j)=concentration(n)%mmix(offsetj+j)+rmf !*fscale(offsetj+j) enddo concentration(n)%n_accumulate=concentration(n)%n_accumulate+1 !! if(myid==root) then ! print*,'DEBUG Check ',myid, concentration(n)%ident,concentration(n)%n_accumulate,concentration(n)%model(:)/concentration(n)%n_accumulate ! endif endif endif ! force/time enddo fcast_loop nullify(m) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(gph) end subroutine get_noaa subroutine evaluate_noaa ! convert concentrations to ppm and average over n_accumulate integer :: n do n=1,ntracet where(concentration%n_accumulate.ge.1) & concentration%mmix(n)=1.e6*concentration%mmix(n)/concentration%n_accumulate enddo end subroutine evaluate_noaa end module user_output_noaa