#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 ! ! 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 ! 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 : vmr1 ! output.pdump.vmr.001.dhour : 3 ! output.pdump.vmr.001.tracers : CO2 ! ! output.pdump.lt.apply : T ! output.pdump.lt.tracers : CO2 ! output.pdump.lt.localtime : 2 ! ! output.pdump.lt2.apply : F ! output.pdump.lt2.tracers : CO2 ! output.pdump.lt2.localtime : 12 ! ! output.pdump.depositions.apply : F ! output.pdump.depositions.dhour : 3 ! output.pdump.depositions.tracers : CO2 ! ! output.pdump.depvels.apply : F ! output.pdump.depvels.dhour : 3 ! output.pdump.depvels.tracers : CO2 ! !\\ !\\ ! !INTERFACE: ! MODULE USER_OUTPUT_PDUMP ! ! !USES: ! use GO, only : gol, goPr, goErr, goLabel use GO, only : TDate use dims, only : nregions, idatee, idatei, okdebug use chem_param, only : ntrace 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' ! ! 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 integer :: dhour 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_temp integer :: varid_ps integer :: varid_a_bnds, varid_b_bnds integer :: ntr integer :: itr(ntrace) character(len=8) :: name_tr(ntrace) integer :: varid_tr(ntrace) character(len=4) :: varid_type(ntrace) real, allocatable :: data3d(:,:,:,:,:) real, allocatable :: data3d_t(:,:,:,:) real, allocatable :: sp(:,:,:) real, allocatable :: time(:) real, allocatable :: date(:,:) 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(:) 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 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(:) integer, 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 ! 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 ! ! !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. ! !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, dh, 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 (date[4] should be 1 unless 1st day) integer, intent(in) :: dh ! time step for timeseries (should divide 24) 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" ! ! !REMARKS: ! - assumed dynamic timestep of 1h. It may be smaller at some steps (CFL ! violation) but that's ok, since it is never zero and we are looking ! at "integer-hours" for output. Only issue is if dynamic timestep is ! LARGER than timestep of timeseries: deemed too exotic, but we could ! add a test in init to avoid that. ! ! !TODO: ! - need to check if Oct 2013 fix works with runs that do NOT stop at ! midnight (labelled "lastday" cases in the code). Deemed exotic for now. ! - 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) if (deps) delta=date(4)+1 ! one DYNAMIC time step done to output something if (modulo(delta,dh)==0) then is=delta/dh else is=(delta+dh)/dh end if ! End index ie = 23 / dh ! 23 = 24 - dynamic time step 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/dh end if ! Case of last day stopping before midnite. (Need testing again - see ! !TODO. Probaly alright for VMR/TP/DEPV, but check required for DEPS) if (lastday) ie=idatee(4)/dh ! 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 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, dhour_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) :: dhour_min 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. 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) ! 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 end if ! time resolution (1 hour by default) call ReadRc( rcF, 'output.pdump.tp.dhour', tp_dhour, status, default=1 ) IF_ERROR_RETURN(status=1) ! 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. ) 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) ! 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 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) 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 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) end if ! no files open yet curr_day = -1 ! lowest time frequency is 1 hour dhour_min = 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 ) ! ! !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 ! --- begin ------------------------------- ! close files: do region = nregions, 1, -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 ! 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 ! 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 ! ! deposition velocities: ! if ( depv_apply ) then ! call RF_DEPV_Write( RF_DEPV(region), region, idate_f, status ) ! IF_NOTOK_RETURN(status=1) ! end if 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 ------------------------------- ! deposition fluxes: ! if ( deps_apply ) then ! call RF_DEPS_Write( RF_DEPS(region), region, idate_f, status ) ! IF_NOTOK_RETURN(status=1) ! end if ! 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 ------------------------------- ! 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 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, 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 CO2 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 : 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) integer, 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 ! !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 ! --- begin ------------------------------------- call goLabel(rname) ! store arguments RF%dhour = dhour RF%tracer_names = tracer_names ! 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, dhour, 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 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) ! convert to tm5 name: select case ( trim(strlowercase(trname)) ) case default ; tmname = trname end select ! -------------------------------- ! NOy and M7 are special cases ... ! -------------------------------- select case ( trim(strlowercase(tmname)) ) 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 ! 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 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' , '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)=='od') then model_meteo='forecast (IFS)' 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, 'time', rf%n_rec, RF%dimid_time , 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, '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_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 ) 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 allocate( RF%sp(i1:i2, j1:j2, rf%n_rec) ) #ifdef tropomi 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 allocate( RF%data3d_t(i1:i2, j1:j2, levi%nlev, rf%n_rec) ) #endif 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 ! loop over tracer to be written: do k = 1, RF%ntr ! ---------------------------- ! setting defaults (gas phase) ! ---------------------------- ! 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' 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 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 ! 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, temper_dat, sp_dat ! ! !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 real, allocatable :: lev(:) integer :: l type(TDate) :: t, t0 real :: time integer :: k, itr integer :: k_comp, itr_comp integer :: ims, ime, jms, jme, lms, lme integer :: gimr, gjmr, glmr real, allocatable :: compo_k(:,:,:) real, allocatable :: field_t(:,:,:) real, allocatable :: field_k(:,:,:) real, allocatable :: pres3d(:,:,:), pmx(:,:,:) integer :: numtrac integer :: listtrac(10) ! --- 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 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 ! next time record: RF%trec = RF%trec + 1 if(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) ! 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) 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 RF%data3d_t(:,:,:,RF%trec) = temper_dat(region)%data(i1:i2,j1:j2,1:lmr) #endif ! loop over all tracers to be written: do k = 1, RF%ntr ! global tracer index: itr = RF%itr(k) ! --------- ! 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/) ) 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 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 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/) ) 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 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) ! ------------------- case default ! ------------------- write (gol,'("strange tracer index requested : ",i6)') itr; call goErr TRACEBACK; status=1; return end select 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) #ifdef tropomi ! temperature 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) #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 !---------------- 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 ) #ifdef tropomi deallocate(rf%data3d_t) #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 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) ! convert to tm5 name: select case ( trim(strlowercase(trname)) ) case default ; tmname = trname end select ! NOy is a special ... select case ( trim(strlowercase(tmname)) ) 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 ! 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 ! 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 ( 'CO2', 'co2' ) varname_spec = 'co2' cf_spec_stnd = 'carbon_dioxide' cf_spec_long = 'CO2' 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 ! 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, itau, 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) ! 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 !-------------------------------------------------------------------------- ! 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