! ********************* DUMMY for TM5-MP ************************** !################################################################# !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 ! real :: weight ! accumulated time weighting, replaces n_accumulate in evaluation ! 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' STATUS=0 ! 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 STATUS=0 ! 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)%weight = 0.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 ! use dims, only : ndyn,ndyn_max ! implicit none ! ! input/output integer,intent(in) :: region ! ! idate_f: date for which output required... integer(kind=8),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 ! real :: this_weight ! ! 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) ! this_weight=real(ndyn)/real(ndyn_max) ! 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 ! concentration(n)%weight=concentration(n)%weight + this_weight ! !! 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%weight ! enddo ! end subroutine evaluate_noaa end module user_output_noaa