Program Regrid_Util
use ESMF
use ESMFL_Mod
use MAPL_ExceptionHandling
use MAPL_Profiler
use MAPL_BaseMod
use MAPL_MemUtilsMod
use MAPL_CFIOMod
use MAPL_CommsMod
use MAPL_ShmemMod
use ESMF_CFIOMod
use ESMF_CFIOUtilMod
use ESMF_CFIOFileMod
use MAPL_NewRegridderManager
use MAPL_AbstractRegridderMod
use mapl_RegridMethods
use MAPL_GridManagerMod
use MAPL_LatLonGridFactoryMod, only: LatLonGridFactory
use MAPL_CubedSphereGridFactoryMod, only: CubedSphereGridFactory
use MAPL_TripolarGridFactoryMod, only: TripolarGridFactory
use MAPL_Constants, only: MAPL_PI_R8
use MAPL_ExceptionHandling
use MAPL_ApplicationSupport
use pFIO
use MAPL_ESMFFieldBundleWrite
use MAPL_ESMFFieldBundleRead
use MAPL_ServerManager
use MAPL_FileMetadataUtilsMod
use gFTL_StringVector
use regrid_util_support_mod
use mpi
implicit NONE
type(DistributedProfiler), target :: t_prof
type (ProfileReporter) :: reporter
call main()
CONTAINS
subroutine main()
type(regrid_support), target :: support
type(ESMF_VM) :: vm ! ESMF Virtual Machine
character(len=ESMF_MAXPATHLEN) :: Filename,OutputFile
integer :: myPET ! The local PET number
integer :: nPET ! The total number of PETs you are running on
integer :: status, rc
type(ESMF_FieldBundle) :: bundle
type(ESMF_Time) :: time
type(ESMF_Time), allocatable :: tSeries(:)
type(ESMF_TimeInterval) :: timeInterval
type(ESMF_Clock) :: clock
logical :: fileCreated,file_exists
integer :: tsteps,i,j,tint
type(VerticalData) :: vertical_data
type(FieldBundleWriter) :: newWriter
logical :: writer_created, has_vertical_level
type(ServerManager) :: io_server
call ESMF_Initialize (LogKindFlag=ESMF_LOGKIND_NONE, vm=vm, _RC)
call ESMF_VMGet(vm, localPET=myPET, petCount=nPet)
call MAPL_Initialize(_RC)
call MAPL_GetNodeInfo (comm=MPI_COMM_WORLD, _RC)
call ESMF_CalendarSetDefault ( ESMF_CALKIND_GREGORIAN, _RC )
call support%process_command_line(_RC)
t_prof=DistributedProfiler('Regrid_Util',MpiTimerGauge(),MPI_COMM_WORLD)
call t_prof%start(_RC)
call io_server%initialize(mpi_comm_world)
filename = support%filenames%at(1)
if (allocated(tSeries)) deallocate(tSeries)
call get_file_times(filename,support%itime,support%allTimes,tseries,timeInterval,tint,tsteps,_RC)
has_vertical_level = support%has_level(_RC)
if (has_vertical_level) then
call get_file_levels(filename,vertical_data,_RC)
end if
Clock = ESMF_ClockCreate ( name="Eric", timeStep=TimeInterval, &
startTime=tSeries(1), _RC )
bundle=ESMF_FieldBundleCreate(name="cfio_bundle",_RC)
call ESMF_FieldBundleSet(bundle,grid=support%new_grid,_RC)
writer_created=.false.
do j=1,support%filenames%size()
filename = support%filenames%at(j)
if (j>1) then
if (allocated(tSeries)) deallocate(tSeries)
call get_file_times(filename,support%itime,support%allTimes,tseries,timeInterval,tint,tsteps,_RC)
end if
outputfile = support%outputfiles%at(j)
inquire(file=trim(outputfile),exist=file_exists)
_ASSERT(.not.file_exists,"output file already exists: exiting!")
fileCreated=.false.
do i=1,tsteps
call t_prof%start("Read")
if (mapl_am_i_root()) write(*,*)'processing timestep from '//trim(filename)
time = tSeries(i)
if (support%onlyvars) then
call MAPL_Read_bundle(bundle,trim(filename),time=time,regrid_method=support%regridMethod,only_vars=support%vars,file_weights=support%use_weights, _RC)
else
call MAPL_Read_bundle(bundle,trim(filename),time=time,regrid_method=support%regridMethod,file_weights=support%use_weights, _RC)
end if
call t_prof%stop("Read")
call MPI_BARRIER(MPI_COMM_WORLD,STATUS)
_VERIFY(status)
call t_prof%start("write")
if (mapl_am_I_root()) write(*,*) "moving on to writing "//trim(outputfile)
call ESMF_ClockSet(clock,currtime=time,_RC)
if (.not. writer_created) then
call newWriter%create_from_bundle(bundle,clock,n_steps=tsteps,time_interval=tint,nbits_to_keep=support%shave,deflate=support%deflate,vertical_data=vertical_data,quantize_algorithm=support%quantize_algorithm,quantize_level=support%quantize_level,_RC)
writer_created=.true.
end if
if (.not.fileCreated) then
call newWriter%start_new_file(outputFile,_RC)
fileCreated=.true.
end if
call newWriter%write_to_file(_RC)
call t_prof%stop("write")
end do
enddo
! All done
! --------
call ESMF_VMBarrier(VM,_RC)
call io_server%finalize()
call t_prof%stop()
call t_prof%reduce()
call t_prof%finalize()
call generate_report()
call MAPL_Finalize(_RC)
call ESMF_Finalize ( _RC )
end subroutine main
subroutine get_file_levels(filename,vertical_data,rc)
character(len=*), intent(in) :: filename
type(VerticalData), intent(inout) :: vertical_data
integer, intent(out), optional :: rc
integer :: status
type(NetCDF4_fileFormatter) :: formatter
type(FileMetadata) :: basic_metadata
type(FileMetadataUtils) :: metadata
character(len=:), allocatable :: lev_name
character(len=ESMF_MAXSTR) :: long_name
character(len=ESMF_MAXSTR) :: standard_name
character(len=ESMF_MAXSTR) :: vcoord
character(len=ESMF_MAXSTR) :: lev_units
real, allocatable, target :: levs(:)
real, pointer :: plevs(:)
call formatter%open(trim(filename),pFIO_Read,_RC)
basic_metadata=formatter%read(_RC)
call metadata%create(basic_metadata,trim(filename))
lev_name = metadata%get_level_name(_RC)
call metadata%get_coordinate_info(lev_name,coords=levs,coordUnits=lev_units,long_name=long_name,&
standard_name=standard_name,coordinate_attr=vcoord,_RC)
plevs => levs
vertical_data = VerticalData(levels=plevs,vunit=lev_units,vcoord=vcoord,standard_name=standard_name,long_name=long_name, &
force_no_regrid=.true.,_RC)
nullify(plevs)
_RETURN(_SUCCESS)
end subroutine get_file_levels
subroutine get_file_times(filename,itime,alltimes,tseries,timeInterval,tint,tsteps,rc)
character(len=*), intent(in) :: filename
integer, intent(in) :: itime(2)
logical, intent(in) :: alltimes
type(ESMF_Time), allocatable, intent(inout) :: tseries(:)
type(ESMF_TimeInterval), intent(inout) :: timeInterval
integer, intent(out) :: tint
integer, intent(out) :: tsteps
integer, intent(out), optional :: rc
integer :: status
integer :: second,minute,hour,day,month,year
type(NetCDF4_fileFormatter) :: formatter
type(FileMetadata) :: basic_metadata
type(FileMetadataUtils) :: metadata
call formatter%open(trim(filename),pFIO_Read,_RC)
basic_metadata=formatter%read(_RC)
call metadata%create(basic_metadata,trim(filename))
call formatter%close(_RC)
tsteps = metadata%get_dimension('time',_RC)
call metadata%get_time_info(timeVector=tSeries,_RC)
if (.not.allTimes) then
tSteps=1
call UnpackDateTIme(itime,year,month,day,hour,minute,second)
deallocate(tSeries)
allocate(tSeries(1))
call ESMF_TimeSet(tSeries(1), yy=year, mm=month, dd=day, h=hour, m=minute, s=second,_RC)
end if
if (tSteps == 1) then
call ESMF_TimeIntervalSet( TimeInterval, h=6, m=0, s=0, _RC )
else
TimeInterval=tSeries(2)-tSeries(1)
end if
call ESMF_TimeIntervalGet(TimeInterval,h=hour,m=minute,s=second,_RC)
tint=hour*10000+minute*100+second
_RETURN(_SUCCESS)
end subroutine get_file_times
subroutine UnpackDateTime(DATETIME, YY, MM, DD, H, M, S)
integer, intent(IN ) :: DATETIME(:)
integer, intent( OUT) :: YY, MM, DD, H, M, S
YY = datetime(1)/10000
MM = mod(datetime(1),10000)/100
DD = mod(datetime(1),100)
H = datetime(2)/10000
M = mod(datetime(2),10000)/100
S = mod(datetime(2),100)
return
end subroutine UnpackDateTime
subroutine generate_report()
character(:), allocatable :: report_lines(:)
integer :: i
character(1) :: empty(0)
reporter = ProfileReporter(empty)
call reporter%add_column(NameColumn(20))
call reporter%add_column(FormattedTextColumn('Inclusive','(f9.6)', 9, InclusiveColumn('MEAN')))
call reporter%add_column(FormattedTextColumn('% Incl','(f6.2)', 6, PercentageColumn(InclusiveColumn('MEAN'),'MAX')))
call reporter%add_column(FormattedTextColumn('Exclusive','(f9.6)', 9, ExclusiveColumn('MEAN')))
call reporter%add_column(FormattedTextColumn('% Excl','(f6.2)', 6, PercentageColumn(ExclusiveColumn('MEAN'))))
call reporter%add_column(FormattedTextColumn(' Max Excl)','(f9.6)', 9, ExclusiveColumn('MAX')))
call reporter%add_column(FormattedTextColumn(' Min Excl)','(f9.6)', 9, ExclusiveColumn('MIN')))
call reporter%add_column(FormattedTextColumn('Max PE)','(1x,i5.5,1x)', 7, ExclusiveColumn('MAX_PE')))
call reporter%add_column(FormattedTextColumn('Min PE)','(1x,i5.5,1x)', 7, ExclusiveColumn('MIN_PE')))
report_lines = reporter%generate_report(t_prof)
if (mapl_am_I_root()) then
write(*,'(a)')'Final profile'
write(*,'(a)')'============='
do i = 1, size(report_lines)
write(*,'(a)') report_lines(i)
end do
write(*,'(a)') ''
end if
end subroutine generate_report
end program Regrid_Util