!################################################################# ! ! MODULE: user_output_column *************** DUMMY TM5-MP VERSION ******************* ! ! PUBLIC SUBROUTINES: ! user_output_column_init ! user_output_column_accum ! user_output_column_evaluate ! user_output_column_write ! user_output_column_reset ! user_output_column_done ! ! DESCRIPTION: ! Write column tracer fields. Called from user_output.F90. ! ! REVISION HISTORY: ! Original module written by: ! Peter Bergamaschi, Maarten Krol, Wouter Peters, others? ! ! Mike Trudeau, Feb 2012 ! Modified to write NetCDF4 output in NOAA CarbonTracker ! release format. ! ! Required rc file keys ! output.column : [T/F] ! output.column.dtsec : [integer] ! output.column.infile : [/path/myColumnList.txt] ! Obsolete rc file keys ! output.column.filename : [column.nc] ! inputdir.column : [/path] ! output.column.infilename : [myColumnList.txt] ! ! Mike Trudeau, Feb 2013 ! 1. Modified pressure units (hectopascals to pascals) ! 2. Removed mean pressure, wind speed, wind direction from ! the meteo output fields. ! 3. CO2 tracers written with background subtracted. ! 4. Total CO2 tracer added to output. ! ! Andy Jacobson, April 2013 ! Change to ndyn weighting ! !### 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_column use GO, only : gol, goErr, goPr ! use dims, only : nregions, itaur ! use dims, only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend, zbeg, zend ! use chem_param, only : ntracet, ntrace, names implicit none character(len=*), parameter :: mname = 'user_output_column' private ! --- public routines ------------------------ public :: user_output_column_init public :: user_output_column_done public :: user_output_column_write public :: user_output_column_accum public :: user_output_column_reset public :: user_output_column_evaluate ! ! --- constants ------------------------------ ! integer, parameter :: nmethods = 3 ! number of column interpolation methods ! integer, parameter :: nmeteoLvl = 4 ! integer, parameter :: nmeteoBnd = 2 ! integer, parameter :: writeMethod = 2 ! interpolation method output to file ! ! 1: no interpolation ('grd') ! ! 2: interpolation based on slopes ('slp') ! ! 3: 3D linear interpolations ('int') ! character(len=70), dimension(5), parameter :: comment_species = & ! (/'background - component due to initial condition (at 2000-01-01)', & ! 'component due to fossil fuel burning (since 2000-01-01)', & ! 'component due to air-sea gas exchange (since 2000-01-01)', & ! 'component due to terrestrial biosphere exchange (since 2000-01-01)', & ! 'component due to direct fire emissions (since 2000-01-01)'/) ! character(len=14), dimension(nmeteoLvl), parameter :: name_meteoLvl = & ! (/'temperature', & ! 'u', & ! 'v', & ! 'blh'/) ! character(len=13), dimension(nmeteoLvl), parameter :: units_meteoLvl = & ! (/'Kelvin', & ! 'meters/second', & ! 'meters/second', & ! 'meters'/) ! character(len=50), dimension(nmeteoLvl), parameter :: standard_name_meteoLvl = & ! (/'temperature at level center', & ! 'zonal component of the wind velocity', & ! 'meridional component of the wind velocity', & ! 'boundary layer height'/) ! character(len=14), dimension(nmeteoLvl), parameter :: long_name_meteoLvl = & ! (/'temperature', & ! 'u', & ! 'v', & ! 'blh'/) ! character(len=14), dimension(nmeteoBnd), parameter :: name_meteoBnd = & ! (/'pressure', 'gph'/) ! character(len=13), dimension(nmeteoBnd), parameter :: units_meteoBnd = & ! (/'pascals', 'meters'/) ! character(len=50), dimension(nmeteoBnd), parameter :: standard_name_meteoBnd = & ! (/'pressure at layer interface', 'geopotential height of layer interface'/) ! character(len=14), dimension(nmeteoBnd), parameter :: long_name_meteoBnd = & ! (/'pressure', 'gph'/) ! ! --- variables ------------------------------ ! integer, dimension(ntrace) :: species ! integer :: nspecies ! integer, save :: nstations ! logical :: evaluated ! type fieldsType1 ! character(len=13) :: id ! column identifier ! character(len=50) :: desc ! column name ! real :: lat ! latitude of column ! real :: lon ! longitude of column ! real :: height ! height of column ! real :: oro ! ground height at station ! integer :: index ! index of column ! integer :: region ! integer :: is ! integer :: js ! integer :: ls ! real,dimension(nmeteoLvl,lm(1)) :: meteoLvl ! meteo (interpolated) at columns ! real,dimension(:,:,:),pointer :: c_avg ! species, interpolation ! real,dimension(:,:,:),pointer :: c_std ! integer :: ndata ! for averages ! real :: accum_weight ! for averages ! real,dimension(nmeteoBnd,lm(1)+1) :: meteoBnd ! meteo (interpolated) at columns ! end type fieldsType1 ! type(fieldsType1), dimension(:), allocatable :: stations ! type(fieldsType1), dimension(:), allocatable :: stationsRegion ! type fieldsType2 ! integer, dimension(6) :: idate_start ! integer, dimension(6) :: idate_end ! end type fieldsType2 ! type(fieldsType2), dimension(nregions), save :: timestamp ! ! --- begin ------------------------------- contains subroutine user_output_column_init( status ) ! use ParTools , only : myid, root ! use dims, only : nregions, idatei ! implicit none ! ! --- in/out ------------------------ ! integer :: region integer, intent(out) :: status ! --- constants --------------------- ! character(len=*), parameter :: rname = mname//'/user_output_column_init' ! ! --- begin ------------------------------- ! if (myid /= root) return ! evaluated = .FALSE. ! ! grab the date ! do region = 1, nregions ! timestamp(region)%idate_start = idatei ! end do ! ! read user supplied station list ! call read_columnlist( status ) ! IF_NOTOK_RETURN( status=1 ) ! ! initialize the column output ! call user_output_column_reset( status ) ! ok status = 0 end subroutine user_output_column_init ! subroutine read_columnlist( status ) ! use toolbox, only : escape_tm ! #ifdef MPI ! use mpi_const, only: myid, root, mpi_integer, mpi_character, & ! my_real, mpi_comm_world, ierr ! #else ! use ParTools, only: myid, root ! #endif ! use GO, only : TrcFile, Init, Done, ReadRc ! use global_data, only : rcfile ! implicit none ! ! --- variables --------------------- ! integer, intent(out) :: status ! character(len=1024) :: infile ! integer :: nchar !string length ! character(len=200) :: dummy ! integer :: q ! integer :: region ! integer,dimension(6) :: idate_sim ! integer :: itau_sim, isim, ispec ! real :: dxr, dyr ! character(len=8) :: name_spec !CMK changed from 6-->8 ! integer :: sunit0=110 ! type(TrcFile) :: rcF ! ! --- constants ----------------------- ! character(len=*), parameter :: rname = mname//'/read_columnlist' ! ! --- begin ----------------------- ! ! start ! if ( myid == root ) then ! ! read rc file ! call Init( rcF, rcfile, status) ! call ReadRc( rcF, 'output.column.infile', infile, status ) ! IF_NOTOK_RETURN( status=1 ) ! write (gol, *) trim(rname)//'/column file: ', trim(infile); call goPr ! call Done( rcF, status) ! ! read columnlist file ! nchar = len_trim(infile) ! if ( nchar > 0 ) then ! !========================================= ! ! open columnFile (with list of columns) ! !========================================= ! open(UNIT=sunit0, FORM='FORMATTED', & ! FILE=trim(infile), & ! status='OLD', ERR=1000) ! read(sunit0, '(a)') dummy ! ! count number of columns ! nstations=0 ! stations_loop: DO ! read(sunit0, '(a)', END=100) dummy ! if (dummy(1:7) /= 'species') then ! nstations = nstations+1 ! else ! nspecies = 0 ! do ! read(sunit0, '(a8)' , end=100) name_spec ! CMK changed a6-->a8 ! !CMK not allowed yet to write short-lived: TODO ! specloop: do ispec = 1, ntracet ! if (name_spec == names(ispec) ) then ! nspecies = nspecies + 1 ! species(nspecies) = ispec ! exit specloop ! endif ! enddo specloop ! enddo ! exit stations_loop ! end if ! end do stations_loop ! 100 continue ! print *, '_________________________________________________' & ! //'_______________________________' ! print *, nstations, ' columns read from file: ', infile ! print *, nspecies, ' concentrations will be gathered:' ! if (nspecies == 0) then ! print *, 'ERROR user_output_column' ! print *, 'ERROR no species for output' ! print *, 'ERROR Set column_output to .F. or give valid species name' ! call escape_tm('Error reading column file') ! endif ! do ispec = 1, nspecies ! print '(2i4,2x,a8)', ispec, species(ispec), & ! CMK change a6--->a8 ! names(species(ispec)) ! enddo ! print *, '_________________________________________________' & ! //'_______________________________' ! allocate(stations(nstations)) ! do q = 1, nstations ! ! wow, this lm(1) assumes column in region 1, or that ! ! child regions all have lm(n) <= lm(1) ! allocate(stations(q)%c_avg(nspecies,nmethods,lm(1))) ! allocate(stations(q)%c_std(nspecies,nmethods,lm(1))) ! enddo ! dyr = dy/yref(1) ! again, as above comment: region 1 assumed. ! dxr = dx/xref(1) ! " " " ! !=============== ! ! read columns ! !=============== ! rewind(sunit0) ! read(sunit0, '(a)') dummy ! do q=1, nstations ! read(sunit0, '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') & ! stations(q)%id, & ! stations(q)%lat, & ! stations(q)%lon, & ! stations(q)%height, & ! stations(q)%desc ! stations(q)%oro = -99999.0 ! stations(q)%index = q ! ! assume global region as default for average mixing ratios ! stations(q)%region = 1 ! stations(q)%is = & ! int((stations(q)%lon-float(xbeg(1)))/dxr + 0.99999) ! stations(q)%js = & ! int((stations(q)%lat-float(ybeg(1)))/dyr + 0.99999) ! end do ! ! close columnfile ! close(sunit0) ! do q=1, nstations ! do region=1, nregions ! dyr = dy/yref(region) ! dxr = dx/xref(region) ! if ( (stations(q)%lon .ge. xbeg(region) .and. & ! stations(q)%lon .le. xend(region)) .and. & ! (stations(q)%lat .ge. ybeg(region) .and. & ! stations(q)%lat .le. yend(region) ) ) then ! !===================== ! ! column is in region ! !===================== ! ! determine region with hightes refinement factor ! if (xref(region) > xref(stations(q)%region)) then ! stations(q)%region = region ! stations(q)%is = & ! int((stations(q)%lon-float(xbeg(region)))/dxr + 0.99999) ! stations(q)%js = & ! int((stations(q)%lat-float(ybeg(region)))/dyr + 0.99999) ! end if ! end if ! end do ! end do ! do q=1, nstations ! print '(a7,1x,f8.2,f8.2,f8.1,1x,i4,1x,a50)', & ! stations(q)%id, & ! stations(q)%lat, & ! stations(q)%lon, & ! stations(q)%height, & ! stations(q)%region, & ! stations(q)%desc ! end do ! end if ! end if ! myid = root ! ! ok ! status = 0 ! return ! ! error handling ! 1000 call escape_tm( 'read_columnlist: could not open '//trim(infile)) ! end subroutine read_columnlist subroutine user_output_column_accum( region, status ) ! #ifdef MPI ! use mpi_const !, only: myid, root ! use mpi_comm, only: gather_tracer_k, gather_tracer_t, gather_tracer_k_3d ! #endif ! use ParTools, only: myid, root, root_t, root_k ! use global_data, only : mass_dat, conv_dat ! use datetime, only : tau2date ! use toolbox, only : escape_tm ! use Meteo, only : gph_dat, oro_dat, temper_dat, pu_dat, pv_dat, m_dat, phlb_dat ! use dims, only : xcyc,ndyn,ndyn_max ! use chem_param, only : fscale, uscale ! use binas, only : ae, twopi ! implicit none ! --- in/out ------------------------ integer, intent(out) :: status ! --- variables ----------------------- integer, intent(in) :: region ! real,dimension(:,:,:,:),pointer :: rm ! tracer mass ! real,dimension(:,:,:,:),pointer :: rxm,rym,rzm ! slopes of tracer mass ! real,dimension(:,:,:) ,pointer :: m ! air mass ! real,dimension(:,:,:) ,pointer :: gph ! geopotential height ! real,dimension(:,:) ,pointer :: oro ! orography ! real, dimension(:,:,:,:), allocatable, target :: rmg ! MPI arrays to gather fields. ! real, dimension(:,:,:,:), allocatable, target :: rxmg, rymg, rzmg ! real, dimension(:,:,: ), allocatable, target :: mg ! integer :: i_interpol ! interpolation procedures: ! ! 1: no interpolation ! ! 2: interpolation based on slopes ! ! 3: 3D linear interpolations ! integer :: i,j,l ! grid cell indices ! integer :: is,js,ls ! i,j,l index of grid cell in which column lis located ! integer :: isn,jsn,lsn ! i,j,l index of grid cell which is taken as neighbour for linear interpolation ! integer :: lst, lstn ! MPI implementation ! integer :: js_north, js_south ! j index of grid cell for neighbours for windspeed_v interpolation ! integer :: is_west ! i index of grid cell for neighbour for windspeed_u interpolation ! integer :: n ! index of tracer ! integer :: q ! index of column ! real :: ris,rjs,rls ! deviation of column from center of the grid cell (-0.5 ... +0.5) ! real :: wcx,wcy,wcz ! x,y,z-weight of grid cell (1.0 ... 0.5) ! real :: dxr, dyr ! x,y resolution (in degrees) for current region ! real,dimension(0:lm(region)) :: height ! height of vertical grid boundaries ! integer, dimension(6) :: idate_f ! current idate for region ! integer :: sunit ! unit number for column output file ! real,dimension(lm(1)) :: rm_interpol ! tracer mass, interpolated to column location ! integer :: ispec ! index over activated tracers ! !============= ! ! meteoLvl output ! !============= ! real,dimension(nstations,nmeteoLvl,lm(1)) :: meteoLvl ! meteo (interpolated) at columns ! real,dimension(:,:,:), pointer :: T ! temperature ! real,dimension(:,:,:), pointer :: phlb ! pressure grid boundaries ! real,dimension(:,:,:), pointer :: pu ! mass flux x-direction [kg/s] ! real,dimension(:,:,:), pointer :: pv ! mass flux y-direction [kg/s] ! real,dimension(:,:), pointer :: blh ! boundary layer height [m] ! real,dimension(lm(1)) :: T_interpol ! temperature, interpolated to column location ! real,dimension(lm(1)) :: p_interpol ! pressure, interpolated to column location ! real,dimension(lm(1)+1) :: p_interpol_bnd ! boundary pressure, interpolated to column location ! real,dimension(lm(1)) :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to column location ! real,dimension(lm(1)) :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to column location ! real,dimension(lm(1)) :: windspeed_interpol ! wind speed [m/s] , interpolated to column location ! real,dimension(lm(1)) :: winddir_interpol ! wind direction , interpolated to column location ! real :: blh_interpol ! boundary layer height , interpolated to column location ! real :: lat_j ! latitude of grid cell center js [degrees] ! real :: lat_jn ! latitude of grid cell center jsn [degrees] ! real :: ddx_j ! x-extension of grid cell js [m] ! real :: ddx_jn ! x-extension of grid cell jsn [m] ! real :: ddy ! y-extension of grid cell js [m] ! real,dimension(lm(1)) :: vmr ! volume mi ! real :: this_weight ! integer :: communicator,root_id,lmr,imr,jmr ! ! start ! ! assign pointers depending on the paralel status ! ! (over levels after chemistry!) ! imr = im(region) ; jmr = jm(region) ; lmr = lm(region) ! this_weight=real(ndyn)/real(ndyn_max) ! #ifdef MPI ! if ( root_k /= root_t ) print *, 'Warning from MPI column output: ', & ! 'root_k and root_t are not equal' ! allocate( rmg (-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rmg = 0.0 ! allocate( rxmg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rxmg = 0.0 ! allocate( rymg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rymg = 0.0 ! allocate( rzmg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rzmg = 0.0 ! allocate( mg (-1:imr+2,-1:jmr+2,lmr ) ) ; mg = 0.0 ! if ( which_par == 'tracer' ) then ! m => m_dat(region)%data ! rm => mass_dat(region)%rm_t ! rxm => mass_dat(region)%rxm_t ! rym => mass_dat(region)%rym_t ! rzm => mass_dat(region)%rzm_t ! call gather_tracer_t(rmg,imr,jmr,lmr,2,2,0,ntracet,rm,.false.) ! call gather_tracer_t(rxmg,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.) ! call gather_tracer_t(rymg,imr,jmr,lmr,2,2,0,ntracet,rym,.false.) ! call gather_tracer_t(rzmg,imr,jmr,lmr,2,2,0,ntracet,rzm,.false.) ! nullify(rm) ! reassign pointers ! nullify(rxm) ! nullify(rym) ! nullify(rzm) ! rm => rmg ! rxm => rxmg ! rym => rymg ! rzm => rzmg ! communicator=com_trac !WP! assign com_trac as communicator ! root_id=root_t !root_t is the PE that knows the result. ! else ! m => m_dat(region)%data ! rm => mass_dat(region)%rm_k ! rxm => mass_dat(region)%rxm_k ! rym => mass_dat(region)%rym_k ! rzm => mass_dat(region)%rzm_k ! call gather_tracer_k(rmg,imr,jmr,lmr,2,2,0,ntracet,rm,.false.) ! call gather_tracer_k(rxmg,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.) ! call gather_tracer_k(rymg,imr,jmr,lmr,2,2,0,ntracet,rym,.false.) ! call gather_tracer_k(rzmg,imr,jmr,lmr,2,2,0,ntracet,rzm,.false.) ! call gather_tracer_k_3d( mg, imr,jmr,lmr, 2,2,0, m, .false.) ! nullify(m) ! reassign pointers ! nullify(rm) ! nullify(rxm) ! nullify(rym) ! nullify(rzm) ! m => mg ! rm => rmg ! rxm => rxmg ! rym => rymg ! rzm => rzmg ! communicator=com_lev !WP! assign com_lev as communicator ! root_id=root_k ! end if ! #else ! m => m_dat(region)%data ! rm => mass_dat(region)%rm_t ! rxm => mass_dat(region)%rxm_t ! rym => mass_dat(region)%rym_t ! rzm => mass_dat(region)%rzm_t ! #endif ! ! assign pointers (METEO) ! gph => gph_dat(region)%data ! oro => oro_dat(region)%data(:,:,1) ! t => temper_dat(region)%data ! pu => pu_dat(region)%data ! pv => pv_dat(region)%data ! blh => conv_dat(region)%blh ! ! point to pressure paralel over tracers (contains all layers!) ! phlb => phlb_dat(region)%data ! if ( myid == root_t ) then ! ! perform only on root: tracers are already gathered here ! ! x and y resolution [degrees] for current region ! dyr = dy/yref(region) ! dxr = dx/xref(region) ! ! idate for current region ! call tau2date(itaur(region),idate_f) ! !==================== ! ! loop over columns ! !==================== ! do q=1, nstations ! IF (stations(q)%region == region)THEN ! !avoid edges! ! ris = (stations(q)%lon-float(xbeg(region)))/dxr + 0.99999 ! rjs = (stations(q)%lat-float(ybeg(region)))/dyr + 0.99999 ! is = int(ris) ! i-index of grid cell in which column is located ! js = int(rjs) ! j-index of grid cell in which column is located ! ! x-deviation from the center of the grid cell (-0.5...+0.5) ! ris = ris-is-0.5 ! ! y-deviation from the center of the grid cell (-0.5...+0.5) ! rjs = rjs-js-0.5 ! !================================= ! !the neighbour for x 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 column 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 == 0) jsn=1 ! if ( jsn == jm(region)+1 ) jsn=jm(region) ! isn-->jsn (wouter Peters) ! if ( xcyc(region) == 0 ) then ! ! non-cyclic boundaries ! if ( isn == 0) isn=1 ! if ( isn == im(region)+1 ) isn=im(region) ! else ! ! cyclic x-boundaries ! if ( isn == 0 ) isn=im(region) ! if ( isn == im(region)+1 ) isn=1 ! end if ! !============================================ ! ! interpolate the gph to xy position of column ! !============================================ ! !ls = 1 !layer ! do l=0,lm(region) ! !CMK note: gph now from 1-->lm+1 (dec 2002) ! 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) ! end do ! !=========================== ! ! determine layer of column ! !=========================== ! !do l=0,lm(region) ! ! if(height(l) .gt. stations(q)%height) exit ! !end do ! !select case(l) ! !case(0) ! ! ! force column to model surface ! ! ls = 1 ! ! rls = -0.5 ! !case default ! ! ls = l ! column's layer ! ! ! the off-set from the center of the layer (-0.5--->+0.5) ! ! ! (interpolation is in (m)) ! ! rls = (stations(q)%height - height(l-1)) / & ! ! (height(l)-height(l-1)) - 0.5 ! !end select ! !stations(q)%ls = ls ! store for output to file ! rls = 0.0 ! no interpolation ! wcz = 1.0 ! no interpolation ! !============================================ ! ! if not yet done, interpolate the orography ! ! to xy position of column ! !============================================ ! if (stations(q)%oro < -10000) then ! stations(q)%oro = wcx * wcy* oro(is,js) + & ! (1.0-wcx)* wcy* oro(isn,js) + & ! wcx *(1.0-wcy)* oro(is,jsn) + & ! (1.0-wcx)*(1.0-wcy)* oro(isn,jsn) ! endif ! do ispec=1,nspecies ! n = species(ispec) ! lst = ls ! lstn = lsn ! !======================================== ! ! value of grid box without interpolation ! !======================================== ! vmr = rm(is,js,:,n) / m(is,js,:) * fscale(n) * uscale(n) ! stations(q)%c_avg(ispec,1,:) = & ! stations(q)%c_avg(ispec,1,:) + this_weight*vmr ! stations(q)%c_std(ispec,1,:) = & ! stations(q)%c_std(ispec,1,:) + this_weight*vmr*vmr ! !======================================== ! ! interpolation using slopes ! !======================================== ! !rm-value is obtained from rm + slopes. ! !slope = rxm = (rm*dX/dx *deltaX/2) ! rm_interpol = rm(is,js,:,n) + 2.0*(ris*rxm(is,js,:,n) + & ! rjs*rym(is,js,:,n) + & ! rls*rzm(is,js,:,n) ) ! rm_interpol = max(0.0,rm_interpol) ! vmr = rm_interpol/m(is,js,:) * fscale(n) * uscale(n) ! stations(q)%c_avg(ispec,2,:) = & ! stations(q)%c_avg(ispec,2,:) + this_weight*vmr ! stations(q)%c_std(ispec,2,:) = & ! stations(q)%c_std(ispec,2,:) + this_weight*vmr*vmr ! !======================================== ! ! 3D linear interpolation ! !======================================== ! ! this PE handles also the neighbour layer ! !PB if ( lstn >= 1 .and. lstn <= lmr) then ! !PBi ! rm_interpol = wcx *wcy *rm(is,js,:,n) /m(is,js,:) + & ! (1.0-wcx)*wcy *rm(isn,js,:,n) /m(isn,js,:) + & ! wcx *(1.0-wcy) *rm(is,jsn,:,n) /m(is,jsn,:) + & ! (1.0-wcx)*(1.0-wcy) *rm(isn,jsn,:,n) /m(isn,jsn,:) ! vmr = rm_interpol * fscale(n) * uscale(n) ! stations(q)%c_avg(ispec,3,:) = & ! stations(q)%c_avg(ispec,3,:) + this_weight*vmr ! stations(q)%c_std(ispec,3,:) = & ! stations(q)%c_std(ispec,3,:) + this_weight*vmr*vmr ! end do !ispec ! !====================== ! !meteoLvl data at columns ! !====================== ! !if (output_columns_meteoLvl) then ! !============================================ ! ! gph already interpolated to xy position of column ! ! Note that height is defined with index 0:lm(region) ! ! whereas meteoLvl(8,:) is 1:lm(region) ! !============================================ ! stations(q)%meteoBnd(2,:) = & ! stations(q)%meteoBnd(2,:) + this_weight*height(0:(lm(region))) ! !============================================ ! ! interpolate the blh to xy position of column ! !============================================ ! blh_interpol = wcx * wcy* blh(is,js) + & ! (1.0-wcx)* wcy* blh(isn,js) + & ! wcx *(1.0-wcy)* blh(is,jsn) + & ! (1.0-wcx)*(1.0-wcy)* blh(isn,jsn) ! stations(q)%meteoLvl(4,:) = & ! stations(q)%meteoLvl(4,:) + this_weight*blh_interpol ! !==================================== ! ! 3D linear interpolation Temperature ! !==================================== ! T_interpol = & ! wcx *wcy *wcz *T(is,js,:) + & ! (1.0-wcx)*wcy *wcz *T(isn,js,:) + & ! wcx *(1.0-wcy)*wcz *T(is,jsn,:) + & ! (1.0-wcx)*(1.0-wcy)*wcz *T(isn,jsn,:) ! stations(q)%meteoLvl(1,:) = & ! stations(q)%meteoLvl(1,:) + this_weight*T_interpol ! !==================================== ! ! 3D linear interpolation pressure ! !==================================== ! p_interpol =(((0.5-rls) * phlb(is,js,1:lm(1)) + (0.5+rls) * phlb(is,js,2:lm(1)+1)) * wcx * wcy + & ! ((0.5-rls) * phlb(isn,js,1:lm(1)) + (0.5+rls) * phlb(isn,js,2:lm(1)+1)) * (1.0-wcx) * wcy + & ! ((0.5-rls) * phlb(is,jsn,1:lm(1)) + (0.5+rls) * phlb(is,jsn,2:lm(1)+1)) * wcx * (1.0-wcy) + & ! ((0.5-rls) * phlb(isn,jsn,1:lm(1))+ (0.5+rls) * phlb(isn,jsn,2:lm(1)+1))* (1.0-wcx) * (1.0-wcy) ) ![Pa] ! !stations(q)%meteoLvl(2,:) = & ! ! stations(q)%meteoLvl(2,:) + p_interpol ! p_interpol_bnd =(phlb(is,js,1:lm(1)+1) * wcx * wcy + & ! phlb(isn,js,1:lm(1)+1) * (1.0-wcx) * wcy + & ! phlb(is,jsn,1:lm(1)+1) * wcx * (1.0-wcy) + & ! phlb(isn,jsn,1:lm(1)+1) * (1.0-wcx) * (1.0-wcy)) ![Pa] ! stations(q)%meteoBnd(1,:) = & ! stations(q)%meteoBnd(1,:) + this_weight*p_interpol_bnd ! !==================================== ! ! 3D windspeed_u_interpol (x-direction) ! !==================================== ! ! pu is eastward flux through east grid box boundary ! ! windspeed_u_interpol wind component from east ! ! latitude of grid cell center js [degrees] ! lat_j =ybeg(region)+(js -0.5)*dyr ! ! latitude of grid cell center jsn [degrees] ! lat_jn=ybeg(region)+(jsn-0.5)*dyr ! ! x-extension of grid cell js [m] ! ddx_j =ae * twopi * dxr / 360.0 * cos(lat_j*twopi/360.0) ! ! x-extension of grid cell jsn [m] ! ddx_jn=ae * twopi * dxr / 360.0 * cos(lat_jn*twopi/360.0) ! is_west=is-1 ! if ( xcyc(region) == 0 ) then ! ! non-cyclic boundaries ! if ( is_west == 0 ) is_west=1 ! else ! ! cyclic x-boundaries ! if ( is_west == 0 ) is_west=im(region) ! end if ! !west border !east border ! windspeed_u_interpol=-(((0.5-ris) * pu(is_west,js,:)/m(is_west,js,:) + & ! (0.5+ris) * pu(is,js,:)/ m(is,js,:) ) * & ! ddx_j * wcy * wcz + & ! ((0.5-ris) * pu(is_west,jsn,:)/m(is_west,jsn,:) + & ! (0.5+ris) * pu(is,jsn,:)/m(is,jsn,:) ) * & ! ddx_jn * (1.0-wcy)* wcz + & ! ((0.5-ris) * pu(is_west,js,:)/m(is_west,js,:) + & ! (0.5+ris) * pu(is,js,:)/ m(is,js,:)) * & ! ddx_j * wcy * (1.0-wcz) + & ! ((0.5-ris) * pu(is_west,jsn,:)/m(is_west,js,:) + & ! (0.5+ris) * pu(is,jsn,:)/ m(is,jsn,:) ) * & ! ddx_jn * (1.0-wcy)* (1.0-wcz) ) ! !==================================== ! ! 3D windspeed_v_interpol (y-direction) ! !==================================== ! ! pv northward flux through south grid box boundary ! ! windspeed_v_interpol wind component from north ! ddy =ae * twopi * dyr / 360.0 ! y-extension of grid cell [m] ! js_south = js-1 ! if (js_south==0) js_south=1 ! js_north = js+1 ! if (js_north==(jm(region)+1)) js_north=jm(region) ! !south border !north border ! windspeed_v_interpol=-(( & ! (0.5-rjs)*pv(is,js,:)/m(is,js_south,:) + & ! (0.5+rjs)*pv(is,js_north,:)/ m(is,js,:) ) * & ! wcx * wcz + & ! ((0.5-rjs)*pv(isn,js,:)/m(isn,js_south,:) + & ! (0.5+rjs)*pv(isn,js_north,:)/ m(isn,js,:) ) * & ! (1.0 - wcx)* wcz + & ! ((0.5-rjs)*pv(is,js,:)/m(is,js_south,:) + & ! (0.5+rjs)*pv(is,js_north,:)/ m(is,js,:) ) * & ! wcx * (1.0-wcz) + & ! ((0.5-rjs)*pv(isn,js,:)/m(isn,js_south,:) + & ! (0.5+rjs)*pv(isn,js_north,:)/ m(isn,js,:) ) * & ! (1.0 - wcx)* (1.0-wcz) ) * & ! ddy ! !windspeed_interpol = & ! ! sqrt(windspeed_u_interpol**2 + windspeed_v_interpol**2) ! stations(q)%meteoLvl(2,:) = & ! stations(q)%meteoLvl(2,:) + this_weight*windspeed_u_interpol ! stations(q)%meteoLvl(3,:) = & ! stations(q)%meteoLvl(3,:) + this_weight*windspeed_v_interpol ! !endif ! meteoLvl out... ! stations(q)%ndata = stations(q)%ndata + 1 ! stations(q)%accum_weight = stations(q)%accum_weight + this_weight ! end if ! end do ! end q loop ! endif ! myid = root_t ! nullify(m) ! nullify(rm) ! nullify(rxm) ! nullify(rym) ! nullify(rzm) ! nullify(gph) ! nullify(oro) ! nullify(t) ! nullify(pu) ! nullify(pv) ! nullify(phlb) ! nullify(blh) ! #ifdef MPI ! deallocate(rmg) ! deallocate(rxmg) ! deallocate(rymg) ! deallocate(rzmg) ! deallocate(mg) ! #endif ! ok status = 0 end subroutine user_output_column_accum subroutine user_output_column_evaluate(status) !================================================================== ! evaluate mean concentration and std !================================================================== ! use dims, only : idate ! use ParTools, only: myid, root ! use toolbox, only : escape_tm ! use binas, only : twopi ! implicit none integer, intent(out) :: status ! integer q, i, ii,iii ! integer ndata, isds, nsds, itp ! real :: accum_weight ! real temp ! real, dimension(:,:), allocatable :: cout ! to get output ! integer, dimension(:,:), allocatable :: iidate ! real :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to column location ! real :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to column location ! real :: windspeed_interpol ! wind speed [m/s] , interpolated to column location ! real :: winddir_interpol ! wind direction , interpolated to column location ! if(myid /= root) return ! if(evaluated) return ! do q=1, nstations ! accum_weight = stations(q)%accum_weight ! if (accum_weight > 0.0) then ! stations(q)%c_avg(:,:,:) = stations(q)%c_avg(:,:,:) / accum_weight ! stations(q)%meteoLvl(:,:) = stations(q)%meteoLvl(:,:) / accum_weight ! stations(q)%meteoBnd(:,:) = stations(q)%meteoBnd(:,:) / accum_weight ! if (ndata > 1) then ! do i=1, nspecies ! do ii= 1, nmethods ! do iii= 1, lm(1) ! temp = (((stations(q)%c_std(i,ii,iii)/accum_weight) - real(ndata) * ((stations(q)%c_avg(i,ii,iii))**2) ) / (real(ndata)-1.0) ) ! if(temp > 0) then ! stations(q)%c_std(i,ii,iii) = sqrt(temp) ! else ! stations(q)%c_std(i,ii,iii) = 0.0 ! endif ! enddo ! enddo ! enddo ! else ! stations(q)%c_std(:,:,:)=-1.0 ! end if ! else ! stations(q)%c_avg(:,:,:)= -1.0 ! stations(q)%c_std(:,:,:) = -1.0 ! stations(q)%meteoLvl(:,:) = -999. ! stations(q)%meteoBnd(:,:) = -999. ! end if ! !do i=1,lm(1) ! ! windspeed_u_interpol = stations(q)%meteoLvl(3,i) ! ! windspeed_v_interpol = stations(q)%meteoLvl(4,i) ! ! windspeed_interpol = stations(q)%meteoLvl(5,i) ! ! !==================================== ! ! ! wind direction ! ! !==================================== ! ! ! definition of winddirection: ! ! ! east : 90 ! ! ! south : 180 ! ! ! west : 270 ! ! ! north : 360 ! ! if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol > 0))then ! wind from north east (0...90) ! ! winddir_interpol = atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi ! ! else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol < 0)) then ! wind from south east (90...180) ! ! winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi ! ! else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol == 0)) then ! wind from east (90) ! ! winddir_interpol = 90.0 ! ! else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol > 0)) then ! wind from north west (270... 360) ! ! winddir_interpol = 360.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi ! ! else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol < 0)) then ! wind from south west (180...270) ! ! winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi ! ! else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol == 0)) then ! wind from west (270) ! ! winddir_interpol = 270.0 ! ! else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol > 0)) then ! wind from north (360) ! ! winddir_interpol = 360.0 ! ! else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol < 0)) then ! wind from south (180) ! ! winddir_interpol = 180.0 ! ! else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol == 0)) then ! no wind ! ! winddir_interpol = -1.0 ! ! else ! ! call escape_tm('output_columnconc: error wind direction') ! ! end if ! ! stations(q)%meteoLvl(6,i) = winddir_interpol ! ! end do ! end do ! evaluated = .TRUE. ! ok status = 0 end subroutine user_output_column_evaluate subroutine user_output_column_write( region, status ) ! use ParTools, only : myid, root ! use global_data, only : outdir ! use dims, only : region_name ! use datetime, only : tau2date, date2tau, idate2ddate ! use MDF, only : MDF_Create, MDF_Close, MDF_EndDef ! use MDF, only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_DOUBLE, MDF_CHAR ! use MDF, only : MDF_Put_Att, MDF_Def_Dim, MDF_Def_Var, MDF_Put_Var ! implicit none ! --- in/out ------------------------ integer, intent(out) :: status integer, intent(in) :: region ! --- variables ----------------------- ! character(len=1024) :: outfile ! integer :: fid ! integer :: i, ii, q ! integer :: dimid_stations, dimid_date, dimid_cal, dimid_nchar ! integer :: dimid_lvl, dimid_bnd ! integer :: varid_desc, varid_id, varid_date_int, varid_column_int ! integer :: varid_elapsed, varid_dec, varid_cal_int, varid_oro ! integer :: varid_lvl, varid_bnd, varid_lon, varid_lat, varid_hgt, varid_total ! integer, dimension(nspecies, nmethods) :: varid_species ! integer, dimension(nmeteoLvl) :: varid_meteoLvl ! integer, dimension(nmeteoBnd) :: varid_meteoBnd ! integer :: itau_start ! integer :: itau_end ! integer :: itau_avg ! integer :: itau_ref ! integer :: itau_elapsed ! real*8 :: ddate_avg ! real*8 :: days_elapsed ! integer, dimension(6) :: idate_ref ! integer, dimension(6) :: idate_avg ! character(len=256) :: notes ! character(len=1024) :: disclaimer ! character(len=256) :: email ! character(len=256) :: url ! character(len=256) :: institution ! character(len=256) :: conventions ! character(len=256) :: history ! character(len=256) :: source ! character(len=256) :: sysdate ! character(len=256) :: progstring ! integer, dimension(8) :: isysdate ! integer :: nstationsRegion ! integer, dimension(:), allocatable :: idxRegion ! real, dimension(:,:,:), allocatable :: ncArr ! collect output for ncdf writes ! real, dimension(:,:,:), allocatable :: ncArr_co2_total ! collect total co2 for ncdf writes ! ! --- constants ----------------------- ! character(len=*), parameter :: rname = mname//'/user_output_column_write' ! real*4 :: fillval_r4 = -1e34 ! ! --- begin ------------------------------- ! if(myid /= root) return ! ! subset column list by region ! nstationsRegion = count( stations%region == region ) ! allocate(idxRegion(nstationsRegion)) ! allocate(stationsRegion(nstationsRegion)) ! ii = 0 ! do i = 1, nstations ! there must be a better way to do this... ! if (stations(i)%region == region) then ! ii = ii + 1 ! idxRegion(ii) = i ! end if ! end do ! stationsRegion = stations(idxRegion) ! call tau2date(itaur(region), timestamp(region)%idate_end) ! call date_and_time(values = isysdate) ! write (sysdate, '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2," UTC")') & ! isysdate(1), isysdate(2), isysdate(3), isysdate(5), isysdate(6), isysdate(7) ! ! Standard file information for CarbonTracker output. ! notes = "This file contains CarbonTracker time averaged mole fractions for vertical profiles at select sampling locations." ! disclaimer = "CarbonTracker is an open product of the NOAA Earth System Research \n"//& ! "Laboratory using data from the Global Monitoring Division greenhouse \n"//& ! "gas observational network and collaborating institutions. Model results \n"//& ! "including figures and tabular material found on the CarbonTracker \n"//& ! "website may be used for non-commercial purposes without restriction,\n"//& ! "but we request that the following acknowledgement text be included \n"//& ! "in documents or publications made using CarbonTracker results: \n\n"//& ! " CarbonTracker results provided by NOAA/ESRL,\n"//& ! " Boulder, Colorado, USA, http://carbontracker.noaa.gov\n\n"//& ! "Since we expect to continuously update the CarbonTracker product, it\n"//& ! "is important to identify which version you are using. To provide\n"//& ! "accurate citation, please include the version of the CarbonTracker\n"//& ! "release in any use of these results.\n\n"//& ! "The CarbonTracker team welcomes special requests for data products not\n"//& ! "offered by default on this website, and encourages proposals for\n"//& ! "collaborative activities. Contact us at carbontracker.team@noaa.gov.\n" ! email = "carbontracker.team@noaa.gov" ! url = "http://carbontracker.noaa.gov" ! institution = "NOAA Earth System Research Laboratory" ! conventions = "CF-1.5" ! call getarg (0, progstring, status) ! write(history,'("Created ",a," by ",a,".")') trim(sysdate),trim(progstring) ! source = "CarbonTracker release CT2012" ! ! date/time conversions ! call date2tau(timestamp(region)%idate_start, itau_start) ! call date2tau(timestamp(region)%idate_end, itau_end) ! itau_avg = nint(5.0D-1 * dble(itau_start + itau_end)) ! call tau2date(itau_avg, idate_avg) ! ddate_avg = idate2ddate(idate_avg) ! idate_ref = (/2000, 1, 1, 0, 0, 0/) ! call date2tau(idate_ref, itau_ref) ! itau_elapsed = itau_avg - itau_ref ! days_elapsed = dble(itau_elapsed) / 86400.0D+0 ! write (outfile, '(a,"/column_",a,"_",i4.4,4i2.2,"_",i4.4,4i2.2,".nc")') & ! trim(outdir), trim(region_name(region)), timestamp(region)%idate_start(1:5), timestamp(region)%idate_end(1:5) ! ! create new file ! call MDF_Create( trim(outfile), MDF_NETCDF, MDF_REPLACE, fid, status ) ! IF_ERROR_RETURN(status=1) ! ! global attributes ! call MDF_Put_Att( fid, MDF_GLOBAL, "notes", values=trim(notes), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "disclaimer", values=trim(disclaimer), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "email", values=trim(email), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "url", values=trim(url), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "institution", values=trim(institution), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "conventions", values=trim(conventions), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "history", values=trim(history), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "source", values=trim(source), status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, MDF_GLOBAL, "region_name", values=trim(region_name(region)), status=status ) ! IF_ERROR_RETURN(status=1) ! ! define dimensions ! call MDF_Def_Dim( fid, 'date', MDF_UNLIMITED, dimid_date, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Dim( fid, 'calendar_components', 6, dimid_cal, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Dim( fid, "level", lm(region), dimid_lvl, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Dim( fid, "boundary", lm(region)+1, dimid_bnd, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Dim( fid, "station", nstationsRegion, dimid_stations, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Dim( fid, "character_length", 50, dimid_nchar, status ) ! must be GE character lengths defined in column list derived type ! IF_ERROR_RETURN(status=1) ! ! dimension variables ! call MDF_Def_Var( fid, 'date', MDF_DOUBLE, (/dimid_date/), varid_elapsed, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_elapsed, "units", values="days since 2000-01-01 00:00:00 UTC", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_elapsed, "long_name", values="date", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Var( fid, 'decimal_date', MDF_DOUBLE, (/dimid_date/), varid_dec, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_dec, "units", values="years", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_dec, "_FillValue", values=-1e34, status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Var( fid, 'calendar_components', MDF_INT, (/dimid_cal/), varid_cal_int, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_cal_int, "units", values="none", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_cal_int, "long_name", values="calendar components", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Var( fid, 'date_components', MDF_INT, (/dimid_cal, dimid_date/), varid_date_int, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_date_int, "units", values="none", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_date_int, "long_name", values="integer value calendar components of UTC date", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_date_int, "comment", values="year, month, day, hour, minute, second", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Var( fid, 'station_index', MDF_INT, (/dimid_stations/), varid_column_int, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_column_int, "units", values="none", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_column_int, "long_name", values="station index", status=status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Def_Var( fid, 'lon', MDF_DOUBLE, (/dimid_stations/), varid_lon, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lon, "units", values="degrees east", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lon, "long_name", values="longitude", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lon, "actual_range", values=(/xbeg(region), xend(region)/), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'lat', MDF_DOUBLE, (/dimid_stations/), varid_lat, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lat, "units", values="degrees north", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lat, "long_name", values="latitude", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lat, "actual_range", values=(/ybeg(region), yend(region)/), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'alt', MDF_DOUBLE, (/dimid_stations/), varid_hgt, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_hgt, "units", values="meters", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_hgt, "long_name", values="altitude", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_hgt, "actual_range", values=(/zbeg(region), zend(region)/), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'level', MDF_INT, (/dimid_lvl/), varid_lvl, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lvl, "units", values="none", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lvl, "long_name", values="level", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_lvl, "positive", values="up", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'boundary', MDF_INT, (/dimid_bnd/), varid_bnd, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_bnd, "units", values="none", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_bnd, "long_name", values="boundary", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_bnd, "positive", values="up", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'station_names', MDF_CHAR, (/dimid_nchar, dimid_stations/), varid_desc, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_desc, "long_name", values="station names", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_desc, "standard_name", values="station description", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'station_id', MDF_CHAR, (/dimid_nchar, dimid_stations/), varid_id, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Att( fid, varid_id, "long_name", values="station id", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_id, "standard_name", values="station code", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'oro', MDF_FLOAT, (/dimid_stations/), varid_oro, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_oro, "units", values="m^s/s^2", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_oro, "long_name", values="orography", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_oro, "standard_name", values="surface geopotential", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Def_Var( fid, 'co2', MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_total, status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_total, "units", values="micromol mol-1", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_total, "long_name", values="mole fraction of carbon dioxide in air", status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_total, "comment", values="total carbon dioxide estimate", status=status ) ! IF_NOTOK_RETURN(status=1) ! do i = 1, nspecies ! call MDF_Def_Var( fid, trim(names(species(i))), MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_species(i, writeMethod), status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_species(i, writeMethod), "units", values='micromol mol-1', status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_species(i, writeMethod), "_FillValue", values=fillval_r4, status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_species(i, writeMethod), "long_name", values='mole fraction of carbon dioxide in air', status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_species(i, writeMethod), "comment", values=comment_species(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! enddo ! do i = 1, nmeteoLvl ! if (i .eq. 4) then ! call MDF_Def_Var( fid, name_meteoLvl(i), MDF_FLOAT, (/dimid_stations, dimid_date/), varid_meteoLvl(i), status ) ! else ! call MDF_Def_Var( fid, name_meteoLvl(i), MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_meteoLvl(i), status ) ! endif ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoLvl(i), "units", values=units_meteoLvl(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoLvl(i), "long_name", values=long_name_meteoLvl(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoLvl(i), "standard_name", values=standard_name_meteoLvl(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! enddo ! do i = 1, nmeteoBnd ! call MDF_Def_Var( fid, name_meteoBnd(i), MDF_FLOAT, (/dimid_stations, dimid_bnd, dimid_date/), varid_meteoBnd(i), status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoBnd(i), "units", values=units_meteoBnd(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoBnd(i), "long_name", values=long_name_meteoBnd(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! call MDF_Put_Att( fid, varid_meteoBnd(i), "standard_name", values=standard_name_meteoBnd(i), status=status ) ! IF_NOTOK_RETURN(status=1) ! enddo ! ! finished definition ! call MDF_EndDef( fid, status ) ! IF_ERROR_RETURN(status=1) ! ! write variables ! call MDF_Put_Var( fid, varid_elapsed, (/(days_elapsed)/), status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_dec, (/(ddate_avg)/), status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_cal_int, (/(i,i=1,size(idate_avg))/), status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_date_int, idate_avg, status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_desc, (/(trim(stationsRegion(i)%desc),i=1,nstationsRegion)/), status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_id, (/(trim(stationsRegion(i)%id),i=1,nstationsRegion)/), status ) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_column_int, (/(i,i=1,size(stationsRegion))/), status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_lat, stationsRegion%lat, status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_lon, stationsRegion%lon, status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_lvl, (/(i,i=1,lm(region))/), status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_bnd, (/(i,i=1,lm(region)+1)/), status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_hgt, stationsRegion%height, status) ! IF_ERROR_RETURN(status=1) ! call MDF_Put_Var( fid, varid_oro, stationsRegion%oro, status) ! IF_ERROR_RETURN(status=1) ! allocate(ncArr(nstationsRegion, lm(region), 1)) ! allocate(ncArr_co2_total(nstationsRegion, lm(region), 1)) ! ncArr_co2_total = 0.0 ! do i = 1, nspecies ! do q = 1, nstationsRegion ! ncArr(q,:,1) = stationsRegion(q)%c_avg(i,writeMethod,:) ! ! subtract the background tracer from the others ! if (i .gt. 1) then ! ncArr(q,:,1) = ncArr(q,:,1) - stationsRegion(q)%c_avg(1,writeMethod,:) ! endif ! ncArr_co2_total(q,:,1) = ncArr_co2_total(q,:,1) + ncArr(q,:,1) ! enddo ! call MDF_Put_Var( fid, varid_species(i,writeMethod), ncArr, status ) ! IF_ERROR_RETURN(status=1) ! enddo ! call MDF_Put_Var( fid, varid_total, ncArr_co2_total, status ) ! IF_ERROR_RETURN(status=1) ! deallocate(ncArr_co2_total) ! do i = 1, nmeteoLvl ! do q = 1, nstationsRegion ! ncArr(q,:,1) = stationsRegion(q)%meteoLvl(i,:) ! enddo ! if (i .eq. 4) then ! call MDF_Put_Var( fid, varid_meteoLvl(i), ncArr(:,1,:), status ) ! else ! call MDF_Put_Var( fid, varid_meteoLvl(i), ncArr, status ) ! endif ! IF_ERROR_RETURN(status=1) ! enddo ! deallocate(ncArr) ! allocate(ncArr(nstationsRegion, lm(region)+1, 1)) ! do i = 1, nmeteoBnd ! do q = 1, nstationsRegion ! ncArr(q,:,1) = stationsRegion(q)%meteoBnd(i,:) ! enddo ! call MDF_Put_Var( fid, varid_meteoBnd(i), ncArr, status ) ! IF_ERROR_RETURN(status=1) ! enddo ! deallocate(ncArr) ! ! close file ! call MDF_Close( fid, status ) ! IF_ERROR_RETURN(status=1) ! write (gol, '("'//trim(rname)//' Done writing arrays and closing output file.")'); call goPr ! timestamp(region)%idate_start = timestamp(region)%idate_end ! end of interval becomes start of next interval ! deallocate(idxRegion) ! deallocate(stationsRegion) ! ok status = 0 end subroutine user_output_column_write subroutine user_output_column_reset( region ) ! use ParTools, only: root_t, myid ! implicit none integer, intent(in) :: region ! integer :: q ! if (myid /= root_t) return ! do q=1, nstations ! stations(q)%c_avg(:,:,:) = 0.0 ! stations(q)%c_std(:,:,:) = 0.0 ! stations(q)%meteoLvl(:,:) = 0.0 ! stations(q)%meteoBnd(:,:) = 0.0 ! stations(q)%ndata = 0 ! stations(q)%accum_weight = 0.0 ! end do ! evaluated = .FALSE. end subroutine user_output_column_reset subroutine user_output_column_done( status ) ! use Partools, only : myid, root ! use dims, only : nregions ! implicit none ! --- in/out ------------------------ integer, intent(out) :: status ! --- variables --------------------------- ! integer :: region ! integer :: q ! ! --- constants --------------------------- ! character(len=*), parameter :: rname = mname//'/user_output_column_done' ! ! --- begin ------------------------------- ! if (myid /= root) return ! do region = 1, nregions ! call user_output_column_evaluate( status ) ! IF_NOTOK_RETURN(status=1) ! call user_output_column_write( region, status ) ! IF_NOTOK_RETURN(status=1) ! end do ! do q = 1, nstations ! deallocate(stations(q)%c_avg) ! deallocate(stations(q)%c_std) ! enddo ! deallocate(stations) ! ok status = 0 end subroutine user_output_column_done end module user_output_column