module tracer_data 27,8
!-----------------------------------------------------------------------
! module used to read (and interpolate) offline tracer data (sources and
! mixing ratios)
! Created by: Francis Vitt -- 2 May 2006
! Modified by : Jim Edwards -- 10 March 2009
!-----------------------------------------------------------------------
use shr_kind_mod
, only : r8 => shr_kind_r8,r4 => shr_kind_r4, shr_kind_cl, SHR_KIND_CS
use time_manager, only : get_curr_date, get_step_size, get_calendar, get_curr_calday
use spmd_utils, only : masterproc
use ppgrid, only : pcols, pver, pverp, begchunk, endchunk
use abortutils, only : endrun
use cam_logfile, only : iulog
use phys_buffer, only : pbuf, pbuf_get_fld_idx
use time_manager, only : set_time_float_from_date, set_date_from_time_float
use pio, only : file_desc_t, var_desc_t, &
pio_seterrorhandling, pio_internal_error, pio_bcast_error, &
pio_setdebuglevel, &
pio_char, pio_noerr, &
pio_inq_dimid, pio_inq_varid, &
pio_def_dim, pio_def_var, &
pio_put_att, pio_put_var, &
pio_get_var, pio_get_att, pio_nowrite, pio_inq_dimlen, &
pio_inq_vardimid, pio_inq_dimlen, pio_closefile
implicit none
private ! all unless made public
save
public :: trfld, input3d, input2d, trfile
public :: trcdata_init
public :: advance_trcdata
public :: get_fld_data
public :: get_fld_ndx
public :: write_trc_restart
public :: read_trc_restart
public :: init_trc_restart
! !PUBLIC MEMBERS
type input3d
real(r8), dimension(:,:,:), pointer :: data
endtype input3d
type input2d
real(r8), dimension(:,:), pointer :: data
endtype input2d
type trfld
real(r8), dimension(:,:,:), pointer :: data
type(input3d), dimension(4) :: input
character(len=32) :: srcnam
character(len=32) :: fldnam
character(len=32) :: units
type(var_desc_t) :: var_id
integer :: chnk_offset
integer :: coords(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM
integer :: order(4) ! LATDIM | LONDIM | LEVDIM | TIMDIM
endtype trfld
type trfile
type(input2d), dimension(4) :: ps_in
character(len=shr_kind_cl) :: pathname = ' '
character(len=shr_kind_cl) :: curr_filename = ' '
character(len=shr_kind_cl) :: next_filename = ' '
type(file_desc_t) :: curr_fileid
type(file_desc_t) :: next_fileid
type(var_desc_t), pointer :: currfnameid ! pio restart file var id
type(var_desc_t), pointer :: nextfnameid ! pio restart file var id
character(len=shr_kind_cl) :: filenames_list = ''
real(r8) :: datatimem = -1.e36_r8 ! time of prv. values read in
real(r8) :: datatimep = -1.e36_r8 ! time of nxt. values read in
real(r8) :: datatimes(4)
integer :: interp_recs
real(r8), pointer, dimension(:) :: curr_data_times
real(r8), pointer, dimension(:) :: next_data_times
logical :: remove_trc_file = .false. ! delete file when finished with it
real(r8) :: offset_time
integer :: cyc_ndx_beg
integer :: cyc_ndx_end
integer :: cyc_yr = 0
real(r8) :: one_yr = 0
real(r8) :: curr_mod_time ! model time - calendar day
real(r8) :: next_mod_time ! model time - calendar day - next time step
integer :: nlon
integer :: nlat
integer :: nlev
integer :: nilev
integer :: ps_coords(3) ! LATDIM | LONDIM | TIMDIM
integer :: ps_order(3) ! LATDIM | LONDIM | TIMDIM
real(r8), pointer, dimension(:) :: lons
real(r8), pointer, dimension(:) :: lats
real(r8), pointer, dimension(:) :: levs
real(r8), pointer, dimension(:) :: ilevs
real(r8), pointer, dimension(:) :: hyam
real(r8), pointer, dimension(:) :: hybm
real(r8), pointer, dimension(:,:) :: ps
real(r8) :: p0
type(var_desc_t) :: ps_id
logical :: has_ps = .false.
logical :: in_pbuf = .false.
logical :: zonal_ave = .false.
logical :: surf_data = .false.
logical :: alt_data = .false.
logical :: cyclical = .false.
logical :: fill_in_months = .false.
logical :: fixed = .false.
logical :: initialized = .false.
logical :: top_bndry = .false.
endtype trfile
integer, public, parameter :: MAXTRCRS = 100
integer, parameter :: LONDIM = 1
integer, parameter :: LATDIM = 2
integer, parameter :: LEVDIM = 3
integer, parameter :: TIMDIM = 4
integer, parameter :: PS_TIMDIM = 3
integer, parameter :: ZA_LATDIM = 1
integer, parameter :: ZA_LEVDIM = 2
integer, parameter :: ZA_TIMDIM = 3
integer, parameter :: nm=1 ! array index for previous (minus) data
integer, parameter :: np=2 ! array indes for next (plus) data
character(len=SHR_KIND_CS) :: calendar
contains
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine trcdata_init( specifier, filename, filelist, datapath, flds, file, & 5,36
rmv_file, data_ymd, data_tod, data_type )
use mo_constants, only : d2r
use cam_control_mod, only : nsrest
use commap, only : clat, clon
use dyn_grid, only : get_dyn_grid_parm
use string_utils, only : to_upper
implicit none
character(len=*), intent(in) :: specifier(:)
character(len=*), intent(in) :: filename
character(len=*), intent(in) :: filelist
character(len=*), intent(in) :: datapath
type(trfld), dimension(:), pointer :: flds
type(trfile), intent(inout) :: file
logical, intent(in) :: rmv_file
integer, intent(in) :: data_ymd
integer, intent(in) :: data_tod
character(len=*), intent(in) :: data_type
integer :: f, mxnflds, astat
integer :: str_yr, str_mon, str_day
integer :: lon_dimid, lat_dimid, lev_dimid, tim_dimid
integer :: dimids(4), did
type(var_desc_t) :: varid
integer :: idx
integer :: plon, plat, ierr
real(r8) :: start_time, time1, time2
calendar = get_calendar()
call specify_fields( specifier, flds )
file%datatimep=-1.e36_r8
file%datatimem=-1.e36_r8
mxnflds = 0
if (associated(flds)) mxnflds = size( flds )
if (mxnflds < 1) return
file%remove_trc_file = rmv_file
file%pathname = trim(datapath)
file%filenames_list = trim(filelist)
file%fill_in_months = .false.
file%cyclical = .false.
! does not work when compiled with pathf90
! select case ( to_upper(data_type) )
select case ( data_type )
case( 'FIXED' )
file%fixed = .true.
case( 'INTERP_MISSING_MONTHS' )
file%fill_in_months = .true.
case( 'CYCLICAL' )
file%cyclical = .true.
file%cyc_yr = data_ymd/10000
case( 'SERIAL' )
case default
write(iulog,*) 'trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename)
write(iulog,*) 'trcdata_init: valid data types: SERIAL | CYCLICAL | FIXED | INTERP_MISSING_MONTHS '
call endrun('trcdata_init: invalid data type: '//trim(data_type)//' file: '//trim(filename))
endselect
if (masterproc) then
write(iulog,*) 'trcdata_init: data type: '//trim(data_type)//' file: '//trim(filename)
endif
if ( len_trim(file%curr_filename)<1 .or. file%fixed ) then ! initial run
file%curr_filename = trim(filename)
call get_model_time(file)
if ( data_ymd > 0 .and. .not. file%cyclical ) then
str_yr = data_ymd/10000
str_mon = (data_ymd - str_yr*10000)/100
str_day = data_ymd - str_yr*10000 - str_mon*100
call set_time_float_from_date( start_time, str_yr, str_mon, str_day, data_tod )
file%offset_time = start_time - file%curr_mod_time
else
file%offset_time = 0
endif
endif
call set_time_float_from_date( time2, 2, 1, 1, 0 )
call set_time_float_from_date( time1, 1, 1, 1, 0 )
file%one_yr = time2-time1
if ( file%cyclical ) then
file%cyc_ndx_beg = -1
file%cyc_ndx_end = -1
if ( file%cyc_yr /= 0 ) then
call set_time_float_from_date( time1, file%cyc_yr , 1, 1, 0 )
call set_time_float_from_date( time2, file%cyc_yr+1, 1, 1, 0 )
file%one_yr = time2-time1
endif
call open_trc_datafile( file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times, &
cyc_ndx_beg=file%cyc_ndx_beg, cyc_ndx_end=file%cyc_ndx_end, cyc_yr=file%cyc_yr )
else
call open_trc_datafile( file%curr_filename, file%pathname, file%curr_fileid, file%curr_data_times )
file%curr_data_times = file%curr_data_times - file%offset_time
endif
call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR)
ierr = pio_inq_dimid( file%curr_fileid, 'lon', idx )
call pio_seterrorhandling(File%curr_fileid, PIO_INTERNAL_ERROR)
file%zonal_ave = (ierr/=PIO_NOERR)
plon = get_dyn_grid_parm('plon')
plat = get_dyn_grid_parm('plat')
if ( .not. file%zonal_ave ) then
call get_dimension( file%curr_fileid, 'lon', file%nlon, dimid=lon_dimid, data=file%lons )
file%lons = file%lons * d2r
endif
ierr = pio_inq_dimid( file%curr_fileid, 'time', tim_dimid )
call get_dimension( file%curr_fileid, 'lat', file%nlat, dimid=lat_dimid, data=file%lats )
file%lats = file%lats * d2r
allocate( file%ps(file%nlon,file%nlat), stat=astat )
if( astat /= 0 ) then
write(iulog,*) 'trcdata_init: file%ps allocation error = ',astat
call endrun('trcdata_init: failed to allocate x array')
end if
call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR)
ierr = pio_inq_varid( file%curr_fileid, 'PS', file%ps_id )
file%has_ps = (ierr==PIO_NOERR)
ierr = pio_inq_dimid( file%curr_fileid, 'altitude', idx )
file%alt_data = (ierr==PIO_NOERR)
if(file%alt_data) then
file%surf_data=.false.
else
ierr = pio_inq_dimid( file%curr_fileid, 'lev', idx )
file%surf_data = (ierr/=PIO_NOERR)
end if
call pio_seterrorhandling(File%curr_fileid, PIO_INTERNAL_ERROR)
if ( file%has_ps) then
ierr = pio_inq_vardimid (file%curr_fileid, file%ps_id, dimids(1:3))
do did = 1,3
if ( dimids(did) == lon_dimid ) then
file%ps_coords(LONDIM) = did
file%ps_order(did) = LONDIM
else if ( dimids(did) == lat_dimid ) then
file%ps_coords(LATDIM) = did
file%ps_order(did) = LATDIM
else if ( dimids(did) == tim_dimid ) then
file%ps_coords(PS_TIMDIM) = did
file%ps_order(did) = PS_TIMDIM
endif
enddo
endif
if (masterproc) then
write(iulog,*) 'trcdata_init: file%has_ps = ' , file%has_ps
endif ! masterproc
if (file%alt_data) then
call get_dimension( file%curr_fileid, 'altitude_int', file%nilev, data=file%ilevs )
call get_dimension( file%curr_fileid, 'altitude', file%nlev, dimid=lev_dimid, data=file%levs )
else
if (.not.file%surf_data) then
call get_dimension( file%curr_fileid, 'lev', file%nlev, dimid=lev_dimid, data=file%levs )
file%levs = file%levs*100._r8 ! mbar->pascals
else
file%nlev = 1
lev_dimid = -1
endif
endif
if (file%has_ps) then
allocate( file%hyam(file%nlev), file%hybm(file%nlev), stat=astat )
if( astat /= 0 ) then
write(iulog,*) 'trcdata_init: file%hyam,file%hybm allocation error = ',astat
call endrun('trcdata_init: failed to allocate file%hyam and file%hybm arrays')
end if
call pio_seterrorhandling(File%curr_fileid, PIO_BCAST_ERROR)
ierr = pio_inq_varid( file%curr_fileid, 'P0', varid)
call pio_seterrorhandling(File%curr_fileid, PIO_INTERNAL_ERROR)
if ( ierr == PIO_NOERR ) then
ierr = pio_get_var( file%curr_fileid, varid, file%p0 )
else
file%p0 = 100000._r8
endif
ierr = pio_inq_varid( file%curr_fileid, 'hyam', varid )
ierr = pio_get_var( file%curr_fileid, varid, file%hyam )
ierr = pio_inq_varid( file%curr_fileid, 'hybm', varid )
ierr = pio_get_var( file%curr_fileid, varid, file%hybm )
allocate( file %ps (pcols,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate file%ps array; error = ',astat
call endrun
end if
allocate( file%ps_in(1)%data(pcols,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(1)%data array; error = ',astat
call endrun
end if
allocate( file%ps_in(2)%data(pcols,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate file%ps_in(2)%data array; error = ',astat
call endrun
end if
if( file%fill_in_months ) then
allocate( file%ps_in(3)%data(pcols,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(*,*) 'trcdata_init: failed to allocate file%ps_in(3)%data array; error = ',astat
call endrun
end if
allocate( file%ps_in(4)%data(pcols,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(*,*) 'trcdata_init: failed to allocate file%ps_in(4)%data array; error = ',astat
call endrun
end if
end if
endif
flds_loop: do f = 1,mxnflds
! allocate memory only if not already in pbuf
if ( .not. file%in_pbuf ) then
if ( file%surf_data .or. file%top_bndry ) then
allocate( flds(f) %data(pcols,1,begchunk:endchunk), stat=astat )
else
allocate( flds(f) %data(pcols,pver,begchunk:endchunk), stat=astat )
endif
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate flds(f)%data array; error = ',astat
call endrun
end if
flds(f)%chnk_offset = 0
else
idx = pbuf_get_fld_idx(flds(f)%fldnam, failcode=-1 )
flds(f)%data => pbuf(idx)%fld_ptr(1,:,:,:,1)
flds(f)%chnk_offset = -begchunk+1
endif
allocate( flds(f)%input(1)%data(pcols,file%nlev,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(1)%data array; error = ',astat
call endrun
end if
allocate( flds(f)%input(2)%data(pcols,file%nlev,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'trcdata_init: failed to allocate flds(f)%input(2)%data array; error = ',astat
call endrun
end if
if( file%fill_in_months ) then
allocate( flds(f)%input(3)%data(pcols,file%nlev,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(*,*) 'trcdata_init: failed to allocate flds(f)%input(3)%data array; error = ',astat
call endrun
end if
allocate( flds(f)%input(4)%data(pcols,file%nlev,begchunk:endchunk), stat=astat )
if( astat/= 0 ) then
write(*,*) 'trcdata_init: failed to allocate flds(f)%input(4)%data array; error = ',astat
call endrun
end if
endif
ierr = pio_inq_varid( file%curr_fileid, flds(f)%srcnam, flds(f)%var_id )
if ( file%zonal_ave ) then
ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids(1:3))
do did = 1,3
if ( dimids(did) == lat_dimid ) then
flds(f)%coords(ZA_LATDIM) = did
flds(f)%order(did) = ZA_LATDIM
else if ( dimids(did) == lev_dimid ) then
flds(f)%coords(ZA_LEVDIM) = did
flds(f)%order(did) = ZA_LEVDIM
else if ( dimids(did) == tim_dimid ) then
flds(f)%coords(ZA_TIMDIM) = did
flds(f)%order(did) = ZA_TIMDIM
endif
enddo
else if ( file%surf_data ) then
ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids(1:3))
do did = 1,3
if ( dimids(did) == lon_dimid ) then
flds(f)%coords(LONDIM) = did
flds(f)%order(did) = LONDIM
else if ( dimids(did) == lat_dimid ) then
flds(f)%coords(LATDIM) = did
flds(f)%order(did) = LATDIM
else if ( dimids(did) == tim_dimid ) then
flds(f)%coords(PS_TIMDIM) = did
flds(f)%order(did) = PS_TIMDIM
endif
enddo
else
ierr = pio_inq_vardimid (file%curr_fileid, flds(f)%var_id, dimids)
do did = 1,4
if ( dimids(did) == lon_dimid ) then
flds(f)%coords(LONDIM) = did
flds(f)%order(did) = LONDIM
else if ( dimids(did) == lat_dimid ) then
flds(f)%coords(LATDIM) = did
flds(f)%order(did) = LATDIM
else if ( dimids(did) == lev_dimid ) then
flds(f)%coords(LEVDIM) = did
flds(f)%order(did) = LEVDIM
else if ( dimids(did) == tim_dimid ) then
flds(f)%coords(TIMDIM) = did
flds(f)%order(did) = TIMDIM
endif
enddo
endif
flds(f)%units = ''
ierr = pio_get_att( file%curr_fileid, flds(f)%var_id, 'units', flds(f)%units)
enddo flds_loop
endsubroutine trcdata_init
!-----------------------------------------------------------------------
! Reads more data if needed and interpolates data to current model time
!-----------------------------------------------------------------------
subroutine advance_trcdata( flds, file, state ) 5,4
use physics_types,only : physics_state
implicit none
type(trfile), intent(inout) :: file
type(trfld), intent(inout) :: flds(:)
type(physics_state), intent(in):: state(begchunk:endchunk)
real(r8) :: data_time
if ( .not.( file%fixed .and. file%initialized ) ) then
call get_model_time(file)
data_time = file%datatimep
if ( file%cyclical ) then
! wrap around
if ( (file%datatimep<file%datatimem) .and. (file%curr_mod_time>file%datatimem) ) then
data_time = data_time + file%one_yr
endif
endif
if ( file%curr_mod_time > data_time ) then
call read_next_trcdata( flds, file )
end if
endif
! need to inperpolate the data, regardless !
! each mpi task needs to interpolate
call interpolate_trcdata( flds, file, state )
file%initialized = .true.
end subroutine advance_trcdata
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine get_fld_data( flds, field_name, data, ncol, lchnk )
implicit none
type(trfld), intent(inout) :: flds(:)
character(len=*), intent(in) :: field_name
real(r8), intent(out) :: data(:,:)
integer, intent(in) :: lchnk
integer, intent(in) :: ncol
integer :: f, nflds
data(:,:) = 0._r8
nflds = size(flds)
do f = 1, nflds
if ( trim(flds(f)%fldnam) == trim(field_name) ) then
data(:ncol,:) = flds(f)%data(:ncol,:,lchnk+flds(f)%chnk_offset )
endif
enddo
end subroutine get_fld_data
!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine get_fld_ndx( flds, field_name, idx )
implicit none
type(trfld), intent(in) :: flds(:)
character(len=*), intent(in) :: field_name
integer, intent(out) :: idx
integer :: f, nflds
idx = -1
nflds = size(flds)
do f = 1, nflds
if ( trim(flds(f)%fldnam) == trim(field_name) ) then
idx = f
return
endif
enddo
end subroutine get_fld_ndx
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine get_model_time(file) 5,9
implicit none
type(trfile), intent(inout) :: file
integer yr, mon, day, ncsec ! components of a date
call get_curr_date(yr, mon, day, ncsec)
if ( file%cyclical ) yr = file%cyc_yr
call set_time_float_from_date( file%curr_mod_time, yr, mon, day, ncsec )
file%next_mod_time = file%curr_mod_time + get_step_size()/86400._r8
end subroutine get_model_time
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine check_files( file ) 2,15
use shr_sys_mod, only: shr_sys_system
use ioFileMod, only: getfil
implicit none
type(trfile), intent(inout) :: file
!-----------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------
character(len=shr_kind_cl) :: ctmp
character(len=shr_kind_cl) :: loc_fname
integer :: istat, astat
if (file%next_mod_time > file%curr_data_times(size(file%curr_data_times))) then
if ( .not. associated(file%next_data_times) ) then
! open next file if not already opened...
file%next_filename = incr_filename( file%curr_filename, filenames_list=file%filenames_list, datapath=file%pathname )
call open_trc_datafile( file%next_filename, file%pathname, file%next_fileid, file%next_data_times )
file%next_data_times = file%next_data_times - file%offset_time
endif
endif
if ( associated(file%next_data_times) ) then
if (file%curr_mod_time >= file%next_data_times(1)) then
! close current file ...
call pio_closefile( file%curr_fileid )
! remove if requested
if( file%remove_trc_file ) then
call getfil( file%curr_filename, loc_fname, 0 )
write(iulog,*) 'check_files: removing file = ',trim(loc_fname)
ctmp = 'rm -f ' // trim(loc_fname)
write(iulog,*) 'check_files: fsystem issuing command - '
write(iulog,*) trim(ctmp)
call shr_sys_system( ctmp, istat )
end if
file%curr_filename = file%next_filename
file%curr_fileid = file%next_fileid
deallocate( file%curr_data_times, stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'check_files: failed to deallocate file%curr_data_times array; error = ',astat
call endrun
end if
allocate( file%curr_data_times( size( file%next_data_times ) ), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'check_files: failed to allocate file%curr_data_times array; error = ',astat
call endrun
end if
file%curr_data_times(:) = file%next_data_times(:)
file%next_filename = ''
deallocate( file%next_data_times, stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'check_files: failed to deallocate file%next_data_times array; error = ',astat
call endrun
end if
nullify( file%next_data_times )
endif
endif
end subroutine check_files
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
function incr_filename( filename, filenames_list, datapath ) 2,20
!-----------------------------------------------------------------------
! ... Increment or decrement a date string withing a filename
! the filename date section is assumed to be of the form
! yyyy-dd-mm
!-----------------------------------------------------------------------
use string_utils, only : incstr
use shr_file_mod, only : shr_file_getunit, shr_file_freeunit
implicit none
character(len=*), intent(in) :: filename ! present dynamical dataset filename
character(len=*), optional, intent(in) :: filenames_list
character(len=*), optional, intent(in) :: datapath
character(len=shr_kind_cl) :: incr_filename ! next filename in the sequence
! set new next_filename ...
!-----------------------------------------------------------------------
! ... local variables
!-----------------------------------------------------------------------
integer :: pos, pos1, istat
character(len=shr_kind_cl) :: fn_new, line, filepath
character(len=6) :: seconds
character(len=5) :: num
integer :: ios,unitnumber
if (( .not. present(filenames_list)) .or.(len_trim(filenames_list) == 0)) then
!-----------------------------------------------------------------------
! ... ccm type filename
!-----------------------------------------------------------------------
pos = len_trim( filename )
fn_new = filename(:pos)
if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(fn_new)
if( fn_new(pos-2:) == '.nc' ) then
pos = pos - 3
end if
istat = incstr( fn_new(:pos), 1 )
if( istat /= 0 ) then
write(iulog,*) 'incr_flnm: incstr returned ', istat
write(iulog,*) ' while trying to decrement ',trim( fn_new )
call endrun
end if
else
! open filenames_list
if ( masterproc ) write(iulog,*) 'incr_flnm: old filename = ',trim(filename)
if ( masterproc ) write(iulog,*) 'incr_flnm: open filenames_list : ',trim(filenames_list)
unitnumber = shr_file_getUnit()
if ( present(datapath) ) then
filepath = trim(datapath) //'/'// trim(filenames_list)
else
filepath = trim(datapath)
endif
open( unit=unitnumber, file=filepath, iostat=ios, status="OLD")
if (ios /= 0) then
call endrun('not able to open filenames_list file: '//trim(filepath))
endif
! read file names
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
endif
do while( trim(line) /= trim(filename) )
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
endif
enddo
read( unit=unitnumber, fmt='(A)', iostat=ios ) line
if (ios /= 0) then
call endrun('not able to increment file name from filenames_list file: '//trim(filenames_list))
endif
fn_new = trim(line)
close(unit=unitnumber)
call shr_file_freeUnit(unitnumber)
endif
incr_filename = trim(fn_new)
if ( masterproc ) write(iulog,*) 'incr_flnm: new filename = ',trim(incr_filename)
end function incr_filename
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
subroutine find_times( itms, fids, time, file, datatimem, datatimep, times_found ) 5,5
implicit none
type(trfile), intent(in) :: file
real(r8), intent(out) :: datatimem, datatimep
integer, intent(out) :: itms(2) ! record numbers that bracket time
type(file_desc_t), intent(out) :: fids(2) ! ids of files that contains these recs
real(r8), intent(in) :: time ! time of interest
logical, intent(inout) :: times_found
integer :: np1 ! current forward time index of dataset
integer :: n,i !
integer :: curr_tsize, next_tsize, all_tsize
integer :: astat
integer :: cyc_tsize
real(r8), allocatable, dimension(:):: all_data_times
curr_tsize = size(file%curr_data_times)
next_tsize = 0
if ( associated(file%next_data_times)) next_tsize = size(file%next_data_times)
all_tsize = curr_tsize + next_tsize
allocate( all_data_times( all_tsize ), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'find_times: failed to allocate all_data_times array; error = ',astat
call endrun
end if
all_data_times(:curr_tsize) = file%curr_data_times(:)
if (next_tsize > 0) all_data_times(curr_tsize+1:all_tsize) = file%next_data_times(:)
if ( .not. file%cyclical ) then
if ( all( all_data_times(:) > time ) ) then
write(iulog,*) 'FIND_TIMES: ALL data times are after ', time
write(iulog,*) 'FIND_TIMES: data times: ',all_data_times(:)
write(iulog,*) 'FIND_TIMES: time: ',time
call endrun('find_times: all(all_data_times(:) > time) '// trim(file%curr_filename) )
endif
! find bracketing times
find_times_loop : do n=1, all_tsize-1
np1 = n + 1
datatimem = all_data_times(n) !+ file%offset_time
datatimep = all_data_times(np1) !+ file%offset_time
if ( (time .ge. datatimem) .and. (time .le. datatimep) ) then
times_found = .true.
exit find_times_loop
endif
enddo find_times_loop
else
cyc_tsize = file%cyc_ndx_end - file%cyc_ndx_beg + 1
if ( cyc_tsize > 1 ) then
call findplb(all_data_times(file%cyc_ndx_beg:file%cyc_ndx_end),cyc_tsize, time, n )
if (n == cyc_tsize) then
np1 = 1
else
np1 = n+1
endif
datatimem = all_data_times(n +file%cyc_ndx_beg-1)
datatimep = all_data_times(np1+file%cyc_ndx_beg-1)
times_found = .true.
endif
endif
if ( .not. times_found ) then
if (masterproc) then
write(iulog,*)'FIND_TIMES: Failed to find dates bracketing desired time =', time
write(iulog,*)' datatimem = ',file%datatimem
write(iulog,*)' datatimep = ',file%datatimep
write(iulog,*)' all_data_times = ',all_data_times
!call endrun()
return
endif
endif
deallocate( all_data_times, stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'find_times: failed to deallocate all_data_times array; error = ',astat
call endrun
end if
if ( .not. file%cyclical ) then
itms(1) = n
itms(2) = np1
else
itms(1) = n +file%cyc_ndx_beg-1
itms(2) = np1 +file%cyc_ndx_beg-1
endif
fids(:) = file%curr_fileid
do i=1,2
if ( itms(i) > curr_tsize ) then
itms(i) = itms(i) - curr_tsize
fids(i) = file%next_fileid
endif
enddo
end subroutine find_times
!------------------------------------------------------------------------
!------------------------------------------------------------------------
subroutine read_next_trcdata( flds, file ) 1,23
implicit none
type (trfile), intent(inout) :: file
type (trfld),intent(inout) :: flds(:)
integer :: recnos(4),i,f,nflds !
integer :: cnt4(4) ! array of counts for each dimension
integer :: strt4(4) ! array of starting indices
integer :: cnt3(3) ! array of counts for each dimension
integer :: strt3(3) ! array of starting indices
type(file_desc_t) :: fids(4)
logical :: times_found
integer :: cur_yr, cur_mon, cur_day, cur_sec, yr1, yr2, mon, date, sec
real(r8) :: series1_time, series2_time
type(file_desc_t) :: fid1, fid2
nflds = size(flds)
times_found = .false.
do while( .not. times_found )
call find_times( recnos, fids, file%curr_mod_time, file,file%datatimem, file%datatimep, times_found )
if ( .not. times_found ) then
call check_files( file )
endif
enddo
file%interp_recs = 2
if ( file%fill_in_months ) then
if( file%datatimep-file%datatimem > file%one_yr ) then
call get_curr_date(cur_yr, cur_mon, cur_day, cur_sec)
call set_date_from_time_float(file%datatimem, yr1, mon, date, sec )
call set_date_from_time_float(file%datatimep, yr2, mon, date, sec )
call set_time_float_from_date( series1_time, yr1, cur_mon, cur_day, cur_sec )
call set_time_float_from_date( series2_time, yr2, cur_mon, cur_day, cur_sec )
fid1 = fids(1)
fid2 = fids(2)
file%cyclical = .true.
call set_cycle_indices( fid1, file%cyc_ndx_beg, file%cyc_ndx_end, yr1)
call find_times( recnos(1:2), fids(1:2), series1_time, file, file%datatimes(1), file%datatimes(2), times_found )
if ( .not. times_found ) then
call endrun('read_next_trcdata: time not found for series1_time')
endif
call set_cycle_indices( fid2, file%cyc_ndx_beg, file%cyc_ndx_end, yr2)
if ( fid1%fh /= fid2%fh ) then
file%cyc_ndx_beg = file%cyc_ndx_beg + size(file%curr_data_times)
file%cyc_ndx_end = file%cyc_ndx_end + size(file%curr_data_times)
endif
call find_times( recnos(3:4), fids(3:4), series2_time, file, file%datatimes(3), file%datatimes(4), times_found )
if ( .not. times_found ) then
call endrun('read_next_trcdata: time not found for series2_time')
endif
file%cyclical = .false.
file%interp_recs = 4
call set_date_from_time_float( file%datatimes(1), yr1, mon, date, sec )
call set_time_float_from_date( file%datatimem, cur_yr, mon, date, sec )
if (file%datatimes(1) > file%datatimes(2) ) then ! wrap around
if ( cur_mon == 1 ) then
call set_time_float_from_date( file%datatimem, cur_yr-1, mon, date, sec )
endif
endif
call set_date_from_time_float( file%datatimes(2), yr1, mon, date, sec )
call set_time_float_from_date( file%datatimep, cur_yr, mon, date, sec )
if (file%datatimes(1) > file%datatimes(2) ) then ! wrap around
if ( cur_mon == 12 ) then
call set_time_float_from_date( file%datatimep, cur_yr+1, mon, date, sec )
endif
endif
endif
endif
!
! Set up hyperslab corners
!
strt4(:) = 1
strt3(:) = 1
do i=1,file%interp_recs
do f = 1,nflds
if ( file%zonal_ave ) then
cnt3(flds(f)%coords(ZA_LATDIM)) = file%nlat
cnt3(flds(f)%coords(ZA_LEVDIM)) = file%nlev
cnt3(flds(f)%coords(ZA_TIMDIM)) = 1
strt3(flds(f)%coords(ZA_TIMDIM)) = recnos(i)
call read_za_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt3, cnt3, file, &
(/ flds(f)%order(ZA_LATDIM),flds(f)%order(ZA_LEVDIM) /) )
else if ( file%surf_data ) then
cnt3( flds(f)%coords(LONDIM)) = file%nlon
cnt3( flds(f)%coords(LATDIM)) = file%nlat
cnt3( flds(f)%coords(PS_TIMDIM)) = 1
strt3(flds(f)%coords(PS_TIMDIM)) = recnos(i)
call read_2d_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data(:,1,:), strt3, cnt3, file, &
(/ flds(f)%order(LONDIM),flds(f)%order(LATDIM) /) )
else
cnt4(flds(f)%coords(LONDIM)) = file%nlon
cnt4(flds(f)%coords(LATDIM)) = file%nlat
cnt4(flds(f)%coords(LEVDIM)) = file%nlev
cnt4(flds(f)%coords(TIMDIM)) = 1
strt4(flds(f)%coords(TIMDIM)) = recnos(i)
call read_3d_trc( fids(i), flds(f)%var_id, flds(f)%input(i)%data, strt4, cnt4, file, &
(/ flds(f)%order(LONDIM),flds(f)%order(LATDIM),flds(f)%order(LEVDIM) /))
endif
enddo
if ( file%has_ps ) then
cnt3(file%ps_coords(LONDIM)) = file%nlon
cnt3(file%ps_coords(LATDIM)) = file%nlat
cnt3(file%ps_coords(PS_TIMDIM)) = 1
strt3(file%ps_coords(PS_TIMDIM)) = recnos(i)
call read_2d_trc( fids(i), file%ps_id, file%ps_in(i)%data, strt3, cnt3, file, &
(/ file%ps_order(LONDIM),file%ps_order(LATDIM) /) )
endif
enddo
end subroutine read_next_trcdata
!------------------------------------------------------------------------
subroutine read_2d_trc( fid, vid, loc_arr, strt, cnt, file, order ) 2,17
use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
use phys_grid, only : pcols, begchunk, endchunk, get_ncols_p, get_rlat_all_p, get_rlon_all_p
use mo_constants, only : pi
use dycore, only: dycore_is
use polar_avg, only: polar_average
implicit none
type(file_desc_t), intent(in) :: fid
type(var_desc_t), intent(in) :: vid
integer, intent(in) :: strt(:), cnt(:), order(2)
real(r8),intent(out) :: loc_arr(:,:)
type (trfile), intent(in) :: file
real(r8) :: to_lats(pcols), to_lons(pcols), wrk(pcols)
real(r8), allocatable, target :: wrk2d(:,:)
real(r8), pointer :: wrk2d_in(:,:)
integer :: tsize, c, i, j, ierr, ncols
real(r8), parameter :: zero=0_r8, twopi=2_r8*pi
type(interp_type) :: lon_wgts, lat_wgts
nullify(wrk2d_in)
allocate( wrk2d(cnt(1),cnt(2)), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr
call endrun
end if
if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlon .or. cnt(2)/=file%nlat) then
allocate( wrk2d_in(file%nlon, file%nlat), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr
call endrun
end if
end if
ierr = pio_get_var( fid, vid, strt, cnt, wrk2d )
if(associated(wrk2d_in)) then
wrk2d_in = reshape( wrk2d(:,:),(/file%nlon,file%nlat/), order=order )
deallocate(wrk2d)
else
wrk2d_in => wrk2d
end if
j=1
do c=begchunk,endchunk
ncols = get_ncols_p(c)
call get_rlat_all_p(c, pcols, to_lats)
call get_rlon_all_p(c, pcols, to_lons)
call lininterp_init(file%lons, file%nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
call lininterp_init(file%lats, file%nlat, to_lats, ncols, 1, lat_wgts)
call lininterp(wrk2d_in, file%nlon, file%nlat, loc_arr(1:ncols,c-begchunk+1), ncols, lon_wgts, lat_wgts)
call lininterp_finish(lon_wgts)
call lininterp_finish(lat_wgts)
end do
if(allocated(wrk2d)) then
deallocate(wrk2d)
else
deallocate(wrk2d_in)
end if
if(dycore_is('LR')) call polar_average(loc_arr)
end subroutine read_2d_trc
!------------------------------------------------------------------------
subroutine read_za_trc( fid, vid, loc_arr, strt, cnt, file, order ) 1,12
use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
use phys_grid, only : pcols, begchunk, endchunk, get_ncols_p, get_rlat_all_p, get_rlon_all_p
use mo_constants, only : pi
use dycore, only : dycore_is
use polar_avg, only : polar_average
implicit none
type(file_desc_t), intent(in) :: fid
type(var_desc_t), intent(in) :: vid
integer, intent(in) :: strt(:), cnt(:), order(2)
real(r8),intent(out) :: loc_arr(:,:,:)
type (trfile), intent(in) :: file
type(interp_type) :: lat_wgts
real(r8) :: to_lats(pcols), to_lons(pcols), wrk(pcols)
real(r8), allocatable, target :: wrk2d(:,:)
real(r8), pointer :: wrk2d_in(:,:)
integer :: tsize, c, i, j, k, ierr, ncols
nullify(wrk2d_in)
allocate( wrk2d(cnt(1),cnt(2)), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'read_2d_trc: wrk2d allocation error = ',ierr
call endrun
end if
if(order(1)/=1 .or. order(2)/=2 .or. cnt(1)/=file%nlat .or. cnt(2)/=file%nlev) then
allocate( wrk2d_in(file%nlat, file%nlev), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'read_2d_trc: wrk2d_in allocation error = ',ierr
call endrun
end if
end if
ierr = pio_get_var( fid, vid, strt, cnt, wrk2d )
if(associated(wrk2d_in)) then
wrk2d_in = reshape( wrk2d(:,:),(/file%nlon,file%nlat/), order=order )
deallocate(wrk2d)
else
wrk2d_in => wrk2d
end if
j=1
do c=begchunk,endchunk
ncols = get_ncols_p(c)
call get_rlat_all_p(c, pcols, to_lats)
call lininterp_init(file%lats, file%nlat, to_lats, ncols, 1, lat_wgts)
do k=1,file%nlev
call lininterp(wrk2d_in(:,k), file%nlat, wrk(1:ncols), ncols, lat_wgts)
loc_arr(1:ncols,k,c-begchunk+1) = wrk(1:ncols)
end do
call lininterp_finish(lat_wgts)
end do
end subroutine read_za_trc
!------------------------------------------------------------------------
subroutine read_3d_trc( fid, vid, loc_arr, strt, cnt, file, order) 1,19
use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
use phys_grid, only : pcols, begchunk, endchunk, get_ncols_p, get_rlat_all_p, get_rlon_all_p
use mo_constants, only : pi
use dycore, only : dycore_is
use polar_avg, only : polar_average
use dycore, only : dycore_is
implicit none
type(file_desc_t), intent(in) :: fid
type(var_desc_t), intent(in) :: vid
integer, intent(in) :: strt(:), cnt(:), order(3)
real(r8),intent(out) :: loc_arr(:,:,:)
type (trfile), intent(in) :: file
integer :: i,j,k, astat, c, ncols
integer :: jlim(2), jl, ju, ierr
integer :: gndx
real(r8), allocatable, target :: wrk3d(:,:,:)
real(r8), pointer :: wrk3d_in(:,:,:)
real(r8) :: to_lons(pcols), to_lats(pcols)
real(r8), parameter :: zero=0_r8, twopi=2_r8*pi
type(interp_type) :: lon_wgts, lat_wgts
loc_arr(:,:,:) = 0._r8
nullify(wrk3d_in)
allocate(wrk3d(cnt(1),cnt(2),cnt(3)), stat=ierr)
if( ierr /= 0 ) then
write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr
call endrun
end if
ierr = pio_get_var( fid, vid, strt, cnt, wrk3d )
if(order(1)/=1 .or. order(2)/=2 .or. order(3)/=3 .or. &
cnt(1)/=file%nlon.or.cnt(2)/=file%nlat.or.cnt(3)/=file%nlev) then
allocate(wrk3d_in(file%nlon,file%nlat,file%nlev),stat=ierr)
if( ierr /= 0 ) then
write(iulog,*) 'read_3d_trc: wrk3d allocation error = ',ierr
call endrun
end if
wrk3d_in = reshape( wrk3d(:,:,:),(/file%nlon,file%nlat,file%nlev/), order=order )
deallocate(wrk3d)
else
wrk3d_in => wrk3d
end if
j=1
do c=begchunk,endchunk
ncols = get_ncols_p(c)
call get_rlat_all_p(c, pcols, to_lats)
call get_rlon_all_p(c, pcols, to_lons)
call lininterp_init(file%lons, file%nlon, to_lons(1:ncols), ncols, 2, lon_wgts, zero, twopi)
call lininterp_init(file%lats, file%nlat, to_lats(1:ncols), ncols, 1, lat_wgts)
call lininterp(wrk3d_in, file%nlon, file%nlat, file%nlev, loc_arr(:,:,c-begchunk+1), ncols, pcols, lon_wgts, lat_wgts)
call lininterp_finish(lon_wgts)
call lininterp_finish(lat_wgts)
end do
if(allocated(wrk3d)) then
deallocate( wrk3d, stat=astat )
else
deallocate( wrk3d_in, stat=astat )
end if
if( astat/= 0 ) then
write(iulog,*) 'read_3d_trc: failed to deallocate wrk3d array; error = ',astat
call endrun
endif
if(dycore_is('LR')) call polar_average(file%nlev, loc_arr)
end subroutine read_3d_trc
!------------------------------------------------------------------------------
subroutine interpolate_trcdata( flds, file, state ) 1,6
use mo_util, only : rebin
use physics_types,only : physics_state
use physconst, only : cday
implicit none
type (trfld),intent(inout) :: flds(:)
type (trfile), intent(inout) :: file
type(physics_state), intent(in):: state(begchunk:endchunk)
real(r8) :: fact1, fact2
real(r8) :: deltat
integer :: f,nflds,c,ncol, i,k
real(r8) :: ps(pcols)
real(r8) :: datain(pcols,file%nlev)
real(r8) :: pin(pcols,file%nlev)
real(r8) :: model_z(pverp)
real(r8), parameter :: m2km = 1.e-3_r8
nflds = size(flds)
if ( file%interp_recs == 4 ) then
deltat = file%datatimes(3) - file%datatimes(1)
fact1 = (file%datatimes(3) - file%datatimem)/deltat
fact2 = 1._r8-fact1
!$OMP PARALLEL DO PRIVATE (C, NCOL)
do c = begchunk,endchunk
ncol = state(c)%ncol
if ( file%has_ps ) then
file%ps_in(1)%data(:ncol,c) = fact1*file%ps_in(1)%data(:ncol,c) + fact2*file%ps_in(3)%data(:ncol,c)
endif
do f = 1,nflds
flds(f)%input(1)%data(:ncol,:,c) = fact1*flds(f)%input(1)%data(:ncol,:,c) + fact2*flds(f)%input(3)%data(:ncol,:,c)
enddo
enddo
deltat = file%datatimes(4) - file%datatimes(2)
fact1 = (file%datatimes(4) - file%datatimep)/deltat
fact2 = 1._r8-fact1
!$OMP PARALLEL DO PRIVATE (C, NCOL)
do c = begchunk,endchunk
ncol = state(c)%ncol
if ( file%has_ps ) then
file%ps_in(2)%data(:ncol,c) = fact1*file%ps_in(2)%data(:ncol,c) + fact2*file%ps_in(4)%data(:ncol,c)
endif
do f = 1,nflds
flds(f)%input(2)%data(:ncol,:,c) = fact1*flds(f)%input(2)%data(:ncol,:,c) + fact2*flds(f)%input(4)%data(:ncol,:,c)
enddo
enddo
endif
file%interp_recs = 2
deltat = file%datatimep - file%datatimem
if ( file%cyclical .and. (deltat < 0._r8) ) then
deltat = deltat+file%one_yr
if ( file%datatimep >= file%curr_mod_time ) then
fact1 = (file%datatimep - file%curr_mod_time)/deltat
else
fact1 = (file%datatimep+file%one_yr - file%curr_mod_time)/deltat
endif
else
fact1 = (file%datatimep - file%curr_mod_time)/deltat
endif
! this assures that FIXED data are b4b on restarts
if ( file%fixed ) then
fact1 = dble(int(fact1*cday+.5_r8))/dble(cday)
endif
fact2 = 1._r8-fact1
!$OMP PARALLEL DO PRIVATE (C, NCOL, PS, I, K, PIN, F, DATAIN, MODEL_Z)
do c = begchunk,endchunk
ncol = state(c)%ncol
if (file%surf_data) then
do f = 1,nflds
do i = 1,ncol
flds(f)%data(i,1,c+flds(f)%chnk_offset) = &
fact1*flds(f)%input(nm)%data(i,1,c) + fact2*flds(f)%input(np)%data(i,1,c)
enddo
enddo
elseif (file%alt_data) then
do f = 1,nflds
datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c) + fact2*flds(f)%input(np)%data(:ncol,:,c)
do i = 1,ncol
model_z(1:pverp) = m2km * state(c)%zi(i,pverp:1:-1)
call rebin( file%nlev, pver, file%ilevs, model_z, datain(i,:), flds(f)%data(i,:,c+flds(f)%chnk_offset) )
enddo
enddo
else
if ( file%has_ps ) then
ps(:ncol) = fact1*file%ps_in(nm)%data(:ncol,c) + fact2*file%ps_in(np)%data(:ncol,c)
do i = 1,ncol
do k = 1,file%nlev
pin(i,k) = file%p0*file%hyam(k) + ps(i)*file%hybm(k)
enddo
enddo
else
do k = 1,file%nlev
pin(:,k) = file%levs(k)
enddo
endif
do f = 1,nflds
datain(:ncol,:) = fact1*flds(f)%input(nm)%data(:ncol,:,c) + fact2*flds(f)%input(np)%data(:ncol,:,c)
if ( file%top_bndry ) then
call vert_interp_ub(ncol, file%nlev, file%levs, datain(:ncol,:), flds(f)%data(:ncol,:,c+flds(f)%chnk_offset) )
else
call vert_interp(ncol, file%nlev, pin, state(c)%pmid, datain, flds(f)%data(:,:,c+flds(f)%chnk_offset) )
endif
enddo
endif
enddo
end subroutine interpolate_trcdata
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine get_dimension( fid, dname, dsize, dimid, data ) 11,2
implicit none
type(file_desc_t), intent(in) :: fid
character(*), intent(in) :: dname
integer, intent(out) :: dsize
integer, optional, intent(out) :: dimid
real(r8), optional, pointer, dimension(:) :: data
integer :: vid, ierr, id
ierr = pio_inq_dimid( fid, dname, id )
ierr = pio_inq_dimlen( fid, id, dsize )
if ( present(dimid) ) then
dimid = id
endif
if ( present(data) ) then
if ( associated(data) ) then
deallocate(data, stat=ierr)
if( ierr /= 0 ) then
write(iulog,*) 'get_dimension: data deallocation error = ',ierr
call endrun('get_dimension: failed to deallocate data array')
end if
endif
allocate( data(dsize), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'get_dimension: data allocation error = ',ierr
call endrun('get_dimension: failed to allocate data array')
end if
ierr = pio_inq_varid( fid, dname, vid )
ierr = pio_get_var( fid, vid, data )
endif
end subroutine get_dimension
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine set_cycle_indices( fileid, cyc_ndx_beg, cyc_ndx_end, cyc_yr ) 2,4
implicit none
type(file_desc_t), intent(in) :: fileid
integer, intent(out) :: cyc_ndx_beg
integer, intent(out) :: cyc_ndx_end
integer, intent(in) :: cyc_yr
integer, allocatable , dimension(:) :: dates, datesecs
integer :: timesize, i, astat, year, ierr
type(var_desc_T) :: dateid
call get_dimension( fileid, 'time', timesize )
cyc_ndx_beg=-1
allocate( dates(timesize), stat=astat )
if( astat/= 0 ) then
write(*,*) 'set_cycle_indices: failed to allocate dates array; error = ',astat
call endrun
end if
ierr = pio_inq_varid( fileid, 'date', dateid )
ierr = pio_get_var( fileid, dateid, dates )
do i=1,timesize
year = dates(i) / 10000
if ( year == cyc_yr ) then
if (cyc_ndx_beg < 0) then
cyc_ndx_beg = i
endif
cyc_ndx_end = i
endif
enddo
deallocate( dates, stat=astat )
if( astat/= 0 ) then
write(*,*) 'set_cycle_indices: failed to deallocate dates array; error = ',astat
call endrun
end if
if (cyc_ndx_beg < 0) then
write(*,*) 'set_cycle_indices: cycle year not found : ' , cyc_yr
call endrun('set_cycle_indices: cycle year not found')
endif
end subroutine set_cycle_indices
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine open_trc_datafile( fname, path, piofile, times, cyc_ndx_beg, cyc_ndx_end, cyc_yr ) 3,13
use ioFileMod, only: getfil
use cam_pio_utils, only: cam_pio_openfile
implicit none
character(*), intent(in) :: fname
character(*), intent(in) :: path
type(file_desc_t), intent(inout) :: piofile
real(r8), pointer :: times(:)
integer, optional, intent(out) :: cyc_ndx_beg
integer, optional, intent(out) :: cyc_ndx_end
integer, optional, intent(in) :: cyc_yr
character(len=shr_kind_cl) :: filen, filepath
integer :: year, month, day, dsize, i, timesize
integer :: dateid,secid
integer, allocatable , dimension(:) :: dates, datesecs
integer :: astat, ierr
filepath = trim(path) // '/' // trim(fname)
!
! open file and get fileid
!
call getfil( filepath, filen, 0 )
call cam_pio_openfile( piofile, filen, PIO_NOWRITE)
if(masterproc) write(iulog,*)'open_trc_datafile: ',trim(filen)
call get_dimension(piofile, 'time', timesize)
if ( associated(times) ) then
deallocate(times, stat=ierr)
if( ierr /= 0 ) then
write(iulog,*) 'open_trc_datafile: data deallocation error = ',ierr
call endrun('open_trc_datafile: failed to deallocate data array')
end if
endif
allocate( times(timesize), stat=ierr )
if( ierr /= 0 ) then
write(iulog,*) 'open_trc_datafile: data allocation error = ',ierr
call endrun('open_trc_datafile: failed to allocate data array')
end if
allocate( dates(timesize), stat=astat )
if( astat/= 0 ) then
if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate dates array; error = ',astat
call endrun
end if
allocate( datesecs(timesize), stat=astat )
if( astat/= 0 ) then
if(masterproc) write(iulog,*) 'open_trc_datafile: failed to allocate datesec array; error = ',astat
call endrun
end if
ierr = pio_inq_varid( piofile, 'date', dateid )
call pio_seterrorhandling( piofile, PIO_BCAST_ERROR)
ierr = pio_inq_varid( piofile, 'datesec', secid )
call pio_seterrorhandling( piofile, PIO_INTERNAL_ERROR)
if(ierr==PIO_NOERR) then
ierr = pio_get_var( piofile, secid, datesecs )
else
datesecs=0
end if
ierr = pio_get_var( piofile, dateid, dates )
do i=1,timesize
year = dates(i) / 10000
month = mod(dates(i),10000)/100
day = mod(dates(i),100)
call set_time_float_from_date( times(i), year, month, day, datesecs(i) )
if ( present(cyc_yr) ) then
if ( year == cyc_yr ) then
if ( present(cyc_ndx_beg) .and. (cyc_ndx_beg < 0) ) then
cyc_ndx_beg = i
endif
if ( present(cyc_ndx_end) ) then
cyc_ndx_end = i
endif
endif
endif
enddo
deallocate( dates, stat=astat )
if( astat/= 0 ) then
if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate dates array; error = ',astat
call endrun
end if
deallocate( datesecs, stat=astat )
if( astat/= 0 ) then
if(masterproc) write(iulog,*) 'open_trc_datafile: failed to deallocate datesec array; error = ',astat
call endrun
end if
if ( present(cyc_yr) .and. present(cyc_ndx_beg) ) then
if (cyc_ndx_beg < 0) then
write(iulog,*) 'open_trc_datafile: cycle year not found : ' , cyc_yr
call endrun('open_trc_datafile: cycle year not found')
endif
endif
end subroutine open_trc_datafile
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
subroutine specify_fields( specifier, fields ) 1,2
implicit none
character(len=*), intent(in) :: specifier(:)
type(trfld), pointer, dimension(:) :: fields
integer :: fld_cnt, astat
integer :: i,j
character(len=shr_kind_cl) :: str1, str2
character(len=32), allocatable, dimension(:) :: fld_name, src_name
integer :: nflds
nflds = size(specifier)
allocate(fld_name(nflds), src_name(nflds), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'specify_fields: failed to allocate fld_name, src_name arrays; error = ',astat
call endrun
end if
fld_cnt = 0
count_cnst: do i = 1, nflds
if ( len_trim( specifier(i) ) == 0 ) then
exit count_cnst
endif
j = scan( specifier(i),':')
if (j > 0) then
str1 = trim(adjustl( specifier(i)(:j-1) ))
str2 = trim(adjustl( specifier(i)(j+1:) ))
fld_name(i) = trim(adjustl( str1 ))
src_name(i) = trim(adjustl( str2 ))
else
fld_name(i) = trim(adjustl( specifier(i) ))
src_name(i) = trim(adjustl( specifier(i) ))
endif
fld_cnt = fld_cnt + 1
enddo count_cnst
if( fld_cnt < 1 ) then
return
end if
!-----------------------------------------------------------------------
! ... allocate field type array
!-----------------------------------------------------------------------
allocate( fields(fld_cnt), stat=astat )
if( astat/= 0 ) then
write(iulog,*) 'specify_fields: failed to allocate fields array; error = ',astat
call endrun
end if
do i = 1,fld_cnt
fields(i)%fldnam = fld_name(i)
fields(i)%srcnam = src_name(i)
enddo
deallocate(fld_name, src_name)
end subroutine specify_fields
!------------------------------------------------------------------------------
subroutine init_trc_restart( whence, piofile, tr_file ) 4
implicit none
character(len=*), intent(in) :: whence
type(file_desc_t), intent(inout) :: piofile
type(trfile), intent(inout) :: tr_file
character(len=32) :: name
integer :: ioerr, mcdimid, maxlen
! Dimension should already be defined in restart file
call pio_seterrorhandling(pioFile, PIO_BCAST_ERROR)
ioerr = pio_inq_dimid(pioFile,'max_chars', mcdimid)
call pio_seterrorhandling(pioFile, PIO_INTERNAL_ERROR)
! but define it if nessasary
if(ioerr/= PIO_NOERR) then
ioerr = pio_def_dim(pioFile, 'max_chars', SHR_KIND_CL, mcdimid)
end if
if(len_trim(tr_file%curr_filename)>1) then
allocate(tr_file%currfnameid)
name = trim(whence)//'_curr_fname'
ioerr = pio_def_var(pioFile, name,pio_char, (/mcdimid/), tr_file%currfnameid)
ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'offset_time', tr_file%offset_time)
maxlen = len_trim(tr_file%curr_filename)
ioerr = pio_put_att(pioFile, tr_file%currfnameid, 'actual_len', maxlen)
else
nullify(tr_file%currfnameid)
end if
if(len_trim(tr_file%next_filename)>1) then
allocate(tr_file%nextfnameid)
name = trim(whence)//'_next_fname'
ioerr = pio_def_var(pioFile, name,pio_char, (/mcdimid/), tr_file%nextfnameid)
maxlen = len_trim(tr_file%next_filename)
ioerr = pio_put_att(pioFile, tr_file%nextfnameid, 'actual_len', maxlen)
else
nullify(tr_file%nextfnameid)
end if
end subroutine init_trc_restart
!-------------------------------------------------------------------------
! writes file names to restart file
!-------------------------------------------------------------------------
subroutine write_trc_restart( piofile, tr_file ) 4
implicit none
type(file_desc_t), intent(inout) :: piofile
type(trfile), intent(inout) :: tr_file
integer :: ioerr, slen ! error status
if(associated(tr_file%currfnameid)) then
ioerr = pio_put_var(pioFile, tr_file%currfnameid, tr_file%curr_filename)
deallocate(tr_file%currfnameid)
nullify(tr_file%currfnameid)
end if
if(associated(tr_file%nextfnameid)) then
ioerr = pio_put_var(pioFile, tr_file%nextfnameid, tr_file%next_filename)
deallocate(tr_file%nextfnameid)
nullify(tr_file%nextfnameid)
end if
end subroutine write_trc_restart
!-------------------------------------------------------------------------
! reads file names from restart file
!-------------------------------------------------------------------------
subroutine read_trc_restart( whence, piofile, tr_file ) 4
implicit none
character(len=*), intent(in) :: whence
type(file_desc_t), intent(inout) :: piofile
type(trfile), intent(inout) :: tr_file
type(var_desc_t) :: vdesc
character(len=64) :: name
integer :: ioerr ! error status
integer :: slen
call PIO_SetErrorHandling(piofile, PIO_BCAST_ERROR)
name = trim(whence)//'_curr_fname'
ioerr = pio_inq_varid(piofile, name, vdesc)
if(ioerr==PIO_NOERR) then
tr_file%curr_filename=' '
ioerr = pio_get_att(piofile, vdesc, 'offset_time', tr_file%offset_time)
ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen)
ioerr = pio_get_var(piofile, vdesc, tr_file%curr_filename)
if(slen<SHR_KIND_CL) tr_file%curr_filename(slen+1:)=' '
end if
name = trim(whence)//'_next_fname'
ioerr = pio_inq_varid(piofile, name, vdesc)
if(ioerr==PIO_NOERR) then
tr_file%next_filename=' '
ioerr = pio_get_att(piofile, vdesc, 'actual_len', slen)
ioerr = pio_get_var(piofile, vdesc, tr_file%next_filename)
if(slen<SHR_KIND_CL) tr_file%next_filename(slen+1:)=' '
end if
call PIO_SetErrorHandling(piofile, PIO_INTERNAL_ERROR)
end subroutine read_trc_restart
!------------------------------------------------------------------------------
subroutine vert_interp( ncol, levsiz, pin, pmid, datain, dataout) 1
!-----------------------------------------------------------------------
!
! Interpolate data from current time-interpolated values to model levels
!--------------------------------------------------------------------------
implicit none
! Arguments
!
integer, intent(in) :: ncol ! number of atmospheric columns
integer, intent(in) :: levsiz
real(r8), intent(in) :: pin(pcols,levsiz)
real(r8), intent(in) :: pmid(pcols,pver) ! level pressures
real(r8), intent(in) :: datain(pcols,levsiz)
real(r8), intent(out) :: dataout(pcols,pver)
!
! local storage
!
integer :: i ! longitude index
integer :: k, kk, kkstart ! level indices
integer :: kupper(pcols) ! Level indices for interpolation
real(r8) :: dpu ! upper level pressure difference
real(r8) :: dpl ! lower level pressure difference
!--------------------------------------------------------------------------
!
! Initialize index array
!
do i=1,ncol
kupper(i) = 1
end do
do k=1,pver
!
! Top level we need to start looking is the top level for the previous k
! for all column points
!
kkstart = levsiz
do i=1,ncol
kkstart = min0(kkstart,kupper(i))
end do
!
! Store level indices for interpolation
!
do kk=kkstart,levsiz-1
do i=1,ncol
if (pin(i,kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(i,kk+1)) then
kupper(i) = kk
end if
end do
end do
! interpolate or extrapolate...
do i=1,ncol
if (pmid(i,k) .lt. pin(i,1)) then
dataout(i,k) = datain(i,1)*pmid(i,k)/pin(i,1)
else if (pmid(i,k) .gt. pin(i,levsiz)) then
dataout(i,k) = datain(i,levsiz)
else
dpu = pmid(i,k) - pin(i,kupper(i))
dpl = pin(i,kupper(i)+1) - pmid(i,k)
dataout(i,k) = (datain(i,kupper(i) )*dpl + &
datain(i,kupper(i)+1)*dpu)/(dpl + dpu)
end if
end do
end do
end subroutine vert_interp
!------------------------------------------------------------------------------
subroutine vert_interp_ub( ncol, nlevs, plevs, datain, dataout ) 1,1
use hycoef, only : hyai, ps0
!-----------------------------------------------------------------------
!
! Interpolate data from current time-interpolated values to top interface pressure
! -- from mo_tgcm_ubc.F90
!--------------------------------------------------------------------------
implicit none
! Arguments
!
integer, intent(in) :: ncol
integer, intent(in) :: nlevs
real(r8), intent(in) :: plevs(nlevs)
real(r8), intent(in) :: datain(ncol,nlevs)
real(r8), intent(out) :: dataout(ncol)
!
! local variables
!
integer :: i,ku,kl,kk
real(r8) :: pinterp, delp
pinterp = ps0 * hyai(1)
if( pinterp <= plevs(1) ) then
kl = 1
ku = 1
delp = 0._r8
else if( pinterp >= plevs(nlevs) ) then
kl = nlevs
ku = nlevs
delp = 0._r8
else
do kk = 2,nlevs
if( pinterp <= plevs(kk) ) then
ku = kk
kl = kk - 1
delp = log( pinterp/plevs(kk) ) / log( plevs(kk-1)/plevs(kk) )
exit
end if
end do
end if
do i = 1,ncol
dataout(i) = datain(i,kl) + delp * (datain(i,ku) - datain(i,kl))
end do
end subroutine vert_interp_ub
!------------------------------------------------------------------------------
end module tracer_data