#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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 #define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if ! #include "tm5.inc" !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: USER_OUTPUT_PDUMP ! ! !DESCRIPTION: ! ! Module to deal with time-series output. Output are in NetCDF-4 and use CF ! conventions. The following output are available: ! ! - one file with grid definition ! - one file with time series of some met fields (pressure, temperature, winds, ...) ! - one or more files with time series of some tracers ! - one or two files with Local Time output for some tracers ! - one file with time series of wet and dry depositions ! - one file with time series of deposition velocity ! ! If the macro (cpp) "tropomi" is used, then the temperature and extra attributes added to the vmr (tracers) datasets. ! ! Activation, tracers to account for, time step of the series, are set in the ! rcfile, following this template : ! ! ! SAMPLE RCFILE ! ! output.pdump : T ! output.pdump.dataset.author : John Doe ! output.pdump.dataset.institution : MyFirm, Anytown, USA ! output.pdump.dataset.version : GEMS GRG; era2003 simulation ! tropomi only: ! output.pdump.tropomi.tm5version : v4 ! output.pdump.tropomi.institution : KNMI ! output.pdump.tropomi.tm5reference : Huijnen et al., ACP ! output.pdump.tropomi.authoremail : Doe@john.com ! output.pdump.tropomi.datasetname : "S5P_AUX_CTMFCT" or "S5P_AUX_CTMANA" ! ! output.pdump.fname.model : TM5 ! output.pdump.fname.expid : V2 ! output.pdump.fname.grid.300x200 : 3x2 ! short name, required if there is zoom regions ! output.pdump.fname.grid.100x100 : 1x1 ! ! output.pdump.griddef.apply : T ! ! output.pdump.tp.apply : T ! output.pdump.tp.dhour : 1 ! ! output.pdump.vmr.n : 3 ! ! output.pdump.vmr.001.apply : T ! output.pdump.vmr.001.fname : vmr3 ! output.pdump.vmr.001.dhour : 3 ! output.pdump.vmr.001.tracers : SO2 NOy CH4 OH HNO3 PAN H2O2 Radon Lead ! ! output.pdump.vmr.002.apply : T ! output.pdump.vmr.002.fname : vmr1 ! output.pdump.vmr.002.dhour : 1 ! output.pdump.vmr.002.tracers : O3 O3s CO NO2 NO CH2O ! ! output.pdump.vmr.003.apply : F ! output.pdump.vmr.003.fname : vmra ! output.pdump.vmr.003.dhour : 3 ! output.pdump.vmr.003.tracers : SO4 NO3_A BC BCS POM SS1_N SS1_M SS2_N SS2_M SS3_N SS3_M DUST2_N DUST2_M DUST3_N DUST3_M ! ! output.pdump.lt.apply : T ! output.pdump.lt.tracers : O3 ! output.pdump.lt.localtime : 2 ! ! output.pdump.lt2.apply : F ! output.pdump.lt2.tracers : ! output.pdump.lt2.localtime : ! ! output.pdump.depositions.apply : F ! output.pdump.depositions.dhour : 3 ! output.pdump.depositions.tracers : O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2 NOy ! ! output.pdump.depvels.apply : F ! output.pdump.depvels.dhour : 3 ! output.pdump.depvels.tracers : O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2 ! !\\ !\\ ! !INTERFACE: ! MODULE USER_OUTPUT_PDUMP ! ! !USES: ! use partools, only : isRoot use GO, only : gol, goPr, goErr, goLabel use GO, only : TDate, IncrDate, NewDate use GO, only : operator(+), SystemDate, Get use dims, only : nregions, idatee, idatei, okdebug, nread use chem_param, only : ntrace use chem_param, only : iNOx, iHNO3, iPAN, iOrgNtr #ifdef with_m7 use chem_param, only : iNO3_a use chem_param, only : iSO4nus, iSO4ais, iSO4acs, iSO4cos use chem_param, only : iBCais, iBCacs, iBCcos, iBCaii use chem_param, only : iPOMais, iPOMacs, iPOMcos, iPOMaii use chem_param, only : iDUacs, iDUcos, iDUaci, iDUcoi use chem_param, only : iSScos, iSSacs #endif USE MDF USE TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Output_PDUMP_Init public :: Output_PDUMP_Step public :: Output_PDUMP_Done ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'user_output_pdump' character(len=*), parameter :: outfileversnr = '0.1' integer, parameter :: time_reftime6(6) = (/1950,01,01,00,00,00/) ! reference time character(len=*), parameter :: time_units = 'days since 1950-01-01 00:00:00' ! ! NOy is not a standard tracer field, but sum of some transported tracers: ! NOx HNO3 PAN orgntr NO3_a ! where NOx is the sum of short lived tracers: ! NOx = NO + NO2 + NO3 + HNO4 + 2*N2O5 ! #ifdef with_m7 integer, parameter :: iNOy = ntrace + 1 integer, parameter :: nNOyt = 5 integer, parameter :: iNOyt(nNOyt) = (/ iNOx, iHNO3, iNO3_a, iPAN, iOrgNtr /) integer, parameter :: iSO4 = ntrace + 2 integer, parameter :: nSO4t = 4 integer, parameter :: iSO4t(nSO4t) = (/ iSO4nus, iSO4ais, iSO4acs, iSO4cos /) integer, parameter :: iBC = ntrace + 3 integer, parameter :: nBCt = 4 integer, parameter :: iBCt(nBCt) = (/ iBCais, iBCacs, iBCcos, iBCaii /) integer, parameter :: iPOM = ntrace + 4 integer, parameter :: nPOMt = 4 integer, parameter :: iPOMt(nPOMt) = (/ iPOMais, iPOMacs, iPOMcos, iPOMaii /) integer, parameter :: iSS = ntrace + 5 integer, parameter :: nSSt = 2 integer, parameter :: iSSt(nSSt) = (/ iSSacs, iSScos /) integer, parameter :: iDU = ntrace + 6 integer, parameter :: nDUt = 4 integer, parameter :: iDUt(nDUt) = (/ iDUacs, iDUcos, iDUaci, iDUcoi /) #else integer, parameter :: iNOy = ntrace + 1 integer, parameter :: nNOyt = 4 integer, parameter :: iNOyt(nNOyt) = (/ iNOx, iHNO3, iPAN, iOrgNtr /) #endif ! ! !PRIVATE TYPES: ! type TPdumpFile_GridDef integer :: trec integer :: ncid integer :: dimid_scalar, dimid_lon, dimid_lat, dimid_lev, dimid_levi integer :: varid_lon, varid_lat, varid_time, varid_date integer :: varid_gridbox_area integer :: varid_a, varid_b integer :: varid_a_bnds, varid_b_bnds integer :: varid_p0 !integer :: varid_ps !integer :: varid_geo_height end type TPdumpFile_GridDef type TPdumpFile_TP integer :: trec integer :: dhour integer :: ncid integer :: dimid_lon, dimid_lat, dimid_lev, dimid_time, dimid_datelen integer :: varid_lon, varid_lat, varid_lev, varid_time, varid_date integer :: varid_ps integer :: varid_surface_temp integer :: varid_orog integer :: varid_geop integer :: varid_pressure integer :: varid_temp integer :: varid_humid integer :: varid_u, varid_v, varid_w real, allocatable :: data3d(:,:,:,:,:) real, allocatable :: data2d(:,:,:,:) real, allocatable :: time(:) real, allocatable :: date(:,:) end type TPdumpFile_TP type TPdumpFile_VMR integer :: trec, n_rec logical :: apply real :: dhour integer :: dsec character(len=256) :: tracer_names integer :: ncid integer :: dimid_lon, dimid_lat, dimid_lev, dimid_levi, dimid_time, dimid_datelen integer :: varid_lon, varid_lat, varid_lev, varid_time, varid_date integer :: varid_ps integer :: varid_a_bnds, varid_b_bnds integer :: ntr integer :: itr(ntrace) character(len=8) :: name_tr(ntrace) #ifdef with_m7 logical :: lpmx(ntrace) real :: sizepmx(ntrace) #endif integer :: varid_tr(ntrace) character(len=4) :: varid_type(ntrace) real, allocatable :: data3d(:,:,:,:,:) real, allocatable :: sp(:,:,:) real, allocatable :: time(:) real, allocatable :: date(:,:) real, allocatable :: data3d_t(:,:,:,:) integer :: varid_temp #ifdef tropomi integer :: varid_hyai, varid_hybi, varid_hyam, varid_hybm integer :: varid_hgt integer :: varid_ltropo real, allocatable :: data2d_hgt(:,:) integer, allocatable:: data2d_ltropo(:,:,:) #endif end type TPdumpFile_VMR type TPdumpFile_LT integer :: trec character(len=256) :: tracer_names integer :: ncid integer :: local_time integer :: dimid_lon, dimid_lat, dimid_lev, dimid_time, dimid_datelen integer :: varid_lon, varid_lat, varid_lev, varid_time, varid_date integer :: varid_ps integer :: ntr integer :: itr(ntrace) character(len=8) :: name_tr(ntrace) integer :: varid_tr(ntrace) real,allocatable :: accu(:,:,:,:) real,allocatable :: naccu(:,:) real,allocatable :: p_accu(:,:) real,allocatable :: np_accu(:) #ifdef with_m7 logical :: laod(ntrace) real :: wavel(ntrace) #endif end type TPdumpFile_LT type TPdumpFile_DEPS integer :: trec integer :: dhour character(len=256) :: tracer_names integer :: ncid integer :: dimid_lon, dimid_lat, dimid_time, dimid_datelen integer :: varid_lon, varid_lat, varid_time, varid_date, varid_accum integer :: ntr integer :: itr(ntrace) character(len=8) :: name_tr(ntrace) integer :: varid_ddep(ntrace) real, pointer :: ddep_budget(:,:,:) logical :: with_wdep(ntrace) integer :: varid_wdep(ntrace) real, pointer :: wdep_budget(:,:,:) type(TDate) :: t0_budget real, allocatable :: data2d_dry(:,:,:,:) real, allocatable :: data2d_wet(:,:,:,:) real, allocatable :: time(:), dt(:) real, allocatable :: date(:,:) end type TPdumpFile_DEPS type TPdumpFile_DEPV integer :: trec integer :: dhour character(len=256) :: tracer_names integer :: ncid integer :: dimid_lon, dimid_lat, dimid_time, dimid_datelen integer :: varid_lon, varid_lat, varid_time, varid_date integer :: ntr integer :: itr(ntrace) character(len=8) :: name_tr(ntrace) integer :: varid_tr(ntrace) real, allocatable :: data2d(:,:,:,:) real, allocatable :: time(:) real, allocatable :: date(:,:) end type TPdumpFile_DEPV ! --- var ----------------------------- integer :: fid ! file id for IF_NOTOK_MDF macro integer :: access_mode ! netcdf-4 access mode integer :: curr_day(nregions,3) logical :: firstday logical :: lastday ! it is last day and not a full day (ie day does not end at 00 of next day) character(len=32) :: fname_model character(len=8) :: fname_expid, meteo_class character(len=32) :: fname_grid(nregions) character(len=256) :: dataset_author, institution, dataset_version #ifdef tropomi character(len=256) :: tropomi_authoremail, tropomi_tm5_reference, tropomi_institution character(len=256) :: tropomi_tm5_version, tropomi_dataset_name character(len=15) :: tropomi_date_start, tropomi_date_stop, tropomi_date_create #endif logical, save :: griddef_apply type(TPdumpFile_GridDef), save :: RF_GridDef(nregions) logical, save :: tp_apply integer :: tp_dhour, n_tp_rec type(TPdumpFile_TP), save :: RF_TP(nregions) integer, save :: nvmr logical, allocatable :: vmr_apply(:) real, allocatable :: vmr_sregbord(:,:) character(len=16), allocatable :: vmr_fname(:) real, allocatable :: vmr_dhour(:) character(len=256), allocatable :: vmr_tracer_names(:) type(TPdumpFile_VMR), allocatable, save :: RF_VMR(:,:) logical, save :: lt_apply character(len=16) :: lt_fname character(len=256) :: lt_tracer_names integer :: lt_localtime type(TPdumpFile_LT), save :: RF_LT(nregions) logical, save :: lt2_apply character(len=16) :: lt2_fname character(len=256) :: lt2_tracer_names integer :: lt2_localtime type(TPdumpFile_LT), save :: RF_LT2(nregions) logical, save :: deps_apply character(len=16) :: deps_fname integer :: deps_dhour, n_deps_rec character(len=256) :: deps_tracer_names type(TPdumpFile_DEPS), save :: RF_DEPS(nregions) logical, save :: depv_apply character(len=16) :: depv_fname integer :: depv_dhour, n_depv_rec character(len=256) :: depv_tracer_names type(TPdumpFile_DEPV), save :: RF_DEPV(nregions) ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - revised older RETRO ouptut : ! add 2nd local time, regional output, ! handle aerosol tracers and M7 ! 10 Jul 2012 - Ph. Le Sager - switch from pnetcdf to netcdf4_par (through ! MDF); get rid of the with_tendencies code. ! 12 Nov 2012 - Ph. Le Sager - adapted for lon-lat MPI decomposition. ! - get rid of unlimited dimensions so we can ! write in collective mode. ! - store series to write them only at end-of-day ! to speed-up code ! 10 Oct 2013 - Ph. Le Sager - fixed GET_N_TIME_RECORDS and several 'init' ! and write' routines. ! 14 Apr 2014 - Ph. Le Sager + JEW - tropomi add-ons in VMR: Temperature, ! As, Bs, better CF ! 8 October 2014 - H. Eskes - changes in tropomi output (based on the "tropomi" macro) ! ! !REMARKS: ! ! (1) Initially called RETRO output for GEMS GRG, the module has been adapted ! for CLIMAQS project and renamed PDUMP. ! (2) Previous remarks "as is": ! - longitudes from [0,360] ? ! this is imposible for zoom area's such as for the heatwave ! - levels from surface to top ! - time from 1950-01-01 00:00 ! (3) This is supposed to work with netcdf4_parallel. You cannot use ! MPI with a non-parallel version of netcdf4 here. ! (4) The parallel writing is done in COLLECTIVE mode, but remain ! expensive on some system. Possible optimization : output one file ! per month (chunk/leg), and/or per field, and/or per processor. ! (5) Switch in nstep for DEPS data should work for full days. Not tested ! for partial days. ! ! !TODO: ! - test with M7 tracers. Which ones? ! - in LT_WRITE : AOD if m7 needs to be coded ! - in RF_VMR_INIT : match tracer with CF standard names for some aerosols ! (dust,...) ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: GET_N_TIME_RECORDS ! ! !DESCRIPTION: return number of time steps for a daily timeseries file !\\ !\\ ! !INTERFACE: ! FUNCTION GET_N_TIME_RECORDS( date, dsec, isDEPS, mess ) ! ! !USES: ! USE GO , only : TDate, NewDate, rTotal, operator(-) ! ! !RETURN VALUE: ! integer :: get_n_time_records ! ! !INPUT PARAMETERS: ! integer, intent(in) :: date(6) ! 1st time step of the day (timestart basically) integer, intent(in) :: dsec ! time step for timeseries in sec (should divide 24*3600, be divided by ndyn/2) logical, optional, intent(in) :: isDEPS ! to differentiate b/w DEPS and others character(len=*), optional, intent(in) :: mess ! message (if okdebug) ! ! !REVISION HISTORY: ! 9 Nov 2012 - Ph. Le Sager - v0 ! 9 Oct 2013 - Ph. Le Sager - fix to work with default "output.after.step: v" ! 15 Jul 2014 - Ph. Le Sager - works with seconds instead of hours ! ! !REMARKS: ! - dynamic timestep cannot be LARGER than timestep of timeseries, with notable exception ! of dynamic timestep = 2*timeseries_timestep. ! ! !TODO: ! - check if anything changes with other possible values of "output.after.step" ! !EOP !------------------------------------------------------------------------ !BOC integer :: is, ie, delta, dynstep logical :: deps type(TDate) :: t, t0 real :: time ! Type of record (standard=vmr, tp, depv) or special (deps) deps=.false. if (present(isDEPS)) deps=isDEPS ! Start index delta=date(4)*3600+date(5)*60+date(6) ! 0, unless start of the run is not at 00:00:00 if (deps) delta=delta + nread ! one DYNAMIC time step done to output something if (modulo(delta,dsec)==0) then is=delta/dsec else is=(delta+dsec)/dsec end if ! End index for daily file (nread=dynamic time step read from rc) ie = (24*3600 - nread/2) / dsec if (deps) then ! there will be an extra step if run goes further than midnight t0 = NewDate( time6=date ) t = NewDate( time6=idatee ) time = rTotal( t - t0, 'day' ) if (time > 1) ie=24*3600/dsec end if ! Case of "last day stopping before midnite". (Need testing for DEPS) if (lastday) ie=(idatee(4)*3600+idatee(5)*60+idatee(6)-nread/2)/dsec ! length get_n_time_records = ie-is+1 if(okdebug)then if (present(mess))then write(gol,*) 'GET_N_TIME_RECORDS -'//trim(mess); call goPr end if write(gol,*) "GET_N_TIME_RECORDS - is, ie, deps, firstday, lastday, get_n_time_records:" ; call goPr write(gol,*) "GET_N_TIME_RECORDS - ", is, ie, deps, firstday, lastday, get_n_time_records ; call goPr write(gol,*) "GET_N_TIME_RECORDS - date, dsec, nread ", date, dsec, nread ; call goPr write(gol,*) "GET_N_TIME_RECORDS - idateE ", idatee ; call goPr end if return END FUNCTION GET_N_TIME_RECORDS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_PDUMP_INIT ! ! !DESCRIPTION: reads rc file keys relevant for pdump !\\ !\\ ! !INTERFACE: ! SUBROUTINE OUTPUT_PDUMP_INIT( rcF, dsec_min, status ) ! ! !USES: ! use GO, only : TrcFile, ReadRc use MeteoData, only : lli, set use MeteoData, only : sp_dat, oro_dat, temper_dat, humid_dat, pu_dat, pv_dat use MeteoData, only : mfw_dat, gph_dat, t2m_dat ! ! !INPUT/OUTPUT PARAMETERS: ! type(TrcFile), intent(inout) :: rcF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: dsec_min ! smallest timeseries period in sec integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - upgrade from RETRO to PDUMP ! 8 Nov 2012 - Ph. Le Sager - added access mode switch ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Output_PDUMP_Init' ! --- local ------------------------------ integer :: region character(len=64) :: key character(len=3) :: nr integer :: ivmr ! --- begin ------------------------------- call goLabel(rname) #ifdef MPI #ifdef with_netcdf4_par access_mode = MDF_COLLECTIVE #else write(gol,'("Time Series output (PDUMP) requires netcdf4 with parallel access enabled")') ; call goErr TRACEBACK status=1; return #endif #else access_mode = MDF_INDEPENDENT #endif ! which day firstday = .true. lastday = .true. ! lowest time frequency in sec dsec_min = 999999 if (any(idatei(1:3)/=idatee(1:3))) lastday=.false. ! i.e. at least one full day ! dataset keys: call ReadRc( rcF, 'output.pdump.dataset.author' , dataset_author , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.dataset.institution', institution , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.dataset.version' , dataset_version , status ) IF_NOTOK_RETURN(status=1) #ifdef tropomi call ReadRc( rcF, 'output.pdump.tropomi.tm5version', tropomi_tm5_version , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.tropomi.institution', tropomi_institution , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.tropomi.tm5reference', tropomi_tm5_reference , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.tropomi.authoremail', tropomi_authoremail , status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.tropomi.datasetname', tropomi_dataset_name , status ) IF_NOTOK_RETURN(status=1) #endif ! filename keys: call ReadRc( rcF, 'output.pdump.fname.model', fname_model, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.fname.expid', fname_expid, status ) IF_NOTOK_RETURN(status=1) ! prefix grid name in case of zooming regions: if ( nregions > 1 ) then ! loop over regions: do region = 1, nregions ! short grid name from rcfile: call ReadRc( rcF, 'output.pdump.fname.grid.'//trim(lli(region)%name), key, status ) IF_NOTOK_RETURN(status=1) ! fill grid extenstion to file names: fname_grid(region) = '-'//trim(key) end do else ! empty fname_grid = '' end if ! griddef file ? call ReadRc( rcF, 'output.pdump.griddef.apply', griddef_apply, status ) IF_NOTOK_RETURN(status=1) ! temperature, pressure, etc file ? call ReadRc( rcF, 'output.pdump.tp.apply', tp_apply, status ) IF_NOTOK_RETURN(status=1) if (tp_apply) then ! ensure that required meteo is loaded do region=1,nregions call Set( sp_dat(region), status, used=.true. ) call Set( oro_dat(region), status, used=.true. ) call Set( temper_dat(region), status, used=.true. ) call Set( t2m_dat(region), status, used=.true. ) call Set( humid_dat(region), status, used=.true. ) call Set( pu_dat(region), status, used=.true. ) call Set( pv_dat(region), status, used=.true. ) call Set( mfw_dat(region), status, used=.true. ) call Set( gph_dat(region), status, used=.true. ) ! used to compute vertical wind end do ! time resolution (1 hour by default) call ReadRc( rcF, 'output.pdump.tp.dhour', tp_dhour, status, default=1 ) IF_ERROR_RETURN(status=1) dsec_min = tp_dhour*3600 end if ! VMR files call ReadRc( rcF, 'output.pdump.vmr.n', nvmr, status ) ! number of vmr files to be written IF_NOTOK_RETURN(status=1) if ( nvmr < 0 ) then write (gol,'("strange specification of number of vmr files : ",i6)') nvmr; call goErr TRACEBACK; status=1; return end if ! meteo call ReadRc( rcF, 'my.meteo.class', meteo_class, status, default='unknown' ) IF_ERROR_RETURN(status=1) ! write any vmr files ? if ( nvmr > 0 ) then ! storage: allocate( vmr_apply(nvmr) ) ; vmr_apply = .false. allocate( vmr_fname(nvmr) ) ; vmr_fname = '' allocate( vmr_dhour(nvmr) ) ; vmr_dhour = -1. allocate( vmr_tracer_names(nvmr) ) ; vmr_tracer_names = '' allocate( vmr_sregbord(nvmr,4) ) ; vmr_sregbord = -999.9 allocate( RF_VMR(nregions,nvmr) ) #ifdef tropomi do region=1,nregions call Set( temper_dat(region), status, used=.true. ) call Set( gph_dat(region), status, used=.true. ) ! used to compute surface altitude end do #endif ! loop over vmr files: do ivmr = 1, nvmr ! number in rc keys: write (nr,'(i3.3)') ivmr ! apply ? call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.apply', vmr_apply(ivmr), status ) IF_NOTOK_RETURN(status=1) rf_vmr(:,ivmr)%apply = vmr_apply(ivmr) ! proceed ? if ( vmr_apply(ivmr) ) then ! first part of filename: call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.fname', vmr_fname(ivmr), status ) IF_NOTOK_RETURN(status=1) ! time resolution: call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.dhour', vmr_dhour(ivmr), status ) IF_NOTOK_RETURN(status=1) ! here is the catch: fractional hour for step should divide 3600 dsec_min = min( dsec_min, int(vmr_dhour(ivmr)*3600) ) ! tracers to be written: call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.tracers', vmr_tracer_names(ivmr), status ) IF_NOTOK_RETURN(status=1) end if ! apply ? end do ! vmr numbers ! required meteo if (any(vmr_apply)) then do region=1,nregions call Set( sp_dat(region), status, used=.true. ) end do end if end if ! nvmr > 0 ! --------------------- ! local time: ! --------------------- ! file 1 lt_fname = 'lt' call ReadRc( rcF, 'output.pdump.lt.apply', lt_apply, status ) IF_NOTOK_RETURN(status=1) if ( lt_apply ) then call ReadRc( rcF, 'output.pdump.lt.tracers', lt_tracer_names, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.lt.localtime', lt_localtime, status ) IF_NOTOK_RETURN(status=1) end if ! file 2 lt2_fname = 'lt2' call ReadRc( rcF, 'output.pdump.lt2.apply', lt2_apply, status ) IF_NOTOK_RETURN(status=1) if ( lt2_apply ) then call ReadRc( rcF, 'output.pdump.lt2.tracers', lt2_tracer_names, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.lt2.localtime', lt2_localtime, status ) IF_NOTOK_RETURN(status=1) end if if (lt_apply .or. lt2_apply) then do region=1,nregions call Set( sp_dat(region), status, used=.true. ) end do end if ! --------------------- ! deposition fluxes: ! --------------------- deps_fname = 'depositions' call ReadRc( rcF, 'output.pdump.depositions.apply', deps_apply, status ) IF_NOTOK_RETURN(status=1) if ( deps_apply ) then #ifdef with_budgets call ReadRc( rcF, 'output.pdump.depositions.dhour', deps_dhour, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.depositions.tracers', deps_tracer_names, status ) IF_NOTOK_RETURN(status=1) dsec_min = min( dsec_min, deps_dhour*3600) #else write(gol,*) "timeseries of deposition fluxes requires using 'with_budget' macro" ; call goErr status=1 ; TRACEBACK ; return #endif end if ! --------------------- ! deposition velocities: ! --------------------- depv_fname = 'depvels' call ReadRc( rcF, 'output.pdump.depvels.apply', depv_apply, status ) IF_NOTOK_RETURN(status=1) if ( depv_apply ) then #ifdef with_budgets call ReadRc( rcF, 'output.pdump.depvels.dhour', depv_dhour, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'output.pdump.depvels.tracers', depv_tracer_names, status ) IF_NOTOK_RETURN(status=1) dsec_min = min( dsec_min, depv_dhour*3600) #else write(gol,*) "timeseries of deposition velocities requires using 'with_budget' macro" ; call goErr status=1 ; TRACEBACK ; return #endif end if ! no files open yet curr_day = -1 call goLabel() ! ok status = 0 END SUBROUTINE OUTPUT_PDUMP_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_PDUMP_STEP ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE OUTPUT_PDUMP_STEP( region, idate_f, status ) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! ! !REMARKS: ! (1) called every hour. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Output_PDUMP_Step' ! --- begin ------------------------------- call goLabel(rname) !---------------------- ! close if necessary !---------------------- ! if a file is open, and it is a new day if ( all(curr_day(region,:) > 0) .and. any(idate_f(1:3) /= curr_day(region,:)) ) then ! write in previous day file end-of-interval data call PDUMP_Files_Write2( region, idate_f, status ) IF_NOTOK_RETURN(status=1) ! close all call PDUMP_Files_Close( region, status ) IF_NOTOK_RETURN(status=1) ! no files open ... curr_day(region,:) = -1 firstday = .false. end if !---------------------- ! open if necessary !---------------------- if ( any(curr_day(region,:) < 0) ) then if (all(idate_f(1:3)==idatee(1:3))) lastday=.true. ! means last day is not a full day write(gol,*) "U_O_Pdump open [idate_f, last day] = ", idate_f, lastday; call goPr call PDUMP_Files_Open( region, idate_f, status ) IF_NOTOK_RETURN(status=1) ! store date of current day curr_day(region,:) = idate_f(1:3) end if !---------------------- ! write !---------------------- call PDUMP_Files_Write( region, idate_f, status ) IF_NOTOK_RETURN(status=1) ! if not midnight, write end-of-interval data if ( any(idate_f(4:6) > 0) ) then call PDUMP_Files_Write2( region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if !---------------------- ! done !---------------------- call goLabel() status = 0 END SUBROUTINE OUTPUT_PDUMP_STEP !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_PDUMP_DONE ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE OUTPUT_PDUMP_DONE( status ) ! ! !USES: ! use dims, only : itaur use datetime, only : tau2date ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 31 Aug 2012 - P. Le Sager - reverse order in which regions are dealt with (MDF issue) ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Output_PDUMP_Done' integer :: region integer,dimension(6) :: idate_f ! --- begin ------------------------------- ! close files: do region = nregions, 1, -1 ! write end of interval DEPS data (requires that DEPS nstep is calculated with .false. -see RF_DEPS_Init-) call tau2date(itaur(region),idate_f) call PDUMP_Files_Write2( region, idate_f, status ) IF_NOTOK_RETURN(status=1) call PDUMP_Files_Close( region, status ) IF_NOTOK_RETURN(status=1) end do ! clear: if ( nvmr > 0 ) then deallocate( vmr_apply ) deallocate( vmr_fname ) deallocate( vmr_dhour ) deallocate( vmr_tracer_names ) deallocate( vmr_sregbord ) deallocate( RF_VMR ) end if ! ok status = 0 END SUBROUTINE OUTPUT_PDUMP_DONE !EOC ! ******************************************************************** ! *** ! *** open/write/close pdump files ! *** ! ******************************************************************** !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PDUMP_FILES_OPEN ! ! !DESCRIPTION: call init method of each output file. !\\ !\\ ! !INTERFACE: ! subroutine PDUMP_Files_Open( region, idate_f, status ) ! ! !USES: ! use global_data, only : outdir ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/PDUMP_Files_Open' ! --- local ------------------------------- integer :: ivmr ! --- begin ------------------------------- ! grid definition: if ( griddef_apply ) then call RF_GridDef_Init( RF_GridDef(region), outdir, fname_model, fname_expid, region, status ) IF_NOTOK_RETURN(status=1) end if ! dynamics: if ( tp_apply ) then call RF_TP_Init ( RF_TP(region) , outdir, fname_model, fname_expid, & region, idate_f, tp_dhour, status ) IF_NOTOK_RETURN(status=1) end if ! tracer concentrations: do ivmr = 1, nvmr if ( .not. vmr_apply(ivmr) ) cycle call RF_VMR_Init( RF_VMR(region,ivmr), outdir, fname_model, fname_expid, & vmr_fname(ivmr), region, idate_f, & vmr_dhour(ivmr), vmr_tracer_names(ivmr), status ) IF_NOTOK_RETURN(status=1) vmr_apply(ivmr) = rf_vmr(region,ivmr)%apply end do ! lt output: if ( lt_apply ) then call RF_LT_Init( RF_LT(region), outdir, fname_model, fname_expid, & lt_fname, region, idate_f, & lt_localtime, lt_tracer_names, status ) IF_NOTOK_RETURN(status=1) end if if ( lt2_apply ) then call RF_LT_Init( RF_LT2(region), outdir, fname_model, fname_expid, & lt2_fname, region, idate_f, & lt2_localtime, lt2_tracer_names, status ) IF_NOTOK_RETURN(status=1) end if #ifdef with_budgets ! deposition fluxes: if ( deps_apply ) then call RF_DEPS_Init( RF_DEPS(region), outdir, fname_model, fname_expid, & deps_fname, region, idate_f, & deps_dhour, deps_tracer_names, status ) IF_NOTOK_RETURN(status=1) end if ! deposition velocities: if ( depv_apply ) then call RF_DEPV_Init( RF_DEPV(region), outdir, fname_model, fname_expid, & depv_fname, region, idate_f, & depv_dhour, depv_tracer_names, status ) IF_NOTOK_RETURN(status=1) end if #endif ! ok status = 0 END SUBROUTINE PDUMP_FILES_OPEN !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PDUMP_FILES_WRITE ! ! !DESCRIPTION: call write method for each output file. !\\ !\\ ! !INTERFACE: ! SUBROUTINE PDUMP_FILES_WRITE( region, idate_f, status ) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/PDUMP_Files_Write' integer :: ivmr ! --- begin ------------------------------- ! grid definition: if ( griddef_apply ) then call RF_GridDef_Write( RF_GridDef(region), region, status ) IF_NOTOK_RETURN(status=1) end if ! dynamics: if ( tp_apply ) then call RF_TP_Write( RF_TP(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if ! tracer fields: do ivmr = 1, nvmr if ( .not. vmr_apply(ivmr) ) cycle call RF_VMR_Write( RF_VMR(region,ivmr), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end do ! lt output: if ( lt_apply ) then call RF_LT_Write( RF_LT(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if if ( lt2_apply ) then call RF_LT_Write( RF_LT2(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if #ifdef with_budgets ! deposition velocities: if ( depv_apply ) then call RF_DEPV_Write( RF_DEPV(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if #endif status = 0 END SUBROUTINE PDUMP_FILES_WRITE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PDUMP_FILES_WRITE2 ! ! !DESCRIPTION: write at end of time interval ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE PDUMP_FILES_WRITE2( region, idate_f, status ) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/PDUMP_Files_Write2' ! --- begin ------------------------------- #ifdef with_budgets ! deposition fluxes: if ( deps_apply ) then call RF_DEPS_Write( RF_DEPS(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if #endif ! lt output: if ( lt_apply ) then call RF_LT_Write( RF_LT(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if if ( lt2_apply ) then call RF_LT_Write( RF_LT2(region), region, idate_f, status ) IF_NOTOK_RETURN(status=1) end if ! ok status = 0 END SUBROUTINE PDUMP_FILES_WRITE2 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: PDUMP_FILES_CLOSE ! ! !DESCRIPTION: call done method of each output file. !\\ !\\ ! !INTERFACE: ! SUBROUTINE PDUMP_FILES_CLOSE( region, status ) ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 31 Aug 2012 - Ph. Le Sager - switch closing order, since it was giving issues on some machine. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/PDUMP_Files_Close' ! --- local ------------------------------- integer :: ivmr ! --- begin ------------------------------- #ifdef with_budgets if ( depv_apply ) then call RF_DEPV_Done( RF_DEPV(region), status ) IF_NOTOK_RETURN(status=1) end if if ( deps_apply ) then call RF_DEPS_Done( RF_DEPS(region), status ) IF_NOTOK_RETURN(status=1) end if #endif if ( lt2_apply ) then call RF_LT_Done( RF_LT2(region), region, status ) IF_NOTOK_RETURN(status=1) end if if ( lt_apply ) then call RF_LT_Done( RF_LT(region), region, status ) IF_NOTOK_RETURN(status=1) end if do ivmr = nvmr, 1, -1 if ( .not. vmr_apply(ivmr) ) cycle call RF_VMR_Done( RF_VMR(region,ivmr), status ) IF_NOTOK_RETURN(status=1) end do if ( tp_apply ) then call RF_TP_Done ( RF_TP(region) , status ) IF_NOTOK_RETURN(status=1) end if if ( griddef_apply ) then call RF_GridDef_Done( RF_GridDef(region), status ) IF_NOTOK_RETURN(status=1) end if status = 0 end subroutine PDUMP_Files_Close !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_GRIDDEF_INIT ! ! !DESCRIPTION: ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE 1: Model horizontal grid definition ! (longitude, latitude, size of gridbox [m2] ). ! For documentation purposes, please also include the native vertical ! grid definition from your model (hybrid level coefficients) and the ! formula used to calculate pressure. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !\\ !\\ ! !INTERFACE: ! subroutine RF_GridDef_Init( RF, fdir, model, expid, region, status ) ! ! !USES: ! use partools, only : MPI_INFO_NULL, localComm use MeteoData, only : global_lli, levi ! ! !OUTPUT PARAMETERS: ! type(TPdumpFile_GridDef), intent(out) :: RF ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid integer, intent(in) :: region integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! 10 Jul 2012 - Ph. Le Sager - switch to MDF_NETCDF4 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_GridDef_Init' character(len=256) :: fname integer :: varid integer :: rtype ! --- begin ------------------------------------- call goLabel(rname) ! o open file ! write filename write (fname,'(a,"/",a,a,"_",a,"_",a,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'griddef' #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title', 'model horizontal definition' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version' , trim(dataset_version) , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call MDF_Def_Dim( RF%ncid, 'scalar', 1, RF%dimid_scalar , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lon', global_lli(region)%nlon, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lat', global_lli(region)%nlat, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lev', levi%nlev, RF%dimid_lev , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'levi', levi%nlev+1, RF%dimid_levi , status) IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Def_Dim( RF%ncid, 'time', NTS, RF%dimid_time , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Def_Dim( RF%ncid, 'datelen', 6, RF%dimid_datelen , status) !IF_NOTOK_MDF(fid=RF%ncid) ! o define variables rtype = MDF_FLOAT call MDF_Def_Var( RF%ncid, 'lon', rtype, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'degrees_east' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call MDF_Def_Var( RF%ncid, 'lat', MDF_FLOAT, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'degrees_north' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid !call MDF_Def_Var( RF%ncid, 'time', MDF_FLOAT, RF%dimid_time, varid , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'time' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'units', 'days since 1950-01-01 00:00:00' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'calender', 'gregorian' , status) !IF_NOTOK_MDF(fid=RF%ncid) !RF%varid_time = varid !call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'date' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status) !IF_NOTOK_MDF(fid=RF%ncid) !RF%varid_date = varid call MDF_Def_Var( RF%ncid, 'area', MDF_FLOAT, (/RF%dimid_lon,RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'grid_cell_area' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'grid-cell area' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'm2' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_gridbox_area = varid call MDF_Def_Var( RF%ncid, 'a', MDF_FLOAT, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a(k)*p0 + b(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'comment', 'bottom-up' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_a = varid call MDF_Def_Var( RF%ncid, 'b', mdf_float, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a(k)*p0 + b(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'comment', 'bottom-up' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_b = varid call MDF_Def_Var( RF%ncid, 'a_bnds', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient for layer bounds' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_a_bnds = varid call MDF_Def_Var( RF%ncid, 'b_bnds', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient for layer bounds' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_b_bnds = varid call MDF_Def_Var( RF%ncid, 'p0', mdf_float, (/RF%dimid_scalar/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'reference pressure value' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_p0 = varid !status = pnf90_def_var( RF%ncid, 'ps', MDF_FLOAT, & ! (/RF%dimid_lon,RF%dimid_lat,RF%dimid_time/), varid ) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'surface pressure' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status) !IF_NOTOK_MDF(fid=RF%ncid) !RF%varid_ps = varid !status = pnf90_def_var( RF%ncid, 'geo_height', MDF_FLOAT, & ! (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid ) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'geopotential height' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'units', 'm' , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Att( RF%ncid, varid, 'comment', 'bottom-up; lower half level; top value implicit infinity' , status) !IF_NOTOK_MDF(fid=RF%ncid) !RF%varid_geo_height = varid ! o end defintion mode call MDF_EndDef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! no records written yet RF%trec = 0 call goLabel() ; status = 0 END SUBROUTINE RF_GRIDDEF_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_GridDef_Write ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_GRIDDEF_WRITE( RF, region, status ) ! ! !USES: ! use GO, only : TDate, NewDate, rTotal, operator(-) use Grid, only : AreaOper use MeteoData, only : global_lli, levi, sp_dat ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_GridDef), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! 10 Jul 2012 - Ph. Le Sager - switch to MDF_NETCDF4 ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_GridDef_Write' integer :: imr, jmr, lmr real, allocatable :: ll(:,:) real :: time ! --- begin ------------------------------------- call goLabel(rname) ! grid size imr = global_lli(region)%nlon jmr = global_lli(region)%nlat lmr = levi%nlev ! next time record: RF%trec = RF%trec + 1 ! o write data if ( RF%trec == 1 ) then ! lat/lon field: allocate( ll(imr,jmr) ) call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg, status) IF_NOTOK_MDF(fid=RF%ncid) ll = 1.0 call AreaOper( global_lli(region), ll, '*', 'm2', status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Var( RF%ncid, RF%varid_gridbox_area, ll , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_a, levi%fa , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_b, levi%fb , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_a_bnds, levi%a(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_b_bnds, levi%b(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_p0, (/1.0/) , status) IF_NOTOK_MDF(fid=RF%ncid) deallocate( ll ) end if !call MDF_Put_Var( RF%ncid, RF%varid_time, time, index=RF%trec , status) !IF_NOTOK_MDF(fid=RF%ncid) !call MDF_Put_Var( RF%ncid, RF%varid_date, reshape(real(idate_f),(/6,1/), status), & ! start=(/1,RF%trec/), count=(/6,1/) ) !IF_NOTOK_MDF(fid=RF%ncid) !status = pnf90_put_var( RF%ncid, RF%varid_ps, & ! reshape(sp_dat(region)%data(1:imr,1:jmr,1:1),(/imr,jmr,1/)), & ! start=(/1,1,RF%trec/), count=(/imr,jmr,1/) ) !IF_NOTOK_MDF(fid=RF%ncid) !status = pnf90_put_var( RF%ncid, RF%varid_geo_height, & ! reshape(gph_dat(region)%data(1:imr,1:jmr,1:lmr),(/imr,jmr,lmr,1/)), & ! start=(/1,1,1,RF%trec/), count=(/imr,jmr,lmr,1/) ) !IF_NOTOK_MDF(fid=RF%ncid) call goLabel() status = 0 END SUBROUTINE RF_GridDef_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_GRIDDEF_DONE ! ! !DESCRIPTION: close file-1 !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_GridDef_Done( RF, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_GridDef), intent(inout) :: RF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_GridDef_Done' ! --- begin ------------------------------------- call goLabel(rname) call MDF_Close( RF%ncid , status) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 END SUBROUTINE RF_GRIDDEF_DONE !EOC ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE2: 3D field of monthly Model pressure [Pa] and temperature [K]. ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_TP_INIT ! ! !DESCRIPTION: file-2 : open and define var/att ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_TP_Init( RF, fdir, model, expid, region, idate_f, dhour, status ) ! ! !USES: ! use partools, only : MPI_INFO_NULL, localComm use MeteoData, only : global_lli, levi ! ! !OUTPUT PARAMETERS: ! type(TPdumpFile_TP), intent(out) :: RF integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid integer, intent(in) :: region integer, intent(in) :: idate_f(6) integer, intent(in) :: dhour ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - switch to netcdf-4 thru MDF ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_TP_Init' ! --- local ------------------------------------ character(len=256) :: fname integer :: varid, i1, i2, j1, j2 ! --- begin ------------------------------------- call goLabel(rname) ! store arguments RF%dhour = dhour call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) n_tp_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, mess='TP_Init' ) if ( n_tp_rec == 0 ) then tp_apply = .false. status=0 return end if ! o open file ! write filename write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'TP', idate_f(1:3) ! open, overwrite existing files (clobber) #ifdef MPI call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes call mdf_put_att( RF%ncid, MDF_GLOBAL, 'title', 'model pressure and temperature', status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author) , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution) , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, MDF_GLOBAL, 'dataset_version', trim(dataset_version) , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call mdf_def_dim( RF%ncid, 'lon', global_lli(region)%nlon, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_def_dim( RF%ncid, 'lat', global_lli(region)%nlat, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_def_dim( RF%ncid, 'lev', levi%nlev, RF%dimid_lev , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_def_dim( RF%ncid, 'time', n_tp_rec, RF%dimid_time , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_def_dim( RF%ncid, 'datelen', 6, RF%dimid_datelen , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define variables call mdf_def_var( RF%ncid, 'lon', MDF_FLOAT, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'degrees_east' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call mdf_def_var( RF%ncid, 'lat', MDF_FLOAT, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'degrees_north' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid call mdf_def_var( RF%ncid, 'lev', MDF_FLOAT, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'level' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lev = varid call mdf_def_var( RF%ncid, 'time', MDF_FLOAT, (/RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'days since 1950-01-01 00:00:00' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'calender', 'gregorian' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_time = varid allocate(RF%time(n_tp_rec)) call mdf_def_var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'date and time' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_date = varid allocate(RF%date(6,n_tp_rec)) call mdf_def_var( RF%ncid, 'ps', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_air_pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'surface pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_ps = varid call mdf_def_var( RF%ncid, 'orog', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_altitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'surface altitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'm' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_orog = varid call mdf_def_var( RF%ncid, 'surface_temp', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'surface temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'K' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', & '2m temperature from MARS archive or IFS model (grib 167, 2T)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_surface_temp = varid allocate( RF%data2d(i1:i2, j1:j2, n_tp_rec, 3) ) call mdf_def_var( RF%ncid, 'geopotential', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), & varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'geopotential' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'geopotential' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'm2 s-2' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'at mid levels' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_geop = varid call mdf_def_var( RF%ncid, 'pressure', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'at mid levels' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_pressure = varid call mdf_def_var( RF%ncid, 'temp', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'air_temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'K' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'bottom-up; full levels' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_temp = varid call mdf_def_var( RF%ncid, 'specific_humidity', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), & varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'specific_humidity' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'specific humidity' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'kg kg-1' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'mass fraction of water vapor in moist air; (kg water)/(kg air)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_humid = varid call mdf_def_var( RF%ncid, 'U', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'eastward_wind' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'zonal wind' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'm s-1' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_u = varid call mdf_def_var( RF%ncid, 'V', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'standard_name', 'northward_wind' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'meridional wind' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'm s-1' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_v = varid call mdf_def_var( RF%ncid, 'W', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'long_name', 'vertical wind velocity' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'units', 'm s-1' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_w = varid allocate( RF%data3d(i1:i2, j1:j2, levi%nlev, n_tp_rec, 7) ) ! o end defintion mode call mdf_enddef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! o ! no records written yet RF%trec = 0 call goLabel() ! ok status = 0 END SUBROUTINE RF_TP_Init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_TP_Write ! ! !DESCRIPTION: store records, and if last time step write data to file !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_TP_Write( RF, region, idate_f, status ) ! ! !USES: ! use Binas , only : grav use Phys , only : GeoPotentialHeight use Grid , only : FPressure, HPressure use GO , only : TDate, NewDate, rTotal, operator(-) use partools , only : myid, root use MeteoData , only : global_lli, lli, levi use MeteoData , only : sp_dat, temper_dat, humid_dat, pu_dat, pv_dat, mfw_dat, gph_dat, oro_dat, t2m_dat use MeteoData , only : m_dat use global_data, only : mass_dat ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_TP), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_TP_Write' ! --- local ------------------------------------ integer :: i, j, l, i1, i2, j1, j2 integer :: imr, jmr, lmr, klm real :: lev(levi%nlev) type(TDate) :: t, t0 real :: time real, allocatable :: field3d(:,:,:) real :: p_hlev(0:levi%nlev) ! --- begin ------------------------------------- ! for multiple of dhour only ... if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0) ) then status=0; return end if call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev ! next time record: RF%trec = RF%trec + 1 ! time since reftime: t0 = NewDate( time6=time_reftime6 ) t = NewDate( time6=idate_f ) time = rTotal( t - t0, 'day' ) if(okdebug)then write(gol,*) "RF_TP_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr end if ! o write data if ( RF%trec == 1 ) then call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status) IF_NOTOK_MDF(fid=RF%ncid) do l = 1, lmr lev(l) = real(l) end do call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status) IF_NOTOK_MDF(fid=RF%ncid) end if ! temporary storage for 3D fields allocate( field3d(i1:i2,j1:j2,1:lmr) ) ; field3d = 0. !-------- FILL DIAGNOSTIC ARRAYS RF%time(RF%trec) = time RF%date(:,RF%trec) = real(idate_f) RF%data2d(:,:,RF%trec,1) = sp_dat(region)%data(i1:i2,j1:j2,1) RF%data2d(:,:,RF%trec,2) = oro_dat(region)%data(i1:i2,j1:j2,1) RF%data2d(:,:,RF%trec,3) = t2m_dat(region)%data(i1:i2,j1:j2,1) ! o geopotential ! fill mid level geopotential: do j = j1, j2 do i = i1, i2 ! half level pressures call HPressure( levi, sp_dat(region)%data(i,j,1), p_hlev, status ) IF_NOTOK_RETURN(status=1) ! mid level gph (m) call GeoPotentialHeight( lmr, p_hlev, temper_dat(region)%data(i,j,:), & humid_dat(region)%data(i,j,:), oro_dat(region)%data(i,j,1)/grav, & field3d(i,j,:) ) ! m end do end do ! multiply with gravity for correct unit: field3d = field3d * grav ! m2/s2 RF%data3d(:,:,:,RF%trec,1) = field3d ! o pressure ! fill mid level pressure call FPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), field3d, status ) IF_NOTOK_RETURN(status=1) RF%data3d(:,:,:,RF%trec,2) = field3d ! o temperature RF%data3d(:,:,:,RF%trec,3) = temper_dat(region)%data(i1:i2,j1:j2,1:lmr) ! o specific humidity RF%data3d(:,:,:,RF%trec,4) = humid_dat(region)%data(i1:i2,j1:j2,1:lmr) ! o wind fields CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ! average U wind field3d = 0.5 * ( pu_dat(region)%data(i1-1:i2-1,j1:j2,1:lmr) + pu_dat(region)%data(i1:i2,j1:j2,1:lmr) ) & / m_dat(region)%data(i1:i2,j1:j2,1:lmr) ! 1/s do j = j1, j2 field3d(:,j,:) = field3d(:,j,:) * lli(region)%dx(j-j1+1) ! m/s end do RF%data3d(:,:,:,RF%trec,5) = field3d ! average V wind: field3d = 0.5 * ( pv_dat(region)%data(i1:i2,j1-1:j2-1,1:lmr) + pv_dat(region)%data(i1:i2,j1:j2,1:lmr) ) & / m_dat(region)%data(i1:i2,j1:j2,1:lmr) ! 1/s field3d = field3d * lli(region)%dy ! m/s RF%data3d(:,:,:,RF%trec,6) = field3d ! from downward massflux to upward average W wind: field3d = 0.5 * ( mfw_dat(region)%data(i1:i2,j1:j2,0:lmr-1) + mfw_dat(region)%data(i1:i2,j1:j2,1:lmr) ) & / m_dat(region)%data(i1:i2,j1:j2,1:lmr) ! 1/s do l = 1, lmr field3d(:,:,l) = - 1.0 * field3d(:,:,l) * & abs( gph_dat(region)%data(i1:i2,j1:j2,l+1) - gph_dat(region)%data(i1:i2,j1:j2,l) ) ! m/s end do RF%data3d(:,:,:,RF%trec,7) = field3d !-------- WRITE ARRAYS if ( RF%trec == n_tp_rec ) then ! time call MDF_Put_Var( RF%ncid, RF%varid_time, RF%time, status)!, start=(/1/), count=(/n_tp_rec/)) IF_NOTOK_MDF(fid=RF%ncid) ! date call MDF_Put_Var( RF%ncid, RF%varid_date, RF%date, status )!, & ! start=(/1,1/), count=(/6,1/) ) IF_NOTOK_MDF(fid=RF%ncid) ! surface pressure call MDF_Put_Var( RF%ncid, RF%varid_ps, RF%data2d(:,:,:,1), status, start=(/i1,j1,1/), count=(/imr,jmr,n_tp_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) ! orography (in m!) call MDF_Put_Var( RF%ncid, RF%varid_orog, RF%data2d(:,:,:,2), status, start=(/i1,j1,1/), count=(/imr,jmr,n_tp_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) ! surface temperature = 2m temperature call MDF_Put_Var( RF%ncid, RF%varid_surface_temp, RF%data2d(:,:,:,3), status, start=(/i1,j1,1/) ) !, count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=RF%ncid) ! geopotential call MDF_Put_Var( RF%ncid, RF%varid_geop, RF%data3d(:,:,:,:,1), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/)) IF_NOTOK_MDF(fid=RF%ncid) ! pressure call MDF_Put_Var( RF%ncid, RF%varid_pressure, RF%data3d(:,:,:,:,2), status, start=(/i1,j1,1,1/), & count=(/imr,jmr,lmr,n_tp_rec/)) IF_NOTOK_MDF(fid=RF%ncid) ! temperature call MDF_Put_Var( RF%ncid, RF%varid_temp, RF%data3d(:,:,:,:,3), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/)) IF_NOTOK_MDF(fid=RF%ncid) ! specific humidity call MDF_Put_Var( RF%ncid, RF%varid_humid, RF%data3d(:,:,:,:,4), status, start=(/i1,j1,1,1/),count=(/imr,jmr,lmr,n_tp_rec/)) IF_NOTOK_MDF(fid=RF%ncid) ! winds call MDF_Put_Var( RF%ncid, RF%varid_u, RF%data3d(:,:,:,:,5), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_v, RF%data3d(:,:,:,:,6), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_w, RF%data3d(:,:,:,:,7), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) end if ! Done deallocate( field3d ) call goLabel() status = 0 END SUBROUTINE RF_TP_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_TP_Done ! ! !DESCRIPTION: close file #2 !\\ !\\ ! !INTERFACE: ! subroutine RF_TP_Done( RF, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_TP), intent(inout) :: RF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_TP_Done' ! --- begin ------------------------------------- call goLabel(rname) call MDF_Close( RF%ncid , status) IF_NOTOK_RETURN(status=1) deallocate( rf%time, rf%date, rf%data2d, rf%data3d ) call goLabel() ; status = 0 end subroutine RF_TP_Done !EOC ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE3: 3D fields for O3, CO, CH4, ... Volume Mixing Ratios ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_VMR_Init ! ! !DESCRIPTION: open and define variables/attribute for file #3 !\\ !\\ ! !INTERFACE: ! subroutine RF_VMR_Init( RF, fdir, model, expid, filetype, region, & idate_f, dhour, tracer_names, status ) ! ! !USES: ! use Binas, only : xmair use GO, only : goReadFromLine, goUpCase use chem_param, only : ntrace, names, ra use partools, only : PAR_BROADCAST, MPI_INFO_NULL, localComm use MeteoData, only : global_lli, lli, levi, sp_dat use dims, only : xbeg, xend, ybeg, yend, dx, dy, dz, xref, yref, zref use dims, only : zbeg, zend ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_VMR), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid character(len=*), intent(in) :: filetype integer, intent(in) :: region integer, intent(in) :: idate_f(6) real, intent(in) :: dhour character(len=*), intent(in) :: tracer_names ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! 15 Apr 2014 - Ph. Le Sager - tropomi add-ons ! 8 Oct 2014 - H. Eskes - tropomi add-ons ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_VMR_Init' ! --- local ------------------------------------ character(len=256) :: fname, history, sysdate, model_meteo integer :: varid, i1, i2, j1, j2 integer, dimension(8) :: isysdate character(len=256) :: trnames character(len=8) :: trname, tmname integer :: k, itr, posend, pospoint integer :: imr, jmr, lmr, si, ei, ix, jy character(len=32) :: varname_spec character(len=5) :: zone character(len=64) :: cf_medium_stnd, cf_medium_long character(len=64) :: cf_enti_stnd, cf_enti_long, cf_enti_unit character(len=64) :: cf_spec_stnd, cf_spec_long character(len=4) :: cf_enti_type character(len=256) :: cf_name_stnd, cf_name_long, cf_name_unit character(len=512) :: comment character(len=6) :: csize integer, dimension(6) :: idate_f_end, idate_create type(TDate) :: date_f_end, date_create ! --- begin ------------------------------------- call goLabel(rname) ! store arguments RF%dhour = dhour RF%dsec = int(dhour*3600.) RF%tracer_names = tracer_names ! Test that dsec is multiple of dynamic-step/2 (nread in sec) if (((RF%dsec*2)/nread < 1).or.(modulo(RF%dsec,nread/2)/=0))then write(gol,*) "timeseries timestep should be a multiple of (dynamic_timestep)/2"; call goErr TRACEBACK; status=1; return end if ! size imr = global_lli(region)%nlon jmr = global_lli(region)%nlat lmr = levi%nlev ! number of time steps rf%n_rec = GET_N_TIME_RECORDS( idate_f, rf%dsec, mess='VMR_Init' ) ! degenerated cases (eg, very short runs) if ( rf%n_rec == 0 ) then rf%apply = .false. status=0 return end if call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! set tracer index for requested tracers: write (gol,'("selected tracers for VMR output:")'); call goPr ! initialise RF RF%ntr = 0 #ifdef with_m7 RF%lpmx = .false. RF%sizepmx = -1.0 #endif RF%itr = -1 trnames = tracer_names do ! empty ? if ( len_trim(trnames) == 0 ) exit ! next number: if ( RF%ntr == ntrace ) then write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr TRACEBACK; status=1; return end if RF%ntr = RF%ntr + 1 ! extract leading name: call goReadFromLine( trnames, trname, status, sep=' ' ) IF_NOTOK_RETURN(status=1) #ifdef with_m7 ! --------------------------- ! check for PMx ! --------------------------- if( strlowercase(trname(1:2)) == 'pm' ) then RF%lpmx(RF%ntr) = .true. RF%itr (RF%ntr) = -1 ! paste size to real read(trname(3:len_trim(trname)), * ) RF%sizepmx(RF%ntr) else #endif ! convert to tm5 name: select case ( trim(strlowercase(trname)) ) case ( 'hcho' ) ; tmname = 'CH2O' case ( 'rn', 'radon' ) ; tmname = 'Rn222' case ( 'pb', 'lead' ) ; tmname = 'Pb210' case default ; tmname = trname end select ! -------------------------------- ! NOy and M7 are special cases ... ! -------------------------------- select case ( trim(strlowercase(tmname)) ) case( 'noy' ) ! defined as ntrace+1 RF%itr(RF%ntr) = iNOy write (gol,'(" * ",a10)') trim(trname); call goPr #ifdef with_m7 case( 'tso4' ) ! defined as ntrace+2 RF%itr(RF%ntr) = iSO4 write (gol,'(" * ",a10)') trim(trname); call goPr case( 'tbc' ) ! defined as ntrace+3 RF%itr(RF%ntr) = iBC write (gol,'(" * ",a10)') trim(trname); call goPr case( 'tpom' ) ! defined as ntrace+4 RF%itr(RF%ntr) = iPOM write (gol,'(" * ",a10)') trim(trname); call goPr case( 'tss' ) ! defined as ntrace+5 RF%itr(RF%ntr) = iSS write (gol,'(" * ",a10)') trim(trname); call goPr case( 'tdu' ) ! defined as ntrace+6 RF%itr(RF%ntr) = iDU write (gol,'(" * ",a10)') trim(trname); call goPr #endif case default ! -------------------------------- ! `regular` constituents ! -------------------------------- ! loop over all names: RF%itr(RF%ntr) = -1 do itr = 1, ntrace ! case indendent match ? if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then write (gol,'(" ",i3," ",a10," (",a10,") ",f12.4)') itr, trim(trname), trim(names(itr)), ra(itr); call goPr RF%itr(RF%ntr) = itr exit end if end do end select ! not found ? if ( RF%itr(RF%ntr) < 0 ) then write (gol,'("tracer name not supported:")'); call goPr write (gol,'(" list all : ",a)') trim(tracer_names); call goPr write (gol,'(" list element : ",i3)') RF%ntr; call goPr write (gol,'(" pdump name : ",a)') trim(trname); call goPr write (gol,'(" tm5 name : ",a)') trim(tmname); call goPr write (gol,'(" tm5 tracers : ")'); call goPr do itr = 1, ntrace write (gol,'(" ",i3," ",a)') itr, trim(names(itr)); call goPr end do TRACEBACK; status=1; return end if ! RF%itr #ifdef with_m7 end if ! pmx #endif ! store pdump name: RF%name_tr(RF%ntr) = tmname end do ! empty file ? if ( RF%ntr < 1 ) then write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr TRACEBACK; status=1; return end if ! o open file ! write filename #ifdef tropomi ! define start/stop of output, and run date date_f_end = NewDate( time6=idate_f ) + IncrDate(hour=24) call Get( date_f_end, time6=idate_f_end ) if ( isRoot ) then date_create = SystemDate() call Get( date_create, time6=idate_create ) endif call PAR_BROADCAST(idate_create, status) IF_NOTOK_RETURN(status=1) date_create = SystemDate() call Get( date_create, time6=idate_create ) write (tropomi_date_start, '(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_f write (tropomi_date_stop, '(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_f_end write (tropomi_date_create,'(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_create ! write filename according to tropomi convention write (fname,'(a,"/",a,"_",a,"_",a,".nc")') & trim(fdir), trim(tropomi_dataset_name), tropomi_date_start, tropomi_date_stop #else write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3) #endif ! open: #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes #ifdef tropomi ! H. Eskes: Extra attributes for TROPOMI ! Conventions = "CF-1.6" call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'Conventions', 'CF-1.6' , status) IF_NOTOK_MDF(fid=RF%ncid) ! validity_start = "20132305T120000" (zoals in filenaam) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'validity_start', tropomi_date_start , status) IF_NOTOK_MDF(fid=RF%ncid) ! validity_stop = "20132405T000000" call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'validity_stop', tropomi_date_stop , status) IF_NOTOK_MDF(fid=RF%ncid) ! creation_date = "20142909T124905" call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'creation_date', tropomi_date_create , status) IF_NOTOK_MDF(fid=RF%ncid) ! version = TM5 version string. call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'version', trim(tropomi_tm5_version) , status) IF_NOTOK_MDF(fid=RF%ncid) ! institution = "KNMI" call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution', trim(tropomi_institution) , status) IF_NOTOK_MDF(fid=RF%ncid) ! reference = TM5 reference (journal article or so) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'reference', trim(tropomi_tm5_reference) , status) IF_NOTOK_MDF(fid=RF%ncid) ! contact = email address of volunteer. call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'contact', trim(tropomi_authoremail) , status) IF_NOTOK_MDF(fid=RF%ncid) ! dataset_name = "S5P_NRTI_AUX_CTMFCT" of "S5P_OFFL_AUX_CTMANA" call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_name', trim(tropomi_dataset_name) , status) IF_NOTOK_MDF(fid=RF%ncid) #endif call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title' , 'mixing ratios & concentrations' , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version' , trim(dataset_version) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'file_version_number', trim(outfileversnr) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'im' , imr , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'jm' , jmr , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'lm' , lmr , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dx' , dx/xref(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dy' , dy/yref(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dz' , dz/zref(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'xbeg' , xbeg(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'xend' , xend(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'ybeg' , ybeg(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'yend' , yend(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'zbeg' , zbeg(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'zend' , zend(region) , status ) IF_NOTOK_MDF(fid=RF%ncid) ! Meteo attribute if (trim(meteo_class)=='ei') then model_meteo='analysis (ERA-Interim)' elseif (trim(meteo_class)=='ea') then model_meteo='reanalysis (ERA5)' elseif (trim(meteo_class)=='od') then model_meteo='forecast (IFS)' elseif (trim(meteo_class)=='ifs10') then model_meteo='EC-Earth (ifs 10L)' elseif (trim(meteo_class)=='ifs34') then model_meteo='EC-Earth (ifs 34L)' elseif (trim(meteo_class)=='ifs62') then model_meteo='EC-Earth (ifs 62L)' elseif (trim(meteo_class)=='ifs91') then model_meteo='EC-Earth (ifs 91L)' else write (gol,'("Meteo Model not known !")'); call goErr TRACEBACK; status=1; return endif call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'meteo_model', trim(model_meteo), status ) IF_NOTOK_MDF(fid=RF%ncid) ! History attribute for audit trail: date, time of day, user name, program name call date_and_time(values=isysdate, zone=zone) write (sysdate, '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2," ",a)') & isysdate(1), isysdate(2), isysdate(3), isysdate(5), isysdate(6), isysdate(7), zone write(history,'("Created ",a," by ",a," with TM5.")') trim(sysdate),trim(dataset_author) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'history', trim(history), status ) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call MDF_Def_Dim( RF%ncid, 'lon', imr, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lat', jmr, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lev', levi%nlev, RF%dimid_lev , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'levi', levi%nlev+1, RF%dimid_levi , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'time', rf%n_rec, RF%dimid_time , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'datelen', 6, RF%dimid_datelen , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define variables call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'degrees_east' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'degrees_north' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid #ifdef tropomi call MDF_Def_Var( RF%ncid, 'hyai', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid A coefficient at layer interfaces' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_hyai = varid #else call MDF_Def_Var( RF%ncid, 'a_bnds', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient for layer bounds' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_a_bnds = varid #endif #ifdef tropomi call MDF_Def_Var( RF%ncid, 'hybi', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid B coefficient at layer interfaces' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_hybi = varid #else call MDF_Def_Var( RF%ncid, 'b_bnds', mdf_float, (/RF%dimid_levi/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid sigma coordinate a coefficient for layer bounds' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_b_bnds = varid #endif #ifdef tropomi call MDF_Def_Var( RF%ncid, 'hyam', mdf_float, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid A coefficient at layer midpoints' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_hyam = varid call MDF_Def_Var( RF%ncid, 'hybm', mdf_float, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid B coefficient at layer midpoints' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_hybm = varid #endif call MDF_Def_Var( RF%ncid, 'lev', mdf_float, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) #ifdef tropomi call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid level at layer midpoints' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'level' , status) IF_NOTOK_MDF(fid=RF%ncid) call mdf_put_att( RF%ncid, varid, 'positive', 'down' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula', 'hyam hybm (mlev=hyam+hybm*ps)' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'ap: hyam b: hybm ps: ps' , status) IF_NOTOK_MDF(fid=RF%ncid) #else call MDF_Put_Att( RF%ncid, varid, 'long_name', 'level' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) #endif RF%varid_lev = varid call MDF_Def_Var( RF%ncid, 'time', mdf_double, (/RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'days since 1950-01-01 00:00:00' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'calender', 'gregorian' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_time = varid allocate(RF%time(rf%n_rec)) call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_date = varid allocate(RF%date(6,rf%n_rec)) call MDF_Def_Var( RF%ncid, 'ps', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'surface_air_pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'surface pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_ps = varid allocate( RF%sp(i1:i2, j1:j2, rf%n_rec) ) #ifndef tropomi call MDF_Def_Var( RF%ncid, 't', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'air_temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'K' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_put_att( RF%ncid, varid, 'comment', 'bottom-up; full levels' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_temp = varid allocate( RF%data3d_t(i1:i2, j1:j2, levi%nlev, rf%n_rec) ) #endif #ifdef tropomi ! Extra temperature field output ! with compression - crash !call MDF_Def_Var( RF%ncid, 't', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status, compression=1, deflate_level=4) call MDF_Def_Var( RF%ncid, 't', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'air_temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'temperature' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'K' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_put_att( RF%ncid, varid, 'comment', 'bottom-up; full levels' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_temp = varid allocate( RF%data3d_t(i1:i2, j1:j2, levi%nlev, rf%n_rec) ) ! Extra surface elevation output, retrieved from GPH (meteo.f90) and g0 (binas.f90) following WGS84? call MDF_Def_Var( RF%ncid, 'surface_altitude', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'surface_altitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'surface altitude of TM5 grid' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'm' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_put_att( RF%ncid, varid, 'comment', 'ECMWF interpolated orography' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_hgt = varid allocate( RF%data2d_hgt(i1:i2, j1:j2) ) ! Extra tropopause level output, retrieved from GPH and temperature (meteo.f90) call MDF_Def_Var( RF%ncid, 'tropopause_layer_index', MDF_INT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'tropopause_layer_index' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'index of the highest model layer in the troposphere' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', '-' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_put_att( RF%ncid, varid, 'comment', 'Based on WMO temperature gradient method' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_ltropo = varid allocate( RF%data2d_ltropo(i1:i2, j1:j2, rf%n_rec) ) #endif ! loop over tracer to be written: do k = 1, RF%ntr #ifdef with_m7 if( RF%lpmx(k) ) then ! get diameter write(csize,'(F5.1)') RF%sizepmx(k) ! remove leading blanks csize = adjustl(csize) pospoint = index(csize,'.') posend = len_trim(csize) ! CF standard name for concentration/mixing ratio/column: RF%varid_type(k) = 'conc' varname_spec = 'pm'//csize(1:pospoint-1)//'p'//csize(pospoint+1:posend) cf_spec_stnd = 'particulate_matter_'//trim(csize) cf_spec_long = 'particulate matter diameter le '//trim(csize)//' micrometers' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' else #endif ! ---------------------------- ! setting defaults (gas phase) ! ---------------------------- ! CF standard name for concentration/mixing ratio/column: cf_enti_stnd = 'mole_fraction' #ifdef tropomi cf_enti_unit = '1' #else cf_enti_unit = 'mole mole-1' #endif cf_enti_long = 'volume mixing ratio' cf_medium_stnd = 'in_air' cf_medium_long = 'in humid air' RF%varid_type(k) = 'mixr' ! global tracer index itr = RF%itr(k) ! no comment yet comment = '' ! standard names from CF conventions: select case ( strlowercase(RF%name_tr(k)) ) case ( 'co2' ) varname_spec = 'co2' cf_spec_stnd = 'carbon_dioxide' cf_spec_long = 'CO2' case ( 'co' ) varname_spec = 'co' cf_spec_stnd = 'carbon_monoxide' cf_spec_long = 'CO' case ( 'o3' ) varname_spec = 'o3' cf_spec_stnd = 'ozone' cf_spec_long = 'O3' case ( 'o3s' ) varname_spec = 'o3s' cf_spec_stnd = 'ozone_from_stratosphere' cf_spec_long = 'O3s' case ( 'no' ) varname_spec = 'no' cf_spec_stnd = 'nitrogen_monoxide' cf_spec_long = 'NO' case ( 'no2' ) varname_spec = 'no2' cf_spec_stnd = 'nitrogen_dioxide' cf_spec_long = 'NO2' case ( 'noy' ) varname_spec = 'noy' cf_spec_stnd = 'nitrogen_oxides' cf_spec_long = 'NOy' comment = 'NOy = NOx + HNO3 + PAN + org.ntr., '// & 'with NOx = NO + NO2 + NO3 + HNO4 + N2O5' case ( 'ch2o', 'choh' ) varname_spec = 'ch2o' cf_spec_stnd = 'formaldehyde' cf_spec_long = 'CH2O' case ( 'so2' ) varname_spec = 'so2' cf_spec_stnd = 'sulfur_dioxide' cf_spec_long = 'SO2' case( 'h2so4' ) varname_spec = 'h2so4' cf_spec_stnd = 'sulfuric_acid_g' cf_spec_long = 'H2SO4 (g)' !!$ case ( 'so4' ) !!$ varname_spec = 'so4' !!$ cf_spec_stnd = 'sulfate_as_sulfate_dry_aerosol' !!$ cf_spec_long = 'SO4' case ( 'ch4' ) varname_spec = 'ch4' cf_spec_stnd = 'methane' cf_spec_long = 'CH4' case ( 'oh' ) varname_spec = 'oh' cf_spec_stnd = 'hydroxyl_radical' cf_spec_long = 'OH' case ( 'h2o2' ) varname_spec = 'h2o2' cf_spec_stnd = 'hydrogen_peroxide' cf_spec_long = 'H2O2' case ( 'hno3' ) varname_spec = 'hno3' cf_spec_stnd = 'nitric_acid' cf_spec_long = 'HNO3' case ( 'hno4' ) varname_spec = 'hno4' cf_spec_stnd = 'peroxonitric_acid' cf_spec_long = 'HNO4' case ( 'n2o5' ) varname_spec = 'n2o5' cf_spec_stnd = 'nitrogen_pentoxide' cf_spec_long = 'N2O5' case ( 'par' ) varname_spec = 'par' cf_spec_stnd = 'paraffinic_carbon_atoms' cf_spec_long = 'PAR' case ( 'eth' ) varname_spec = 'eth' cf_spec_stnd = 'ethylene' cf_spec_long = 'ETH' case ( 'ole' ) varname_spec = 'ole' cf_spec_stnd = 'olefinic_carbon_bonds' cf_spec_long = 'OLE' case ( 'ald2' ) varname_spec = 'ald2' cf_spec_stnd = 'acetaldehyde_and_higher_aldehydes' cf_spec_long = 'ALD2' case ( 'mgly' ) varname_spec = 'mgly' cf_spec_stnd = 'methylglyoxal' cf_spec_long = 'MGLY' case ( 'isop' ) varname_spec = 'isop' cf_spec_stnd = 'isoprene' cf_spec_long = 'ISOP' case ( 'nh3' ) varname_spec = 'nh3' cf_spec_stnd = 'ammonia' cf_spec_long = 'NH3' case ( 'ORGNTR','orgntr' ) varname_spec = 'orgntr' cf_spec_stnd = 'organic_nitrate' cf_spec_long = 'ORGNTR' case ( 'pan' ) varname_spec = 'pan' cf_spec_stnd = 'peroxyacetyl_nitrate' cf_spec_long = 'PAN' case ( 'terp' ) varname_spec = 'terp' cf_spec_stnd = 'terpene' cf_spec_long = 'TERP' case ( 'elvoc' ) varname_spec = 'elvoc' cf_spec_stnd = 'extremely low volatile OC' cf_spec_long = 'ELVOC' case ( 'svoc' ) varname_spec = 'svoc' cf_spec_stnd = 'semi volatile OC' cf_spec_long = 'SVOC' case ( 'rn', 'radon', 'rn222' ) varname_spec = 'rn' cf_spec_stnd = 'radon' cf_spec_long = 'Rn' case ( 'pb', 'lead', 'pb210' ) varname_spec = 'pb' cf_spec_stnd = 'lead' cf_spec_long = 'Pb' #ifdef with_m7 ! Sulphate case( 'tso4' ) RF%varid_type(k) = 'conc' varname_spec = 'so4' cf_spec_stnd = 'total_sulphate_aerosol' cf_spec_long = 'SO4' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Black Carbon case( 'tbc' ) RF%varid_type(k) = 'conc' varname_spec = 'bc' cf_spec_stnd = 'total_black_carbon_aerosol' cf_spec_long = 'BC' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Particulate Organic Matter case( 'tpom' ) RF%varid_type(k) = 'conc' varname_spec = 'pom' cf_spec_stnd = 'total_particulate_organic_matter_aerosol' cf_spec_long = 'POM' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Sea Salt case( 'tss' ) RF%varid_type(k) = 'conc' varname_spec = 'ss' cf_spec_stnd = 'total_sea_salt_aerosol' cf_spec_long = 'SS' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Dust case( 'tdu' ) RF%varid_type(k) = 'conc' varname_spec = 'du' cf_spec_stnd = 'total_dust_aerosol' cf_spec_long = 'SS' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Nucleation Soluble (nus): number, SO4 case ( 'nus_n' ) RF%varid_type(k) = 'numb' varname_spec = 'nus_n' cf_spec_stnd = 'number_wet_nucleation' cf_spec_long = 'Number_nus' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'so4nus' ) RF%varid_type(k) = 'conc' varname_spec = 'so4nus' cf_spec_stnd = 'sulphate_wet_nucleation' cf_spec_long = 'SO4_nus' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'soanus' ) RF%varid_type(k) = 'conc' varname_spec = 'soanus' cf_spec_stnd = 'SOA_wet_nucleation' cf_spec_long = 'SOA_nus' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Aitken Soluble (ais): number, SO4, BC, POM case ( 'ais_n' ) RF%varid_type(k) = 'numb' varname_spec = 'ais_n' cf_spec_stnd = 'number_wet_aitken' cf_spec_long = 'Number_ais' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'so4ais' ) RF%varid_type(k) = 'conc' varname_spec = 'so4ais' cf_spec_stnd = 'sulphate_wet_aitken' cf_spec_long = 'SO4_ais' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'bcais' ) RF%varid_type(k) = 'conc' varname_spec = 'bcais' cf_spec_stnd = 'black_carbon_wet_aitken' cf_spec_long = 'BC_ais' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'pomais' ) RF%varid_type(k) = 'conc' varname_spec = 'pomais' cf_spec_stnd = 'particulate_organic_matter_wet_aitken' cf_spec_long = 'POM_ais' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'soaais' ) RF%varid_type(k) = 'conc' varname_spec = 'soaais' cf_spec_stnd = 'SOA_dry_Aitken' cf_spec_long = 'SOA_ais' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Accumulation Soluble (acs): number, SO4, BC, POM, SS, DU case ( 'acs_n' ) RF%varid_type(k) = 'numb' varname_spec = 'acs_n' cf_spec_stnd = 'number_wet_accumulation' cf_spec_long = 'Number_acs' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'so4acs' ) RF%varid_type(k) = 'conc' varname_spec = 'so4acs' cf_spec_stnd = 'sulphate_wet_accumulation' cf_spec_long = 'SO4_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'bcacs' ) RF%varid_type(k) = 'conc' varname_spec = 'bcacs' cf_spec_stnd = 'black_carbon_wet_accumulation' cf_spec_long = 'BC_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'pomacs' ) RF%varid_type(k) = 'conc' varname_spec = 'pomacs' cf_spec_stnd = 'particulate_organic_matter_wet_accumulation' cf_spec_long = 'POM_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'ssacs' ) RF%varid_type(k) = 'conc' varname_spec = 'ssacs' cf_spec_stnd = 'seasalt_wet_accumulation' cf_spec_long = 'SS_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'duacs' ) RF%varid_type(k) = 'conc' varname_spec = 'duacs' cf_spec_stnd = 'dust_wet_accumulation' cf_spec_long = 'DU_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'soaacs' ) RF%varid_type(k) = 'conc' varname_spec = 'soaacs' cf_spec_stnd = 'SOA_dry_Accumulation' cf_spec_long = 'SOA_acs' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Coarse Soluble (cos): number, SO4, BC, POM, SS, DU case ( 'cos_n' ) RF%varid_type(k) = 'numb' varname_spec = 'cos_n' cf_spec_stnd = 'number_wet_coarse' cf_spec_long = 'Number_cos' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'so4cos' ) RF%varid_type(k) = 'conc' varname_spec = 'so4cos' cf_spec_stnd = 'sulphate_wet_coarse' cf_spec_long = 'SO4_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'bccos' ) RF%varid_type(k) = 'conc' varname_spec = 'bccos' cf_spec_stnd = 'black_carbon_wet_coarse' cf_spec_long = 'BC_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'pomcos' ) RF%varid_type(k) = 'conc' varname_spec = 'pomcos' cf_spec_stnd = 'particulate_organic_matter_wet_coarse' cf_spec_long = 'POM_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'sscos' ) RF%varid_type(k) = 'conc' varname_spec = 'sscos' cf_spec_stnd = 'seasalt_wet_coarse' cf_spec_long = 'SS_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'ducos' ) RF%varid_type(k) = 'conc' varname_spec = 'ducos' cf_spec_stnd = 'dust_wet_coarse' cf_spec_long = 'DU_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'soacos' ) RF%varid_type(k) = 'conc' varname_spec = 'soacos' cf_spec_stnd = 'SOA_dry_coarse' cf_spec_long = 'SOA_cos' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Aitken Insoluble (aii): number, BC, POM case ( 'aii_n' ) RF%varid_type(k) = 'numb' varname_spec = 'aii_n' cf_spec_stnd = 'number_dry_aitken' cf_spec_long = 'Number_aii' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'bcaii' ) RF%varid_type(k) = 'conc' varname_spec = 'bcaii' cf_spec_stnd = 'black_carbon_dry_aitken' cf_spec_long = 'BC_aii' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'pomaii' ) RF%varid_type(k) = 'conc' varname_spec = 'pomaii' cf_spec_stnd = 'particulate_organic_matter_dry_aitken' cf_spec_long = 'POM_aii' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'soaaii' ) RF%varid_type(k) = 'conc' varname_spec = 'soaaii' cf_spec_stnd = 'SOA_dry_Aitken' cf_spec_long = 'SOA_aii' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Accumulation Insoluble (aci): number, DU case ( 'aci_n' ) RF%varid_type(k) = 'numb' varname_spec = 'aci_n' cf_spec_stnd = 'number_dry_accumulation' cf_spec_long = 'Number_aci' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'duaci' ) RF%varid_type(k) = 'conc' varname_spec = 'duaci' cf_spec_stnd = 'dust_dry_accumulation' cf_spec_long = 'DU_aci' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' ! Coarse Insoluble (coi): number, DU case ( 'coi_n' ) RF%varid_type(k) = 'numb' varname_spec = 'coi_n' cf_spec_stnd = 'number_dry_coarse' cf_spec_long = 'Number_coi' cf_enti_stnd = 'number' cf_enti_unit = '1.' cf_enti_long = '' case ( 'ducoi' ) RF%varid_type(k) = 'conc' varname_spec = 'ducoi' cf_spec_stnd = 'dust_dry_coarse' cf_spec_long = 'DU_coi' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' #endif case ( 'nh4' ) RF%varid_type(k) = 'conc' varname_spec = 'nh4' cf_spec_stnd = 'ammonium_as_ammonium_dry_aerosol' cf_spec_long = 'NH4' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' case ( 'no3_a' ) RF%varid_type(k) = 'conc' varname_spec = 'no3' cf_spec_stnd = 'nitrate_as_nitrate_dry_aerosol' cf_spec_long = 'NO3' cf_enti_stnd = 'concentration' cf_enti_unit = 'kg m-3 ' cf_enti_long = 'mass per volume' !!$ case ( 'bc' ) !!$ varname_spec = 'bc' !!$ cf_spec_stnd = 'black_carbon_dry_aerosol' !!$ cf_spec_long = 'BC' !!$ case ( 'BCS', 'bcs' ) !!$ varname_spec = 'bcs' !!$ cf_spec_stnd = 'hydrophilic_black_carbon_dry_aerosol' !!$ cf_spec_long = 'BC(aq)' !!$ case ( 'POM', 'pom' ) !!$ varname_spec = 'om' !!$ cf_spec_stnd = 'organic_carbon_as_particulate_organic_matter_dry_aerosol' !!$ cf_spec_long = 'OM' !!$ case ( 'SS1_N', 'ss1_n' ) !!$ varname_spec = 'ss1_n' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode1_number' !!$ cf_spec_long = 'SS1_n' !!$ case ( 'SS1_M', 'ss1_m' ) !!$ varname_spec = 'ss1_m' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode1_mass' !!$ cf_spec_long = 'SS1_m' !!$ case ( 'SS2_N', 'ss2_n' ) !!$ varname_spec = 'ss2_n' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode2_number' !!$ cf_spec_long = 'SS2_n' !!$ case ( 'SS2_M', 'ss2_m' ) !!$ varname_spec = 'ss2_m' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode2_mass' !!$ cf_spec_long = 'SS2_m' !!$ case ( 'SS3_N', 'ss3_n' ) !!$ varname_spec = 'ss3_n' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode3_number' !!$ cf_spec_long = 'SS3_n' !!$ case ( 'SS3_M', 'ss3_m' ) !!$ varname_spec = 'ss3_m' !!$ cf_spec_stnd = 'seasalt_dry_aerosol_mode3_mass' !!$ cf_spec_long = 'SS3_m' !!$ case ( 'DUST2_N', 'dust2_n' ) !!$ varname_spec = 'dust2_n' !!$ cf_spec_stnd = 'dust_dry_aerosol_mode2_number' !!$ cf_spec_long = 'DUST2_n' !!$ case ( 'DUST2_M', 'dust2_m' ) !!$ varname_spec = 'dust2_m' !!$ cf_spec_stnd = 'dust_dry_aerosol_mode2_madust' !!$ cf_spec_long = 'DUST2_m' !!$ case ( 'DUST3_N', 'dust3_n' ) !!$ varname_spec = 'dust3_n' !!$ cf_spec_stnd = 'dust_dry_aerosol_mode3_number' !!$ cf_spec_long = 'DUST3_n' !!$ case ( 'DUST3_M', 'dust3_m' ) !!$ varname_spec = 'dust3_m' !!$ cf_spec_stnd = 'dust_dry_aerosol_mode3_madust' !!$ cf_spec_long = 'DUST3_m' case default write (gol,'("do not know how to match tracer with CF standard names : ",a)') RF%name_tr(k); call goErr TRACEBACK; status=1; return end select #ifdef with_m7 end if ! RF%lpmx(k) #endif ! define variable: call MDF_Def_Var( RF%ncid, trim(varname_spec), MDF_FLOAT, & (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) ! total names: cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)//'_'//trim(cf_medium_stnd) cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long)//' '//trim(cf_medium_long) cf_name_unit = trim(cf_enti_unit) ! write attributes: call MDF_Put_Att( RF%ncid, varid, 'standard_name', trim(cf_name_stnd) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', trim(cf_name_long) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', trim(cf_name_unit) , status) IF_NOTOK_MDF(fid=RF%ncid) ! moleweights; ra from chem_param is in g/mol . if ( itr <= ntrace .and. itr > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) else call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0 , status) IF_NOTOK_MDF(fid=RF%ncid) end if call MDF_Put_Att( RF%ncid , varid, 'moleweight_air' , xmair*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_unit' , 'kg mole-1' , status) IF_NOTOK_MDF(fid=RF%ncid) if ( len_trim(comment) > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'comment' , trim(comment), status) IF_NOTOK_MDF(fid=RF%ncid) end if ! store varid RF%varid_tr(k) = varid end do ! storage allocate(rf%data3d(i1:i2,j1:j2,lmr,rf%n_rec,rf%ntr)) ! o end defintion mode call MDF_EndDef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! o ! no records written yet RF%trec = 0 call goLabel() status = 0 END SUBROUTINE RF_VMR_Init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_VMR_Write ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_VMR_Write( RF, region, idate_f, status ) ! ! !USES: ! use Binas, only : xmair use GO, only : TDate, NewDate, rTotal, operator(-) use binas, only : Rgas use chem_param, only : ntrace, ntracet, fscale, ra use tracer_data, only : mass_dat, chem_dat use Grid, only : FPressure use MeteoData, only : global_lli, levi, m_dat, sp_dat, temper_dat #ifdef tropomi use MeteoData, only : gph_dat use toolbox, only : ltropo_ifs, lvlpress #endif #ifdef with_m7 use calc_pm, only : PMx_Integrate_3d #endif ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_VMR), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! 2 Oct 2012 - Ph. Le Sager - adapted for lat-lon mpi decomp ! - no more sub-regions available ! ! !REMARKS: ! (1) ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_VMR_Write' ! --- local ------------------------------------ integer :: imr, jmr, lmr, i1, i2, j1, j2, i, j real, allocatable :: lev(:) integer :: l type(TDate) :: t, t0 real :: time integer :: k, itr, dsec integer :: k_comp, itr_comp integer :: ims, ime, jms, jme, lms, lme integer :: gimr, gjmr, glmr real, allocatable :: compo_k(:,:,:) real, allocatable :: field_k(:,:,:) real, allocatable :: pres3d(:,:,:), pmx(:,:,:) integer :: numtrac integer :: listtrac(10) ! --- begin ------------------------------------- ! for multiple of timestep only ... dsec = idate_f(4)*3600 + idate_f(5)*60 + idate_f(6) if ( modulo(dsec,RF%dsec) /= 0 ) then status=0; return end if call goLabel(rname) ! grid sizes call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev gimr = global_lli(region)%nlon gjmr = global_lli(region)%nlat ! yet to change ?? lms = 1 lme = levi%nlev lmr = levi%nlev glmr = levi%nlev #ifdef with_m7 ! get helping pressure field in 3d allocate( pres3d(i1:i2,j1:j2,lmr) ) ! fill mid level pressure call FPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), pres3d, status ) IF_NOTOK_RETURN(status=1) #endif ! next time record: RF%trec = RF%trec + 1 if(isRoot.and.okdebug)then write(gol,*) "RF_VMR_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr end if ! time since 1950-1-1 00:00 t0 = NewDate( time6=time_reftime6 ) t = NewDate( time6=idate_f ) time = rTotal( t - t0, 'day' ) ! only once ... if ( RF%trec == 1 ) then ! write longitudes: call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! write latitudes: call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! write level indices: allocate( lev(lmr) ) do l = lms, lme lev(l) = real(l) end do call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status) IF_NOTOK_MDF(fid=RF%ncid) deallocate(lev) #ifdef tropomi ! As and Bs interfaces call MDF_Put_Var( RF%ncid, RF%varid_hyai, levi%a(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_hybi, levi%b(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) ! As and Bs mid-level (full level) call MDF_Put_Var( RF%ncid, RF%varid_hyam, levi%fa(1:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_hybm, levi%fb(1:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) #else ! As and Bs call MDF_Put_Var( RF%ncid, RF%varid_a_bnds, levi%a(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_b_bnds, levi%b(0:levi%nlev) , status) IF_NOTOK_MDF(fid=RF%ncid) #endif end if ! first record RF%time(RF%trec) = time RF%date(:,RF%trec) = real(idate_f) RF%sp(:,:,RF%trec) = sp_dat(region)%data(i1:i2,j1:j2,1) #ifdef tropomi ! copy of temperature field RF%data3d_t(:,:,:,RF%trec) = temper_dat(region)%data(i1:i2,j1:j2,1:lmr) ! orography: copy of lowest interface gph field. gph in the model is in "m", at interfaces, and gph(1)=oro ! only once ... if ( RF%trec == 1 ) then RF%data2d_hgt(:,:) = gph_dat(region)%data(i1:i2,j1:j2,1) end if ! compute highest tropopause layer index do i = i1, i2 do j = j1, j2 RF%data2d_ltropo(i,j,RF%trec) = ltropo_ifs(region,i,j,temper_dat(region)%data(i,j,1:lmr),lmr) end do end do #endif ! loop over all tracers to be written: do k = 1, RF%ntr ! global tracer index: itr = RF%itr(k) #ifdef with_m7 ! --------------------- ! particulate matter ! --------------------- if( RF%lpmx(k) ) then allocate( pmx( i1:i2, j1:j2, 1:lmr ) ) ; pmx = 0.0 call PMx_Integrate_3d( region, RF%sizepmx(k), pmx, status ) IF_NOTOK_RETURN(status=1) rf%data3d(:,:,:, rf%trec, k) = pmx ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( pmx(i1:i2,j1:j2,lms:lme), (/imr,jmr,lmr,1/) ), status & ! start=(/i1,j1,1,RF%trec/), count=(/imr,jmr,lmr,1/) ) deallocate( pmx ) else #endif ! --------- ! transported or chemistry only ? ! --------- select case( itr ) case( 1:ntracet ) ! ---------------------------------------------------- ! distinguish between mixing ratios and concentrations ! ---------------------------------------------------- select case( RF%varid_type(k) ) case( 'conc' ) ! write slab of concentrations ! m(trace) pressure xm(trace) ! C = -------- * fscale * ----------- * --------- ! m(air) temperature Rgas ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr) / & ! m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & ! pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & ! Rgas, (/imr,jmr,lmr,1/) ), & ! status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) ) #ifdef with_m7 rf%data3d(:,:,:, rf%trec, k) = mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr) / & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & Rgas #else write(gol,*)"Not using m7 - did not expected to be here."; call goErr write(gol,*)" - make pres3d available"; call goErr status=1; TRACEBACK; return #endif case( 'mixr' ) ! write slab of volume mixing ratios ! m(trace) ! X = -------- * fscale ! m(air) ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr)/ & ! m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr), & ! (/imr,jmr,lmr,1/) ), & ! status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) ) rf%data3d(:,:,:, rf%trec, k) = mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr)/ & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr) case( 'numb' ) ! write slab of concentrations ! number(trace) pressure #/gridcell Pa*K*mol ! C = ------------- * molmass_air * ---------------- = ------------- * kg/mol *----------- ! m(air) temperature*Rgas kg/gridcell K*J #ifdef with_m7 rf%data3d(:,:,:, rf%trec, k) = mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr) / & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & Rgas #else write(gol,*)"Not using m7 - did not expected to be here."; call goErr write(gol,*)" - make pres3d available"; call goErr status=1; TRACEBACK; return #endif case default write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr status=1 end select ! IF_NOTOK_MDF(fid=RF%ncid) ! --------- case( ntracet+1:ntrace ) ! --------- ! ---------------------------------------------------- ! distinguish between mixing ratios and concentrations ! ---------------------------------------------------- select case( RF%varid_type(k) ) case( 'conc' ) ! write slab of concentrations ! m(trace) pressure xm(trace) ! C = -------- * fscale * ----------- * --------- ! m(air) temperature Rgas ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr) / & ! m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & ! pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & ! Rgas, (/imr,jmr,lmr,1/) ), & ! status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) ) #ifdef with_m7 rf%data3d(:,:,:, rf%trec, k) = chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr) / & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & Rgas #else write(gol,*)"Not using m7 - did not expected to be here."; call goErr write(gol,*)" - make pres3d available"; call goErr status=1; TRACEBACK; return #endif case( 'mixr' ) ! write slab of volume mixing ratios ! m(trace) ! X = -------- * fscale ! m(air) ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr)/ & ! m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr), & ! (/imr,jmr,lmr,1/) ), & ! status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) ) rf%data3d(:,:,:, rf%trec, k) = chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr)/ & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr) case default write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr status=1 end select IF_NOTOK_MDF(fid=RF%ncid) ! --------- ! NOy and others (M7) ! --------- #ifdef with_m7 case( iNOy, iSO4, iBC, iPOM, iSS, iDU ) #else case( iNOy ) #endif listtrac(:) = -999 select case( itr ) case( iNOy ); numtrac = nNOyt; listtrac(1:nNOyt) = iNOyt #ifdef with_m7 case( iSO4 ); numtrac = nSO4t; listtrac(1:nSO4t) = iSO4t case( iBC ); numtrac = nBCt ; listtrac(1:nBCt ) = iBCt case( iPOM ); numtrac = nPOMt; listtrac(1:nPOMt) = iPOMt case( iSS ); numtrac = nSSt ; listtrac(1:nSSt ) = iSSt case( iDU ); numtrac = nDUt ; listtrac(1:nDUt ) = iDUt #endif end select ! mole fraction = sum of mole fractions of components ! storage for sum of components (distributed over levels): allocate( Compo_k(i1:i2,j1:j2,lmr) ) ! 3d fields with all levels or local levels only: allocate( field_k(i1:i2,j1:j2,lmr) ) ! loop over transported components: Compo_k = 0.0 do k_comp = 1, numtrac ! no more components?? if( listtrac(k_comp) < 0 ) exit ! global tracer index: itr_comp = listtrac(k_comp) ! check ... if ( itr_comp > ntracet ) then write (gol,'("index of NOy component does not represent a transported tracer : ",i3)') itr_comp; call goErr TRACEBACK; status=1; return end if ! ---------------------------------------------------- ! distinguish between mixing ratios and concentrations ! ---------------------------------------------------- select case( RF%varid_type(k) ) case( 'conc' ) ! calculate concentrations ! m(trace) pressure xm(trace) ! C = -------- * fscale * ----------- * --------- ! m(air) temperature Rgas #ifdef with_m7 field_k = mass_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr_comp) / & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * & xmair * 1.E-03 * pres3d(i1:i2,j1:j2,1:lmr) / & temper_dat(region)%data(i1:i2,j1:j2,1:lmr) / Rgas #else write(gol,*)"Not using m7 - did not expected to be here."; call goErr write(gol,*)" - make pres3d available"; call goErr status=1; TRACEBACK; return #endif case( 'mixr' ) ! m(trace) ! X = -------- * fscale ! m(air) field_k = mass_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr_comp) / & m_dat(region)%data(i1:i2,j1:j2,lms:lme) * & fscale(itr_comp) case default write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr TRACEBACK; status=1; return end select ! add contribution of this component: Compo_k = Compo_k + field_k end do ! write slab of volume mixing ratio's: ! call MDF_Put_Var( RF%ncid, RF%varid_tr(k), & ! reshape( Compo_k, (/imr,jmr,lmr,1/) ), & ! status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) ) ! IF_NOTOK_MDF(fid=RF%ncid) rf%data3d(:,:,:, rf%trec, k) = Compo_k ! clear deallocate( Compo_k ) deallocate( field_k ) ! ------------------- case default ! ------------------- write (gol,'("strange tracer index requested : ",i6)') itr; call goErr TRACEBACK; status=1; return end select #ifdef with_m7 endif #endif end do ! tracer !---------------- ! WRITE !---------------- if ( RF%trec == rf%n_rec ) then call MDF_Put_Var( RF%ncid, RF%varid_time, rf%time, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_date, rf%date, status) IF_NOTOK_MDF(fid=RF%ncid) ! surface presure call MDF_Put_Var( RF%ncid, RF%varid_ps, rf%sp, status, start=(/i1,j1,1/) ) IF_NOTOK_MDF(fid=RF%ncid) ! temperature (3d) call MDF_Put_Var( RF%ncid, RF%varid_temp, RF%data3d_t(:,:,:,:), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,RF%n_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) #ifdef tropomi if ( isRoot ) then write (gol,'(a,2i4)') 'PDUMP - writing fields T, hgt, ltropo, no2, so2, hcho; trec, n_rec ', RF%trec, rf%n_rec call goPr end if ! temperature (3d) call MDF_Put_Var( RF%ncid, RF%varid_temp, RF%data3d_t(:,:,:,:), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,RF%n_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) ! surface altitude (2d) call MDF_Put_Var( RF%ncid, RF%varid_hgt, RF%data2d_hgt(:,:), status, start=(/i1,j1/), count=(/imr,jmr/) ) IF_NOTOK_MDF(fid=RF%ncid) ! highest tropopause level (2d) call MDF_Put_Var( RF%ncid, RF%varid_ltropo, RF%data2d_ltropo(:,:,:), status, start=(/i1,j1,1/), count=(/imr,jmr,RF%n_rec/) ) IF_NOTOK_MDF(fid=RF%ncid) #endif ! vmr do k = 1, RF%ntr call MDF_Put_Var( RF%ncid, RF%varid_tr(k), RF%data3d(:,:,:,:,k), status, start=(/i1,j1,1,1/) ) IF_NOTOK_MDF(fid=RF%ncid) end do end if !---------------- ! DONE !---------------- #ifdef with_m7 deallocate(pres3d) #endif call goLabel() status = 0 END SUBROUTINE RF_VMR_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_VMR_Done ! ! !DESCRIPTION: close file #3 !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_VMR_Done( RF, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_VMR), intent(inout) :: RF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_VMR_Done' ! --- begin ------------------------------------- call goLabel(rname) call MDF_Close( RF%ncid, status ) IF_NOTOK_RETURN(status=1) deallocate(rf%date, rf%time, rf%sp, rf%data3d ) deallocate(rf%data3d_t) #ifdef tropomi deallocate(rf%data2d_hgt) deallocate(rf%data2d_ltropo) #endif call goLabel() ; status = 0 END SUBROUTINE RF_VMR_Done !EOC ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE: 2D LT output ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_LT_Init ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! subroutine RF_LT_Init( RF, fdir, model, expid, filetype, region, & idate_f, local_time, tracer_names, status ) ! ! !USES: ! use Binas, only : xmair use GO, only : goReadFromLine, goUpCase use GO, only : NewDate use dims, only : im, jm use chem_param, only : ntrace, names, ra use partools, only : MPI_INFO_NULL, localComm use MeteoData, only : global_lli, levi, sp_dat, Set ! ! !OUTPUT PARAMETERS: ! type(TPdumpFile_LT), intent(out) :: RF ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid character(len=*), intent(in) :: filetype integer, intent(in) :: region integer, intent(in) :: idate_f(6) integer, intent(in) :: local_time character(len=*), intent(in) :: tracer_names integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_LT_Init' ! --- local ------------------------------------ character(len=256) :: fname integer :: varid integer :: imr, jmr, lmr character(len=256) :: trnames character(len=8) :: trname, tmname character(len=3) :: cwavel integer :: k, itr, i1, i2, j1, j2 character(len=32) :: varname, varname_enti, varname_spec character(len=64) :: cf_medium_stnd, cf_medium_long character(len=64) :: cf_enti_stnd, cf_enti_long, cf_enti_unit character(len=64) :: cf_spec_stnd, cf_spec_long character(len=256) :: cf_name_stnd, cf_name_long, cf_name_unit character(len=512) :: comment ! --- begin ------------------------------------- call goLabel(rname) ! store arguments RF%local_time = local_time RF%tracer_names = tracer_names ! set tracer index for requested tracers: write (gol,'("selected tracers for LT output:")'); call goPr RF%ntr = 0 #ifdef with_m7 RF%laod = .false. RF%wavel = -1.0 #endif RF%itr = -1 trnames = tracer_names do ! empty ? if ( len_trim(trnames) == 0 ) exit ! next number: if ( RF%ntr == ntrace ) then write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr TRACEBACK; status=1; return end if RF%ntr = RF%ntr + 1 ! extract leading name: call goReadFromLine( trnames, trname, status, sep=' ' ) IF_NOTOK_RETURN(status=1) #ifdef with_m7 ! --------------------------- ! check for AOD ! --------------------------- if( strlowercase(trname(1:3)) == 'aod' ) then RF%laod(RF%ntr) = .true. RF%itr (RF%ntr) = -1 ! paste size to real read(trname(5:len_trim(trname)), * ) RF%wavel(RF%ntr) else #endif ! convert to tm5 name: select case ( trim(strlowercase(trname)) ) case ( 'hcho' ) ; tmname = 'CH2O' case ( 'rn', 'radon' ) ; tmname = 'Rn222' case ( 'pb', 'lead' ) ; tmname = 'Pb210' case default ; tmname = trname end select ! NOy is a special ... select case ( trim(strlowercase(tmname)) ) case ( 'NOy' ) ! defined as ntrace+1 RF%itr(RF%ntr) = iNOy write (gol,'(" * ",a10)') trim(trname); call goPr case default ! loop over all names: RF%itr(RF%ntr) = -1 do itr = 1, ntrace ! case indendent match ? if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then write (gol,'(" ",i3," ",a10," (",a10,") ",f12.4)') itr, trim(trname), trim(names(itr)), ra(itr); call goPr RF%itr(RF%ntr) = itr exit end if end do end select ! not found ? if ( RF%itr(RF%ntr) < 0 ) then write (gol,'("tracer name not supported:")'); call goPr write (gol,'(" list all : ",a)') trim(tracer_names); call goPr write (gol,'(" list element : ",i3)') RF%ntr; call goPr write (gol,'(" pdump name : ",a)') trim(trname); call goPr write (gol,'(" tm5 name : ",a)') trim(tmname); call goPr write (gol,'(" tm5 tracers : ")'); call goPr do itr = 1, ntrace write (gol,'(" ",i3," ",a)') itr, trim(names(itr)); call goPr end do TRACEBACK; status=1; return end if #ifdef with_m7 end if ! aod #endif ! store pdump name: RF%name_tr(RF%ntr) = trname end do ! empty file ? if ( RF%ntr < 1 ) then write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr TRACEBACK; status=1; return end if ! grid size imr = global_lli(region)%nlon jmr = global_lli(region)%nlat lmr = levi%nlev ! o open file ! write filename write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3) ! open: #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title' , 'local time output' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version' , trim(dataset_version) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'file_version_number', trim(outfileversnr) , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call MDF_Def_Dim( RF%ncid, 'lon' , global_lli(region)%nlon, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lat' , global_lli(region)%nlat, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lev' , levi%nlev , RF%dimid_lev , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'time' , 1 , RF%dimid_time , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'datelen', 6 , RF%dimid_datelen, status) IF_NOTOK_MDF(fid=RF%ncid) ! o define variables call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_east', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_north', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid call MDF_Def_Var( RF%ncid, 'lev', mdf_float, (/RF%dimid_lev/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'level' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , '1' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lev = varid call MDF_Def_Var( RF%ncid, 'time', mdf_float, (/RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'days since 1950-01-01 00:00:00', status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'calender' , 'gregorian' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_time = varid call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_date = varid call MDF_Def_Var( RF%ncid, 'ps', MDF_FLOAT, & (/RF%dimid_lon,RF%dimid_lat,RF%dimid_time/), varid, status ) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'surface_air_pressure', status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'surface pressure' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'Pa' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_ps = varid ! CF standard name for medium: cf_medium_stnd = 'in_air' ; cf_medium_long = 'in humid air' ! loop over tracer to be written: do k = 1, RF%ntr #ifdef with_m7 if( RF%laod(k) ) then ! get diameter write(cwavel,'(I3)') RF%wavel(k) ! Aerosol Optical Depth (AOD): varname_spec = 'AOD@'//trim(cwavel) cf_spec_stnd = 'AOD at '//trim(cwavel)//'nm' cf_spec_long = 'aerosol optical depth at '//trim(cwavel)//' nanometer' cf_enti_stnd = 'aerosol_optical_depth' cf_enti_unit = '1' cf_enti_long = 'aerosol optical depth' else #endif ! global tracer index itr = RF%itr(k) ! ~~ local time species info ! CF standard name for concentration/mixing ratio/column: cf_enti_stnd = 'mole_fraction' cf_enti_unit = 'mole mole-1' cf_enti_long = 'volume mixing ratio' ! start of dataset name: varname_enti = 'dry' ! no comment yet comment = '' ! standard names from CF conventions: select case ( RF%name_tr(k) ) case ( 'CO', 'co' ) varname_spec = 'co' cf_spec_stnd = 'carbon_monoxide' cf_spec_long = 'CO' case ( 'O3', 'o3' ) varname_spec = 'o3' cf_spec_stnd = 'ozone' cf_spec_long = 'O3' case ( 'O3s', 'o3s' ) varname_spec = 'o3s' cf_spec_stnd = 'ozone_from_stratosphere' cf_spec_long = 'O3s' case ( 'NO', 'no' ) varname_spec = 'no' cf_spec_stnd = 'nitrogen_monoxide' cf_spec_long = 'NO' case ( 'NO2', 'no2' ) varname_spec = 'no2' cf_spec_stnd = 'nitrogen_dioxide' cf_spec_long = 'NO2' case ( 'NOy', 'noy' ) varname_spec = 'noy' cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen' cf_spec_long = 'NOy' comment = 'NOy = NOx + HNO3 + PAN + org.ntr., '// & 'with NOx = NO + NO2 + NO3 + HNO4 + N2O5' case ( 'CH2O', 'ch2o', 'CHOH', 'choh' ) varname_spec = 'ch2o' cf_spec_stnd = 'formaldehyde' cf_spec_long = 'CH2O' case ( 'SO2', 'so2' ) varname_spec = 'so2' cf_spec_stnd = 'sulfur_dioxide' cf_spec_long = 'SO2' case ( 'CH4', 'ch4' ) varname_spec = 'ch4' cf_spec_stnd = 'methane' cf_spec_long = 'CH4' case ( 'OH', 'oh' ) varname_spec = 'oh' cf_spec_stnd = 'hydroxyl_radical' cf_spec_long = 'OH' case ( 'H2O2', 'h2o2' ) varname_spec = 'h2o2' cf_spec_stnd = 'hydrogen_peroxide' cf_spec_long = 'H2O2' case ( 'HNO3', 'hno3' ) varname_spec = 'hno3' cf_spec_stnd = 'nitric_acid' cf_spec_long = 'HNO3' case ( 'NH3', 'nh3' ) varname_spec = 'nh3' cf_spec_stnd = 'ammonia' cf_spec_long = 'NH3' case ( 'NH4', 'nh4' ) varname_spec = 'nh4' cf_spec_stnd = 'ammonium' cf_spec_long = 'NH4' case ( 'ORGNTR','orgntr' ) varname_spec = 'orgntr' cf_spec_stnd = 'organic_nitrate' cf_spec_long = 'ORGNTR' case ( 'PAN', 'pan' ) varname_spec = 'pan' cf_spec_stnd = 'peroxyacetyl_nitrate' cf_spec_long = 'PAN' case ( 'Rn', 'rn', 'Radon', 'radon' ) varname_spec = 'rn' cf_spec_stnd = 'radon' cf_spec_long = 'Rn' case ( 'Pb', 'pb', 'Lead', 'lead' ) varname_spec = 'pb' cf_spec_stnd = 'lead' cf_spec_long = 'Pb' case default write (gol,'("do not know how to match tracer with CF standard names : ",a)') RF%name_tr(k); call goPr TRACEBACK; status=1; return end select #ifdef with_m7 end if ! RF%laod(k) #endif ! define variable: call MDF_Def_Var( RF%ncid, trim(varname_spec), MDF_FLOAT, & (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) ! total names: cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)//'_'//trim(cf_medium_stnd) cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long)//' '//trim(cf_medium_long) cf_name_unit = trim(cf_enti_unit) ! write attributes: call MDF_Put_Att( RF%ncid, varid, 'standard_name', trim(cf_name_stnd) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', trim(cf_name_long) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units', trim(cf_name_unit) , status) IF_NOTOK_MDF(fid=RF%ncid) if ( itr <= ntrace .and. itr > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) else call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0 , status) IF_NOTOK_MDF(fid=RF%ncid) end if call MDF_Put_Att( RF%ncid, varid, 'moleweight_air', xmair*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'moleweight_unit', 'kg mole-1' , status) IF_NOTOK_MDF(fid=RF%ncid) if ( len_trim(comment) > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'comment', trim(comment) , status) IF_NOTOK_MDF(fid=RF%ncid) end if ! store varid RF%varid_tr(k) = varid end do ! o end defintion mode call MDF_EndDef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! no records written yet RF%trec = 0 call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate(RF%accu (i1:i2, j1:j2, 1:lmr, RF%ntr)) ; RF%accu = 0 allocate(RF%naccu (i1:i2, RF%ntr )) ; RF%naccu = 0 allocate(RF%p_accu (i1:i2, j1:j2 )) ; RF%p_accu = 0 allocate(RF%np_accu(i1:i2 )) ; RF%np_accu = 0 call goLabel() status = 0 END SUBROUTINE RF_LT_Init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_LT_Write ! ! !DESCRIPTION: does not write anything, just get !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_LT_Write( RF, region, idate_f, status ) ! ! !USES: ! use GO, only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol use chem_param, only : ntrace, ntracet, fscale use tracer_data, only : mass_dat, chem_dat use MeteoData, only : global_lli, levi, m_dat, sp_dat ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_LT), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_LT_Write' ! --- local ------------------------------------ integer :: imr, jmr, lmr, gimr, i1, i2, j1, j2 real, allocatable :: lev(:) real, allocatable :: field_out(:,:,:) real, allocatable :: field_out_b(:,:) integer :: l, ls, le type(TDate) :: t, t0 real :: time real :: dt_sec integer :: i, j, k, itr integer(kind=8) :: itau integer :: loctim, gridboxtimestep integer :: iloctim,itautoday,ilon integer :: icomp, itr_loc, ncells, window ! --- begin ------------------------------------- ! for multiple of dhour only ... ! if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0) ) then ! status=0; return ! end if call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 gimr = global_lli(region)%nlon ! gjmr = global_lli(region)%nlat lmr = levi%nlev ! next time record: RF%trec = RF%trec + 1 if(okdebug)then write(gol,*) "RF_LT_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr end if ! grid index offsets for GMT and local time loctim=RF%local_time if( loctim < 0 ) loctim=loctim+24*3600 ! time since 1950-1-1 00:00 t0 = NewDate( time6=time_reftime6 ) t = NewDate( time6=idate_f ) call SET( t, hour=0, min=0, sec=0 ) time = rTotal( t - t0, 'day' ) + loctim / 86400. ! ! ~~ time, grid ! ! only once ... if ( RF%trec == 1 ) then ! write longitudes: call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! write latitudes: call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! write level indices: allocate( lev(lmr) ) do l = 1, lmr lev(l) = real(l) end do call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status) IF_NOTOK_MDF(fid=RF%ncid) deallocate(lev) ! time: call MDF_Put_Var( RF%ncid, RF%varid_time, (/time/) , status, start=(/RF%trec/)) IF_NOTOK_MDF(fid=RF%ncid) ! date: call MDF_Put_Var( RF%ncid, RF%varid_date, reshape(real(idate_f),(/6,1/)), status, & start=(/1,1/), count=(/6,1/) ) IF_NOTOK_MDF(fid=RF%ncid) end if ! first record ! ! local time ! if ( RF%trec > 1 ) then ! do not accumulate fields on 00:00 ! grid index offsets for GMT and local time loctim=RF%local_time if( loctim < 0 ) loctim=loctim+24*3600 gridboxtimestep=24*3600/gimr itau = idate_f(4)*3600+idate_f(5)*60+idate_f(6) itautoday= nint(real(mod(itau,24*3600)*gimr)/real(24*3600)) iloctim = nint(real(loctim *gimr)/real(24*3600)) ! determine longitude index wrt Greenwich from difference (local time - GMT) ! also process neigboring longitudes (i-2 , i-1 , i , i+1 , i+2) depending on ! number of longitudinal grid cells ncells = ceiling( gimr / 24. ) window = ceiling( ncells / 2. ) do ilon = 1, ncells i = 1 + mod( gimr + gimr/2 + iloctim - itautoday + (ilon - window),gimr ) if (i .ge. i1 .and. i.le. i2) then RF%p_accu(i,j1:j2)= RF%p_accu(i,j1:j2)+sp_dat(region)%data(i,j1:j2,1) RF%np_accu(i)= RF%np_accu(i)+1 ! loop over tracers to be written: do k = 1, RF%ntr ! global tracer index: itr = RF%itr(k) !!$#ifdef with_m7 !!$ !!$ ! --------------------- !!$ ! AOD !!$ ! --------------------- !!$ if( RF%laod(k) ) then !!$ !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ allocate( ....... ) ) !!$ !!$ call PMx_Integrate_3d( region, RF%sizepmx(k), pmx, status ) !!$ IF_NOTOK_RETURN(status=1) !!$ !!$ ! root only: !!$ if ( myid == root ) then !!$ !!$ status = pnf90_put_var( RF%ncid, RF%varid_tr(k), & !!$ reshape( pmx(ims:ime,jms:jme,lms:lme), (/imr,jmr,lmr,1/) ), & !!$ start=(/1,1,1,RF%trec/), count=(/imr,jmr,lmr,1/) ) !!$ !!$ end if !!$ !!$ deallocate( ............. ) !!$ !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ ! .............. PUT HERE THE CODE FOR AOD ..................! !!$ else !!$ !!$#endif ! transported or chemistry only ? if ( (itr >= 1) .and. (itr <= ntracet) ) then RF%accu(i,j1:j2,1:lmr,k)= RF%accu(i,j1:j2,1:lmr,k)+& (mass_dat(region)%rm(i,j1:j2,1:lmr,itr)/ & m_dat(region)%data(i,j1:j2,1:lmr))*fscale(itr) RF%naccu(i,k)=RF%naccu(i,k)+1 else if ( (itr >= ntracet+1) .and. (itr <= ntrace) ) then RF%accu(i,j1:j2,1:lmr,k)= RF%accu(i,j1:j2,1:lmr,k)+& (chem_dat(region)%rm(i,j1:j2,1:lmr,itr)/ & m_dat(region)%data(i,j1:j2,1:lmr))*fscale(itr) RF%naccu(i,k)=RF%naccu(i,k)+1 end if enddo endif enddo endif ! do not accumulate fields on 00:00 call goLabel(); status = 0 END SUBROUTINE RF_LT_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_LT_Done ! ! !DESCRIPTION: write final data, then close file #4 !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_LT_Done( RF, region, status ) ! ! !USES: ! use MeteoData, only : global_lli, levi ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_LT), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! - move averaging & writing here ! !EOP !------------------------------------------------------------------------ !BOC character(len =*), parameter :: rname = mname//'/RF_LT_Done' integer :: imr, jmr real, allocatable :: field_out(:,:,:) real, allocatable :: field_out_b(:,:) integer :: i, ls, le, k, itr, i1, i2, j1, j2, lmr ! --- begin ------------------------------------- call goLabel(rname) !--------------------- ! average & write data !--------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev allocate(field_out_b(i1:i2,j1:j2)); field_out_b = 0.0 do i = i1, i2 if (RF%np_accu(i).gt.0) then field_out_b(i,:) =RF%p_accu(i,:)/RF%np_accu(i) endif enddo call MDF_Put_Var( RF%ncid, RF%varid_ps, reshape(field_out_b(i1:i2,j1:j2), & (/imr,jmr,1/) ), status, start=(/i1,j1,1/), count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=RF%ncid) deallocate(field_out_b) TRACERS: do k = 1, RF%ntr ! global tracer index: itr = RF%itr(k) if ( (itr >= 1) .and. (itr <= ntrace) ) then ! normalize fields, if necessary allocate(field_out(i1:i2,j1:j2,1:lmr)); field_out = 0.0 do i = i1,i2 if (RF%naccu(i,k).gt.0) then field_out(i,:,1:lmr) =RF%accu(i,:,1:lmr,k)/RF%naccu(i,k) endif enddo ! write fields: call MDF_Put_Var( RF%ncid, RF%varid_tr(k) , & reshape(field_out(i1:i2,j1:j2,1:lmr) , & (/imr,jmr,lmr,1/) ) , & status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,1/) ) IF_NOTOK_MDF(fid=RF%ncid) deallocate(field_out) endif end do TRACERS !--------------------- ! DONE !--------------------- call MDF_Close( RF%ncid , status) IF_NOTOK_RETURN(status=1) deallocate(RF%accu) deallocate(RF%naccu) deallocate(RF%p_accu) deallocate(RF%np_accu) call goLabel() ; status = 0 END SUBROUTINE RF_LT_Done !EOC #ifdef with_budgets ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE ##5 : 2D dry and wet deposition fields ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPS_Init ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! subroutine RF_DEPS_Init( RF, fdir, model, expid, filetype, region, & idate_f, dhour, tracer_names, status ) ! ! !USES: ! use Binas, only : xmair use GO, only : goReadFromLine, goUpCase use GO, only : NewDate use dims, only : im, jm use chem_param, only : ntrace, names, ra use partools, only : MPI_INFO_NULL, localComm use MeteoData, only : global_lli, levi ! ! !OUTPUT PARAMETERS: ! type(TPdumpFile_DEPS), intent(out) :: RF integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid character(len=*), intent(in) :: filetype integer, intent(in) :: region integer, intent(in) :: idate_f(6) integer, intent(in) :: dhour character(len=*), intent(in) :: tracer_names ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retor -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPS_Init' ! --- local ------------------------------------ character(len=256) :: fname integer :: varid character(len=256) :: trnames character(len=8) :: trname, tmname integer :: k, itr character(len=32) :: varname, varname_enti, varname_spec character(len=64) :: cf_enti_stnd, cf_enti_long, cf_enti_unit character(len=64) :: cf_spec_stnd, cf_spec_long character(len=256) :: cf_name_stnd, cf_name_long, cf_name_unit character(len=512) :: comment logical :: with_wdep integer :: imr, jmr, i1, i2, j1, j2 ! --- begin ------------------------------------- call goLabel(rname) ! -- store arguments, init var RF%dhour = dhour RF%tracer_names = tracer_names RF%ntr = 0 trnames = tracer_names ! -- get dims call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 ! Switch to default .false., requires an extra call to PDUMP_Files_Write2 in OUTPUT_PDUMP_DONE n_deps_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, mess='DPS_Init' ) !n_deps_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, .true., 'DPS_Init' ) if ( n_deps_rec == 0 ) then ! degenerated case deps_apply = .false. status=0 return end if ! -- tracer index for requested tracers: if ( len_trim(trnames) == 0 ) then deps_apply = .false. write (gol,'("WARNING - NO tracers selected for depositions output!")') ; call goPr write (gol,'(" - deps_apply set to False.")' ) ; call goPr status=0 return else write (gol,'("selected tracers for depositions output:")'); call goPr end if do if ( len_trim(trnames) == 0 ) exit ! next number: if ( RF%ntr == ntrace ) then write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr TRACEBACK; status=1; return end if RF%ntr = RF%ntr + 1 ! extract leading name: call goReadFromLine( trnames, trname, status, sep=' ' ) IF_NOTOK_RETURN(status=1) ! store pdump name: RF%name_tr(RF%ntr) = trname ! convert to tm5 name: select case ( trname ) case ( 'HCHO' ) ; tmname = 'CH2O' case ( 'Rn', 'Radon' ) ; tmname = 'Rn222' case ( 'Pb', 'Lead' ) ; tmname = 'Pb210' case default ; tmname = trname end select ! wet deposition ? with_wdep = .false. select case ( trname ) case ( 'HNO3' ) ; with_wdep = .true. case ( 'NOy' ) ; with_wdep = .true. case ( 'NH3' ) ; with_wdep = .true. case ( 'NH4' ) ; with_wdep = .true. case ( 'SO4' ) ; with_wdep = .true. end select RF%with_wdep(RF%ntr) = with_wdep ! NOy is a special ... select case ( tmname ) case ( 'NOy' ) ! defined as ntrace+1 RF%itr(RF%ntr) = iNOy write (gol,'(" ",i3," ",a10," (",a10,") ",f12.4,"; wdep : ",l1)') & -1,trim(trname), '*', -1.0, with_wdep; call goPr case default ! loop over all names: RF%itr(RF%ntr) = -1 do itr = 1, ntrace ! case indendent match ? if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then write (gol,'(" ",i3," ",a10," (",a10,") ",f12.4," ; wdep : ",l1)') & itr, trim(trname), trim(names(itr)), ra(itr), with_wdep; call goPr RF%itr(RF%ntr) = itr exit end if end do end select ! not found ? if ( RF%itr(RF%ntr) < 0 ) then write (gol,'("tracer name not supported:") ') ; call goPr write (gol,'(" list all : ",a) ') trim(tracer_names) ; call goPr write (gol,'(" list element : ",i3) ') RF%ntr ; call goPr write (gol,'(" pdump name : ",a) ') trim(trname) ; call goPr write (gol,'(" tm5 name : ",a) ') trim(tmname) ; call goPr write (gol,'(" tm5 tracers : ") ') ; call goPr do itr = 1, ntrace write (gol,'(" ",i3," ",a)') itr, trim(names(itr)); call goPr end do TRACEBACK; status=1; return end if end do ! empty file ? if ( RF%ntr < 1 ) then write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr TRACEBACK; status=1; return end if ! allocate storage: allocate( RF%ddep_budget(imr,jmr,RF%ntr) ) ; RF%ddep_budget = 0.0 allocate( RF%wdep_budget(imr,jmr,RF%ntr) ) ; RF%wdep_budget = 0.0 ! store current time (when budgets are reset): RF%t0_budget = NewDate(time6=idate_f) ! o open file write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3) #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title' , 'dry and wet deposition' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version' , trim(dataset_version) , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call MDF_Def_Dim( RF%ncid, 'lon' , global_lli(region)%nlon, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lat' , global_lli(region)%nlat, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'time' , n_deps_rec , RF%dimid_time , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'datelen', 6 , RF%dimid_datelen, status) IF_NOTOK_MDF(fid=RF%ncid) ! o define variables call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_east', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_north', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid call MDF_Def_Var( RF%ncid, 'time', mdf_float, (/RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'days since 1950-01-01 00:00:00', status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'calender' , 'gregorian' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_time = varid call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid , 'long_name', 'date and time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid , 'units' , 'year, month, day, hour, minute, second', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_date = varid call MDF_Def_Var( RF%ncid, 'accum', mdf_float , (/RF%dimid_time/) , varid, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid , 'long_name', 'length of accumulated time interval', status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid , 'units' , 'second' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_accum = varid allocate( RF%time(n_deps_rec) ) allocate( RF%date(6,n_deps_rec) ) allocate( RF%dt(n_deps_rec) ) ! loop over tracer to be written: do k = 1, RF%ntr ! global tracer index itr = RF%itr(k) ! ~~ dry deposition ! CF standard name for concentration/mixing ratio/column: cf_enti_stnd = 'surface_dry_deposition_mole_flux' cf_enti_unit = 'mole m-2 s-1' cf_enti_long = 'dry deposition of ' ! start of dataset name: varname_enti = 'dry' ! no comment yet comment = '' ! standard names from CF conventions: select case ( RF%name_tr(k) ) case ( 'CO', 'co' ) varname_spec = 'co' cf_spec_stnd = 'carbon_monoxide' cf_spec_long = 'CO' case ( 'O3', 'o3' ) varname_spec = 'o3' cf_spec_stnd = 'ozone' cf_spec_long = 'O3' case ( 'O3s', 'o3s' ) varname_spec = 'o3s' cf_spec_stnd = 'ozone_from_stratosphere' cf_spec_long = 'O3s' case ( 'NO', 'no' ) varname_spec = 'no' cf_spec_stnd = 'nitrogen_monoxide' cf_spec_long = 'NO' case ( 'NO2', 'no2' ) varname_spec = 'no2' cf_spec_stnd = 'nitrogen_dioxide' cf_spec_long = 'NO2' case ( 'NOy', 'noy' ) varname_spec = 'noy' cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen' cf_spec_long = 'NOy' comment = 'NOy = NOx + HNO3 + PAN + org.ntr., '// & 'with NOx = NO + NO2 + NO3 + HNO4 + N2O5' case ( 'CH2O', 'ch2o', 'CHOH', 'choh' ) varname_spec = 'ch2o' cf_spec_stnd = 'formaldehyde' cf_spec_long = 'CH2O' case ( 'SO2', 'so2' ) varname_spec = 'so2' cf_spec_stnd = 'sulfur_dioxide' cf_spec_long = 'SO2' case ( 'CH4', 'ch4' ) varname_spec = 'ch4' cf_spec_stnd = 'methane' cf_spec_long = 'CH4' case ( 'OH', 'oh' ) varname_spec = 'oh' cf_spec_stnd = 'hydroxyl_radical' cf_spec_long = 'OH' case ( 'H2O2', 'h2o2' ) varname_spec = 'h2o2' cf_spec_stnd = 'hydrogen_peroxide' cf_spec_long = 'H2O2' case ( 'HNO3', 'hno3' ) varname_spec = 'hno3' cf_spec_stnd = 'nitric_acid' cf_spec_long = 'HNO3' case ( 'NH3', 'nh3' ) varname_spec = 'nh3' cf_spec_stnd = 'ammonia' cf_spec_long = 'NH3' case ( 'ORGNTR','orgntr' ) varname_spec = 'orgntr' cf_spec_stnd = 'organic_nitrate' cf_spec_long = 'ORGNTR' case ( 'NH4', 'nh4' ) varname_spec = 'nh4' cf_spec_stnd = 'ammonium' cf_spec_long = 'NH4' case ( 'PAN', 'pan' ) varname_spec = 'pan' cf_spec_stnd = 'peroxyacetyl_nitrate' cf_spec_long = 'PAN' case ( 'Rn', 'rn', 'Radon', 'radon' ) varname_spec = 'rn' cf_spec_stnd = 'radon' cf_spec_long = 'Rn' case ( 'Pb', 'pb', 'Lead', 'lead' ) varname_spec = 'pb' cf_spec_stnd = 'lead' cf_spec_long = 'Pb' case default write (gol,'("do not know how to match tracer with CF standard names : ",a)') RF%name_tr(k); call goErr TRACEBACK; status=1; return end select write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec) ! define variable: call MDF_Def_Var( RF%ncid, trim(varname), MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) ! total names: cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd) cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long) cf_name_unit = trim(cf_enti_unit) ! write attributes: call MDF_Put_Att( RF%ncid , varid, 'standard_name' , trim(cf_name_stnd), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'long_name' , trim(cf_name_long), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'units' , trim(cf_name_unit), status) IF_NOTOK_MDF(fid=RF%ncid) if ( itr <= ntrace ) then call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) else call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0 , status) IF_NOTOK_MDF(fid=RF%ncid) end if call MDF_Put_Att( RF%ncid , varid, 'moleweight_air' , xmair*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_unit' , 'kg mole-1' , status) IF_NOTOK_MDF(fid=RF%ncid) if ( len_trim(comment) > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'comment' , trim(comment) , status) IF_NOTOK_MDF(fid=RF%ncid) end if ! store varid RF%varid_ddep(k) = varid ! ~~ wet deposition if ( RF%with_wdep(k) ) then ! CF standard name for concentration/mixing ratio/column: cf_enti_stnd = 'surface_wet_deposition_mole_flux' cf_enti_unit = 'mole m-2 s-1' cf_enti_long = 'wet deposition of ' ! start of dataset name: varname_enti = 'wet' ! by default no comment: comment = '' ! standard names from CF conventions: select case ( RF%name_tr(k) ) case ( 'NOy', 'noy' ) varname_spec = 'noy' cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen' cf_spec_long = 'NOy' comment = 'NOy = NOx + HNO3 + PAN + org.ntr., '// & 'with NOx = NO + NO2 + NO3 + HNO4 + N2O5' case ( 'HNO3', 'hno3' ) varname_spec = 'hno3' cf_spec_stnd = 'nitric_acid' cf_spec_long = 'HNO3' case ( 'NH3', 'nh3' ) varname_spec = 'nh3' cf_spec_stnd = 'ammonia' cf_spec_long = 'NH3' case ( 'NH4', 'nh4' ) varname_spec = 'nh4' cf_spec_stnd = 'ammonium' cf_spec_long = 'NH4' case ( 'SO2', 'so2' ) varname_spec = 'so2' cf_spec_stnd = 'sulfur_dioxide' cf_spec_long = 'SO2' case default write (gol,'("unsupported tracer name for CF standard name : ",a)') RF%name_tr(k); call goPr TRACEBACK; status=1; return end select write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec) ! define variable: call MDF_Def_Var( RF%ncid, trim(varname), MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) ! total names: cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd) cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long) cf_name_unit = trim(cf_enti_unit) ! write attributes: call MDF_Put_Att( RF%ncid , varid, 'standard_name' , trim(cf_name_stnd), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'long_name' , trim(cf_name_long), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'units' , trim(cf_name_unit), status) IF_NOTOK_MDF(fid=RF%ncid) if ( itr <= ntrace ) then call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) else call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0 , status) IF_NOTOK_MDF(fid=RF%ncid) end if call MDF_Put_Att( RF%ncid , varid, 'moleweight_air' , xmair*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_unit' , 'kg mole-1' , status) IF_NOTOK_MDF(fid=RF%ncid) if ( len_trim(comment) > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'comment' , trim(comment) , status) IF_NOTOK_MDF(fid=RF%ncid) end if ! store varid RF%varid_wdep(k) = varid end if end do allocate( RF%data2d_dry(i1:i2, j1:j2, n_deps_rec, RF%ntr) ) allocate( RF%data2d_wet(i1:i2, j1:j2, n_deps_rec, RF%ntr) ) ! RF%data2d_dry = 0. ! RF%data2d_wet = 0. ! o end defintion mode call MDF_EndDef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! o ! no records written yet RF%trec = 0 call goLabel() ! ok status = 0 end subroutine RF_DEPS_Init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPS_Write ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_DEPS_Write( RF, region, idate_f, status ) ! ! !USES: ! use GO, only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol use Grid, only : AreaOper use MeteoData, only : global_lli, levi, lli #ifndef without_chemistry use ebischeme, only : buddrydep_dat => buddep_dat #endif #ifndef without_wet_deposition use wet_deposition, only : buddep_dat #endif ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_DEPS), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPS_Write' ! --- local ------------------------------------ integer :: imr, jmr, lmr type(TDate) :: t, t0 real :: time real :: dt_sec integer :: k, itr, i1, i2, j1, j2 real, allocatable :: budget(:,:) real, allocatable :: budget_loc(:,:) real, allocatable :: depflux(:,:) integer :: icomp ! --- begin ------------------------------------- ! for multiple of dhour only ... if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0) ) then status=0; return end if call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr = i2-i1+1 jmr = j2-j1+1 lmr = levi%nlev ! temporary storage: allocate( budget_loc(imr,jmr) ) allocate( depflux (imr,jmr) ) ! next time record: RF%trec = RF%trec + 1 if(okdebug)then write(gol,*) "RF_DEPS_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr end if ! time since 1950-1-1 00:00 t0 = NewDate( time6=time_reftime6 ) t = NewDate( time6=idate_f ) time = rTotal( t - t0, 'day' ) ! length of time interval: dt_sec = iTotal( t - RF%t0_budget, 'sec' ) ! zero time interval ? routine should not have been called ... if ( dt_sec == 0 ) then write (gol,'("routine called after zero lenght time interval:")'); call goErr call wrtgol( ' t0_budget : ', RF%t0_budget ); call goErr call wrtgol( ' t : ', t ); call goErr !status=1 TRACEBACK end if ! reset timer: call Set( RF%t0_budget, time6=idate_f ) !--------------- ! Write GRID !--------------- if ( RF%trec == 1 ) then ! longitudes call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! latitudes call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status) IF_NOTOK_MDF(fid=RF%ncid) end if !--------------- ! FILL DIAGNOSTIC ARRAYS !--------------- !--------------- time rf%time(rf%trec) = time rf%date(:,rf%trec) = real(idate_f) rf%dt(rf%trec) = dt_sec !--------------- dry deposition do k = 1, RF%ntr ! global tracer index: itr = RF%itr(k) ! extract current budget #ifndef without_chemistry if ( itr == iNOy ) then ! add contributions of all NOy components: budget_loc = 0.0 do icomp = 1, nNOyt budget_loc = budget_loc + buddrydep_dat(region)%dry(:,:,iNOyt(icomp)) end do else ! extract budget for requested tracer: budget_loc = buddrydep_dat(region)%dry(:,:,itr) end if #else budget_loc = 0.0 #endif ! deposition flux ~ (current budget - previous budget)/dt depflux = ( budget_loc - RF%ddep_budget(:,:,k) ) / dt_sec ! mole/s call AreaOper( lli(region), depflux, '/', 'm2', status ) ! mole/m2/s IF_NOTOK_RETURN(status=1) ! save current budget & store record RF%ddep_budget(:,:,k) = budget_loc rf%data2d_dry(:,:,RF%trec,k)= depflux end do !--------------- wet deposition do k = 1, RF%ntr ! skip ? if ( .not. RF%with_wdep(k) ) cycle ! global tracer index: itr = RF%itr(k) ! extract current budget #ifndef without_wet_deposition if ( itr == iNOy ) then ! add contributions of all NOy components: budget_loc = 0.0 do icomp = 1, nNOyt ! add wet depositions for large scale and convective precip; total column: budget_loc = budget_loc + sum(buddep_dat(region)%lsp(:,:,:,iNOyt(icomp)),3) + & sum(buddep_dat(region)% cp(:,:,:,iNOyt(icomp)),3) end do else ! extract budget for requested tracer; ! add wet depositions for large scale and convective precip; total column: budget_loc = sum(buddep_dat(region)%lsp(:,:,:,itr),3) + & sum(buddep_dat(region)% cp(:,:,:,itr),3) end if #else budget_loc = 0.0 #endif ! deposition flux ~ (current budget - previous budget)/dt depflux = ( budget_loc - RF%wdep_budget(:,:,k) ) / dt_sec ! mole/s call AreaOper( lli(region), depflux, '/', 'm2', status ) ! mole/m2/s IF_NOTOK_RETURN(status=1) ! save current budget & store record RF%wdep_budget(:,:,k) = budget_loc RF%data2d_wet(:,:,RF%trec,k)= depflux end do !---------------- ! WRITE !---------------- if ( RF%trec == n_deps_rec ) then call MDF_Put_Var( RF%ncid, RF%varid_time, rf%time, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_date, rf%date, status) IF_NOTOK_MDF(fid=RF%ncid) ! accumulation interval call MDF_Put_Var( RF%ncid, RF%varid_accum, rf%dt , status) IF_NOTOK_MDF(fid=RF%ncid) ! deposition flux do k = 1, RF%ntr call MDF_Put_Var( RF%ncid, RF%varid_ddep(k), rf%data2d_dry(:,:,:,k), status, start=(/i1,j1,1/) ) IF_NOTOK_MDF(fid=RF%ncid) if ( .not. RF%with_wdep(k) ) cycle call MDF_Put_Var( RF%ncid, RF%varid_wdep(k), rf%data2d_wet(:,:,:,k), status, start=(/i1,j1,1/) ) IF_NOTOK_MDF(fid=RF%ncid) end do end if !---------------- ! DONE !---------------- deallocate( budget_loc ) deallocate( depflux ) call goLabel() status = 0 END SUBROUTINE RF_DEPS_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPS_Done ! ! !DESCRIPTION: close file #5 !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_DEPS_Done( RF, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_DEPS), intent(inout) :: RF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPS_Done' ! --- begin ------------------------------------- call goLabel(rname) ! close file call MDF_Close( RF%ncid , status) IF_NOTOK_RETURN(status=1) ! clear deallocate( RF%ddep_budget ) deallocate( RF%wdep_budget ) deallocate( rf%time, rf%date, rf%dt, rf%data2d_dry, rf%data2d_wet ) call goLabel() ; status = 0 END SUBROUTINE RF_DEPS_Done !EOC ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! FILE #6 : deposition velocities ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPV_Init ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! subroutine RF_DEPV_Init( RF, fdir, model, expid, filetype, region, & idate_f, dhour, tracer_names, status ) ! ! !USES: ! use Binas, only : xmair use GO, only : goReadFromLine, goUpCase use GO, only : NewDate use dims, only : im, jm use chem_param, only : ntrace, names, ra use partools, only : MPI_INFO_NULL, localComm use MeteoData, only : global_lli, levi ! ! !OUTPUT PARAMETERS: ! type(TPdumpFile_DEPV), intent(out) :: RF ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fdir character(len=*), intent(in) :: model character(len=*), intent(in) :: expid character(len=*), intent(in) :: filetype integer, intent(in) :: region integer, intent(in) :: idate_f(6) integer, intent(in) :: dhour character(len=*), intent(in) :: tracer_names integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPV_Init' ! --- local ------------------------------------ character(len=256) :: fname integer :: varid, i1, i2, j1, j2 character(len=256) :: trnames character(len=8) :: trname, tmname integer :: k, itr character(len=32) :: varname, varname_enti, varname_spec character(len=64) :: cf_enti_stnd, cf_enti_long, cf_enti_unit character(len=64) :: cf_spec_stnd, cf_spec_long character(len=256) :: cf_name_stnd, cf_name_long, cf_name_unit character(len=512) :: comment ! --- begin ------------------------------------- call goLabel(rname) ! store arguments RF%dhour = dhour RF%tracer_names = tracer_names RF%ntr = 0 trnames = tracer_names ! get dims call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) n_depv_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, mess='DEPV_Init' ) if ( n_depv_rec == 0 ) then ! degenerated cases depv_apply = .false. status=0 return end if ! tracer index for requested tracers if ( len_trim(trnames) == 0 ) then depv_apply = .false. write (gol,'("WARNING - NO tracers selected for depositions velocity output!")') ; call goPr write (gol,'(" - depv_apply set to False.")' ) ; call goPr status=0 return else write (gol,'("selected tracers for deposition velocity output:")'); call goPr end if do if ( len_trim(trnames) == 0 ) exit ! next number: if ( RF%ntr == ntrace ) then write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr TRACEBACK; status=1; return end if RF%ntr = RF%ntr + 1 ! extract leading name: call goReadFromLine( trnames, trname, status, sep=' ' ) IF_NOTOK_RETURN(status=1) ! store pdump name: RF%name_tr(RF%ntr) = trname ! convert to tm5 name: select case ( trname ) case ( 'HCHO' ) ; tmname = 'CH2O' case ( 'Rn', 'Radon' ) ; tmname = 'Rn222' case ( 'Pb', 'Lead' ) ; tmname = 'Pb210' case default ; tmname = trname end select ! loop over all names: RF%itr(RF%ntr) = -1 do itr = 1, ntrace ! case indendent match ? if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then write (gol,'(" ",i3," ",a10," (",a10,")",f12.4)') & itr, trim(trname), trim(names(itr)), ra(itr); call goPr RF%itr(RF%ntr) = itr exit end if end do ! not found ? if ( RF%itr(RF%ntr) < 0 ) then write (gol,'("tracer name not supported:") ') ; call goPr write (gol,'(" list all : ",a) ') trim(tracer_names) ; call goPr write (gol,'(" list element : ",i3) ') RF%ntr ; call goPr write (gol,'(" pdump name : ",a) ') trim(trname) ; call goPr write (gol,'(" tm5 name : ",a) ') trim(tmname) ; call goPr write (gol,'(" tm5 tracers : ") ') ; call goPr do itr = 1, ntrace write (gol,'(" ",i3," ",a)') itr, trim(names(itr)); call goPr end do TRACEBACK; status=1; return end if end do ! empty file ? if ( RF%ntr < 1 ) then write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr TRACEBACK; status=1; return end if ! o open file ! write filename write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') & trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3) ! open: #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status ) IF_NOTOK_RETURN(status=1) #endif ! o global attributes call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title' , 'volume mixing ratios' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author) , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution' , trim(institution), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version' , trim(dataset_version) , status) IF_NOTOK_MDF(fid=RF%ncid) ! o define dimensions call MDF_Def_Dim( RF%ncid, 'lon' , global_lli(region)%nlon, RF%dimid_lon , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'lat' , global_lli(region)%nlat, RF%dimid_lat , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'time' , n_depv_rec , RF%dimid_time , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Def_Dim( RF%ncid, 'datelen', 6 , RF%dimid_datelen, status) IF_NOTOK_MDF(fid=RF%ncid) ! o define variables call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'longitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_east', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lon = varid call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'latitude' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'degrees_north', status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_lat = varid call MDF_Def_Var( RF%ncid, 'time', mdf_float, (/RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name' , 'time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'days since 1950-01-01 00:00:00', status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'calender' , 'gregorian' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_time = varid allocate( rf%time(n_depv_rec) ) call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen, RF%dimid_time/), varid , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid, varid, 'units' , 'year, month, day, hour, minute, second' , status) IF_NOTOK_MDF(fid=RF%ncid) RF%varid_date = varid allocate( rf%date(6,n_depv_rec) ) ! loop over tracer to be written: do k = 1, RF%ntr ! global tracer index itr = RF%itr(k) ! CF standard name for concentration/mixing ratio/column: cf_enti_stnd = 'surface_dry_deposition_velocity_due_to_turbulence' cf_enti_unit = 'mole m-2 s-1' cf_enti_long = 'dry deposition of ' ! start of dataset name: varname_enti = 'ddepvel' ! by default no comment: comment = '' ! standard names from CF conventions: select case ( RF%name_tr(k) ) case ( 'CO', 'co' ) varname_spec = 'co' cf_spec_stnd = 'carbon_monoxide' cf_spec_long = 'CO' case ( 'O3', 'o3' ) varname_spec = 'o3' cf_spec_stnd = 'ozone' cf_spec_long = 'O3' case ( 'O3s', 'o3s' ) varname_spec = 'o3s' cf_spec_stnd = 'ozone_from_stratosphere' cf_spec_long = 'O3s' case ( 'NO', 'no' ) varname_spec = 'no' cf_spec_stnd = 'nitrogen_monoxide' cf_spec_long = 'NO' case ( 'NO2', 'no2' ) varname_spec = 'no2' cf_spec_stnd = 'nitrogen_dioxide' cf_spec_long = 'NO2' case ( 'NOy', 'noy' ) varname_spec = 'noy' cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen' cf_spec_long = 'NOy' comment = 'NOy = NOx + HNO3 + PAN + org.ntr., '// & 'with NOx = NO + NO2 + NO3 + HNO4 + N2O5' case ( 'CH2O', 'ch2o', 'CHOH', 'choh' ) varname_spec = 'ch2o' cf_spec_stnd = 'formaldehyde' cf_spec_long = 'CH2O' case ( 'SO2', 'so2' ) varname_spec = 'so2' cf_spec_stnd = 'sulfur_dioxide' cf_spec_long = 'SO2' case ( 'CH4', 'ch4' ) varname_spec = 'ch4' cf_spec_stnd = 'methane' cf_spec_long = 'CH4' case ( 'OH', 'oh' ) varname_spec = 'oh' cf_spec_stnd = 'hydroxyl_radical' cf_spec_long = 'OH' case ( 'H2O2', 'h2o2' ) varname_spec = 'h2o2' cf_spec_stnd = 'hydrogen_peroxide' cf_spec_long = 'H2O2' case ( 'HNO3', 'hno3' ) varname_spec = 'hno3' cf_spec_stnd = 'nitric_acid' cf_spec_long = 'HNO3' case ( 'PAN', 'pan' ) varname_spec = 'pan' cf_spec_stnd = 'peroxyacetyl_nitrate' cf_spec_long = 'PAN' case ( 'Rn', 'rn', 'Radon', 'radon' ) varname_spec = 'rn' cf_spec_stnd = 'radon' cf_spec_long = 'Rn' case ( 'Pb', 'pb', 'Lead', 'lead' ) varname_spec = 'pb' cf_spec_stnd = 'lead' cf_spec_long = 'Pb' case ( 'NH3', 'nh3' ) varname_spec = 'nh3' cf_spec_stnd = 'ammonia' cf_spec_long = 'NH3' case ( 'NH4', 'nh4' ) varname_spec = 'nh4' cf_spec_stnd = 'ammonium' cf_spec_long = 'NH4' case default write (gol,'("unsupported tracer name for CF standard name : ",a)') RF%name_tr(k); call goPr TRACEBACK; status=1; return end select write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec) write (gol,'(" varname : ",a)') trim(varname); call goPr ! define variable: call MDF_Def_Var( RF%ncid, trim(varname), MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status ) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status ) IF_NOTOK_MDF(fid=RF%ncid) ! total names: cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd) cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long) cf_name_unit = trim(cf_enti_unit) ! write attributes: call MDF_Put_Att( RF%ncid , varid, 'standard_name' , trim(cf_name_stnd), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'long_name' , trim(cf_name_long), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'units' , trim(cf_name_unit), status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_tracer', ra(itr)*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_air' , xmair*1e3 , status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Att( RF%ncid , varid, 'moleweight_unit' , 'kg mole-1' , status) IF_NOTOK_MDF(fid=RF%ncid) if ( len_trim(comment) > 0 ) then call MDF_Put_Att( RF%ncid, varid, 'comment' , trim(comment) , status) IF_NOTOK_MDF(fid=RF%ncid) end if ! store varid RF%varid_tr(k) = varid end do allocate( rf%data2d(i1:i2, j1:j2, n_depv_rec, rf%ntr) ) ! o end defintion mode call MDF_EndDef( RF%ncid , status) IF_NOTOK_MDF(fid=RF%ncid) ! o ! no records written yet RF%trec = 0 call goLabel() ; status = 0 END SUBROUTINE RF_DEPV_Init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPV_Write ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_DEPV_Write( RF, region, idate_f, status ) ! ! !USES: ! use GO, only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol use Grid, only : AreaOper use MeteoData, only : global_lli #ifndef without_dry_deposition use dry_deposition, only : vd #endif ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_DEPV), intent(inout) :: RF ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: idate_f(6) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPV_Write' ! --- local ------------------------------------ integer :: imr, jmr type(TDate) :: t, t0 real :: time integer :: k, itr, i1, i2, j1, j2 real, allocatable :: depvel(:,:) ! --- begin ------------------------------------- ! for multiple of dhour only ... if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0) ) then status=0; return end if call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 ! next time record: RF%trec = RF%trec + 1 if(okdebug)then write(gol,*) "RF_DEPV_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr end if ! time since 1950-1-1 00:00 t0 = NewDate( time6=time_reftime6 ) t = NewDate( time6=idate_f ) time = rTotal( t - t0, 'day' ) ! Only once : Dimensions if ( RF%trec == 1 ) then ! write longitudes: call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status) IF_NOTOK_MDF(fid=RF%ncid) ! write latitudes: call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status) IF_NOTOK_MDF(fid=RF%ncid) end if !-------- FILL DIAGNOSTIC ARRAYS rf%time(rf%trec) = time rf%date(:,rf%trec) = real(idate_f) ! loop over tracers to be written: do k = 1, RF%ntr itr = RF%itr(k) ! global tracer index #ifndef without_dry_deposition rf%data2d(:,:,RF%trec,k) = vd(region,itr)%surf ! deposition velocity #else rf%data2d(:,:,RF%trec,k) = 0.0 #endif end do !-------- WRITE if ( RF%trec == n_depv_rec ) then call MDF_Put_Var( RF%ncid, RF%varid_time, rf%time, status) IF_NOTOK_MDF(fid=RF%ncid) call MDF_Put_Var( RF%ncid, RF%varid_date, rf%date, status) IF_NOTOK_MDF(fid=RF%ncid) ! loop over tracers to be written: do k = 1, RF%ntr call MDF_Put_Var( RF%ncid, RF%varid_tr(k), rf%data2d(:,:,:,k), status, start=(/i1,j1,1/)) IF_NOTOK_MDF(fid=RF%ncid) end do end if call goLabel() status = 0 END SUBROUTINE RF_DEPV_Write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RF_DEPV_Done ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RF_DEPV_Done( RF, status ) ! ! !INPUT/OUTPUT PARAMETERS: ! type(TPdumpFile_DEPV), intent(inout) :: RF ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - retro -> pdump ! 7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/RF_DEPV_Done' ! --- begin ------------------------------------- call goLabel(rname) ! close file call MDF_Close( RF%ncid , status) IF_NOTOK_RETURN(status=1) deallocate( rf%time, rf%date, rf%data2d ) call goLabel() ; status = 0 END SUBROUTINE RF_DEPV_Done !EOC #endif !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: strlowercase ! ! !DESCRIPTION: ! ! This function returns a copy of the input string 'struppercase' with all ! letters changed to lowercase. All other characters remain unchanged. !\\ !\\ ! !INTERFACE: ! FUNCTION strlowercase(struppercase) ! ! !USES: ! IMPLICIT NONE ! ! !INPUT PARAMETERS: ! CHARACTER(LEN=*), INTENT(IN) :: struppercase ! ! !RETURN VALUE: ! CHARACTER(LEN=LEN(struppercase)) :: strlowercase ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! !EOP !------------------------------------------------------------------------ !BOC CHARACTER(LEN=1) :: u INTEGER :: i,j strlowercase = struppercase DO i=1,LEN(struppercase) u = struppercase(i:i) j = IACHAR(u) IF(j < 65 .OR. j > 90) CYCLE strlowercase(i:i) = ACHAR(j+32) END DO !------------------------------------------------------------------------------- END FUNCTION STRLOWERCASE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: struppercase ! ! !DESCRIPTION: ! ! This function returns a copy of the input string 'struppercase' with all ! letters changed to lowercase. All other characters remain unchanged. !\\ !\\ ! !INTERFACE: ! FUNCTION STRUPPERCASE(strlowercase) ! ! !USES: ! IMPLICIT NONE ! ! !INPUT PARAMETERS: ! CHARACTER(LEN=*), INTENT(IN) :: strlowercase ! ! !RETURN VALUE: ! CHARACTER(LEN=LEN(strlowercase)) :: struppercase ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - ! !EOP !------------------------------------------------------------------------ !BOC CHARACTER(LEN=1) :: u INTEGER :: i,j struppercase = strlowercase DO i=1,LEN(strlowercase) u = strlowercase(i:i) j = IACHAR(u) IF(j < 97 .OR. j > 122) CYCLE struppercase(i:i) = ACHAR(j-32) END DO !------------------------------------------------------------------------------- END FUNCTION STRUPPERCASE !EOC END MODULE USER_OUTPUT_PDUMP