FIXME : Add documentation
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| complex, | intent(out), | dimension (-ntgrid:,:,g_lo%llim_proc:) | :: | g | ||
| real, | intent(in) | :: | scale | |||
| integer, | intent(out) | :: | istatus | |||
| logical, | intent(in) | :: | has_phi | |||
| logical, | intent(in) | :: | has_apar | |||
| logical, | intent(in) | :: | has_bpar | |||
| character(len=*), | intent(in), | optional | :: | fileopt |
subroutine gs2_restore_many (g, scale, istatus, has_phi, has_apar, has_bpar, fileopt)
!MR, 2007: restore kx_shift array if already allocated
# ifdef NETCDF
use mp, only: iproc, mp_abort, proc0, broadcast
use fields_arrays, only: phinew, aparnew, bparnew
use fields_arrays, only: phi, apar, bpar
use dist_fn_arrays, only: gexp_1, gexp_2, gexp_3, kx_shift, theta0_shift
use kt_grids, only: naky, ntheta0
use gs2_layouts, only: layout
# endif
use theta_grid, only: ntgrid
use gs2_layouts, only: g_lo
use file_utils, only: error_unit
use array_utils, only: zero_array
implicit none
complex, dimension (-ntgrid:,:,g_lo%llim_proc:), intent (out) :: g
logical :: this_proc_reads_fields, need_to_broadcast_fields
real, intent (in) :: scale
integer, intent (out) :: istatus
logical, intent (in) :: has_phi, has_apar, has_bpar
character (len=*), intent(in), optional :: fileopt
# ifdef NETCDF
!> This parameter controls whether or not gs2
!> aborts if it cannot read the restart file (and
!> it has been told to load g from a restart).
!> It is set to .true..
!> I can think of no conceivable reason why this
!> would ever need to be .false., but comments welcome. EGH
logical, parameter :: abort_on_restart_fail = .true.
integer, dimension(3) :: counts, start_pos
character (run_name_size) :: file_proc
character (len = 5) :: restart_layout
integer :: i, n_elements
real :: fac
logical :: has_explicit_source_terms
logical :: is_one_file_per_processor
logical :: has_gexp_1
logical, save :: explicit_warning_given = .false.
n_elements = max(g_lo%ulim_proc-g_lo%llim_proc+1, 0)
is_one_file_per_processor = read_many .or. (.not. has_netcdf_parallel)
! Decide if we are reading the fields
need_to_broadcast_fields = proc_to_save_fields >= 0
this_proc_reads_fields = (iproc == proc_to_save_fields) .or. (.not. need_to_broadcast_fields)
file_proc = get_file_proc(is_one_file_per_processor, fileopt)
if (is_one_file_per_processor) then
istatus = nf90_open (file_proc, NF90_NOWRITE, ncid)
else
# ifdef NETCDF_PARALLEL
! If using netcdf version 4.1.2 deleted NF90_MPIIO and the associated IOR
istatus = nf90_open (file_proc, IOR(NF90_NOWRITE, NF90_MPIIO), ncid, comm=mp_comm, info=mp_info)
# endif
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, file=file_proc)
call check_netcdf_file_precision (ncid)
istatus = nf90_inq_dimid (ncid, "theta", thetaid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='theta',&
abort=abort_on_restart_fail)
istatus = nf90_inq_dimid (ncid, "sign", signid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='sign',&
abort=abort_on_restart_fail)
istatus = nf90_inq_dimid (ncid, "glo", gloid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='glo',&
abort=abort_on_restart_fail)
istatus = nf90_inq_dimid (ncid, "aky", kyid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='aky',&
abort=abort_on_restart_fail)
istatus = nf90_inq_dimid (ncid, "akx", kxid)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, dim='akx',&
abort=abort_on_restart_fail)
istatus = nf90_inquire_dimension (ncid, thetaid, len=i)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=thetaid,&
abort=abort_on_restart_fail, dim='theta')
if (i /= 2*ntgrid + 1) write(*,*) 'Restart error: ntgrid=? ',i,' : ',ntgrid,' : ',iproc
istatus = nf90_inquire_dimension (ncid, signid, len=i)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=signid,&
abort=abort_on_restart_fail)
if (i /= 2) write(*,*) 'Restart error: sign=? ',i,' : ',iproc
istatus = nf90_inquire_dimension (ncid, gloid, len=i)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=gloid,&
abort=abort_on_restart_fail)
if (is_one_file_per_processor) then
if (i /= g_lo%ulim_proc-g_lo%llim_proc+1) write(*,*) 'Restart error: glo=? ',i,' : ',iproc
else
if (i /= g_lo%ulim_world+1) write(*,*) 'Restart error: glo=? ',i,' : ',iproc
endif
istatus = nf90_inquire_dimension (ncid, kyid, len=i)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=kyid,&
abort=abort_on_restart_fail)
if (i /= naky) write(*,*) 'Restart error: naky=? ',i,' : ',naky,' : ',iproc
istatus = nf90_inquire_dimension (ncid, kxid, len=i)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, dimid=kxid,&
abort=abort_on_restart_fail)
if (i /= ntheta0) write(*,*) 'Restart error: ntheta0=? ',i,' : ',ntheta0,' : ',iproc
if (this_proc_reads_fields) then
if (has_phi) then
istatus = nf90_inq_varid (ncid, "phi_r", phir_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='phi_r')
istatus = nf90_inq_varid (ncid, "phi_i", phii_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='phi_i')
end if
if (has_apar) then
istatus = nf90_inq_varid (ncid, "apar_r", aparr_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='apar_r')
istatus = nf90_inq_varid (ncid, "apar_i", apari_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='apar_i')
end if
if (has_bpar) then
istatus = nf90_inq_varid (ncid, "bpar_r", bparr_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='bpar_r')
istatus = nf90_inq_varid (ncid, "bpar_i", bpari_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='bpar_i')
end if
end if
if (allocated(kx_shift)) then ! MR begin
istatus = nf90_inq_varid (ncid, "kx_shift", kx_shift_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='kx_shift')
endif ! MR end
if (allocated(theta0_shift)) then
istatus = nf90_inq_varid (ncid, "theta0_shift", theta0_shift_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='theta0_shift')
endif
if (include_explicit_source_in_restart) then
has_explicit_source_terms = .true.
#ifndef SHMEM
has_gexp_1 = allocated(gexp_1)
#else
has_gexp_1 = associated(gexp_1)
#endif
if (has_gexp_1) then
istatus = nf90_inq_varid (ncid, "gexp_1_r", gexp_1_r_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
istatus = nf90_inq_varid (ncid, "gexp_1_i", gexp_1_i_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
endif
if (allocated(gexp_2)) then
istatus = nf90_inq_varid (ncid, "gexp_2_r", gexp_2_r_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
istatus = nf90_inq_varid (ncid, "gexp_2_i", gexp_2_i_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
endif
if (allocated(gexp_3)) then
istatus = nf90_inq_varid (ncid, "gexp_3_r", gexp_3_r_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
istatus = nf90_inq_varid (ncid, "gexp_3_i", gexp_3_i_id)
if (istatus /= NF90_NOERR) has_explicit_source_terms = .false.
endif
if ((.not. explicit_warning_given).and. proc0 .and. .not.has_explicit_source_terms) then
write(*, '("Warning: At least one explicit source term absent in restart file.")')
write(*, '(" Will not load these terms. This is probably OK unless you expect the restart file to contain these terms.")')
write(*, '(" This warning will not be repeated.")')
! Update this flag to ensure we don't give this warning again
explicit_warning_given = .true.
end if
else
! Always skip trying to get the explicit source terms if we're not
! including them
has_explicit_source_terms = .false.
end if
istatus = nf90_inq_varid (ncid, "gr", gr_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='gr',&
abort=abort_on_restart_fail)
istatus = nf90_inq_varid (ncid, "gi", gi_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='gi',&
abort=abort_on_restart_fail)
istatus = nf90_inq_varid (ncid, "layout", layout_id)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, var='layout_r')
istatus = nf90_get_var (ncid, layout_id, restart_layout)
if (istatus /= NF90_NOERR) then
call netcdf_error (istatus, ncid, layout_id, message=' in gs2_restore_many')
restart_layout = layout
endif
! Abort if the layouts don't match
if(restart_layout /= layout) then
call mp_abort("Incompatible layouts in restart file ("//restart_layout//") and simulation ("//layout//")",.true.)
end if
if (.not. allocated(tmpr)) &
allocate (tmpr(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc))
if (.not. allocated(tmpi)) &
allocate (tmpi(2*ntgrid+1,2,g_lo%llim_proc:g_lo%ulim_alloc))
call zero_array(tmpr); call zero_array(tmpi)
if (is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gr_id, tmpr)
else
start_pos = (/1,1,g_lo%llim_proc+1/)
counts = (/2*ntgrid+1, 2, n_elements/)
istatus = nf90_get_var (ncid, gr_id, tmpr, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gr_id,&
abort=abort_on_restart_fail, var='gr')
if(is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gi_id, tmpi)
else
istatus = nf90_get_var (ncid, gi_id, tmpi, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gi_id,&
abort=abort_on_restart_fail, var='gi')
g = cmplx(tmpr, tmpi)
if (include_explicit_source_in_restart .and. has_explicit_source_terms) then
! Explicit source term
#ifndef SHMEM
has_gexp_1 = allocated(gexp_1)
#else
has_gexp_1 = associated(gexp_1)
#endif
if (has_gexp_1) then
if (is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_1_r_id, tmpr)
else
start_pos = (/1,1,g_lo%llim_proc+1/)
counts = (/2*ntgrid+1, 2, n_elements/)
istatus = nf90_get_var (ncid, gexp_1_r_id, tmpr, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_1_r_id,&
var='gexp_1_r')
if(is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_1_i_id, tmpi)
else
istatus = nf90_get_var (ncid, gexp_1_i_id, tmpi, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_1_i_id,&
var='gexp_1_i')
gexp_1 = cmplx(tmpr, tmpi)
endif
! Explicit source term
if(allocated(gexp_2)) then
if (is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_2_r_id, tmpr)
else
start_pos = (/1,1,g_lo%llim_proc+1/)
counts = (/2*ntgrid+1, 2, n_elements/)
istatus = nf90_get_var (ncid, gexp_2_r_id, tmpr, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_2_r_id,&
var='gexp_2_r')
if(is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_2_i_id, tmpi)
else
istatus = nf90_get_var (ncid, gexp_2_i_id, tmpi, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_2_i_id,&
var='gexp_2_i')
gexp_2 = cmplx(tmpr, tmpi)
endif
! Explicit source term
if(allocated(gexp_3)) then
if (is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_3_r_id, tmpr)
else
start_pos = (/1,1,g_lo%llim_proc+1/)
counts = (/2*ntgrid+1, 2, n_elements/)
istatus = nf90_get_var (ncid, gexp_3_r_id, tmpr, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_3_r_id,&
var='gexp_3_r')
if (is_one_file_per_processor) then
istatus = nf90_get_var (ncid, gexp_3_i_id, tmpi)
else
istatus = nf90_get_var (ncid, gexp_3_i_id, tmpi, start=start_pos, count=counts)
end if
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, gexp_3_i_id,&
var='gexp_3_i')
gexp_3 = cmplx(tmpr, tmpi)
endif
end if
if (allocated(kx_shift)) then ! MR begin
if (.not. allocated(stmp)) allocate (stmp(naky)) ! MR
istatus = nf90_get_var (ncid, kx_shift_id, stmp)
if (istatus /= NF90_NOERR) &
call netcdf_error (istatus, ncid, kx_shift_id, var='kx_shift')
kx_shift = stmp
endif ! MR end
if (allocated(theta0_shift)) then
if (.not. allocated(stmp)) allocate (stmp(naky))
istatus = nf90_get_var (ncid, theta0_shift_id, stmp)
if (istatus /= NF90_NOERR) &
call netcdf_error (istatus, ncid, theta0_shift_id, var='theta0_shift')
theta0_shift = stmp
endif
if (this_proc_reads_fields) then
if (.not. allocated(ftmpr)) allocate (ftmpr(2*ntgrid+1,ntheta0,naky))
if (.not. allocated(ftmpi)) allocate (ftmpi(2*ntgrid+1,ntheta0,naky))
if (has_phi) then
istatus = nf90_get_var (ncid, phir_id, ftmpr)
if (istatus /= NF90_NOERR) &
call netcdf_error (istatus, ncid, phir_id, var= 'phir')
istatus = nf90_get_var (ncid, phii_id, ftmpi)
if (istatus /= NF90_NOERR) &
call netcdf_error (istatus, ncid, phii_id, var='phii')
phinew = cmplx(ftmpr, ftmpi)
end if
if (has_apar) then
istatus = nf90_get_var (ncid, aparr_id, ftmpr)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, aparr_id)
istatus = nf90_get_var (ncid, apari_id, ftmpi)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, apari_id)
aparnew = cmplx(ftmpr, ftmpi)
end if
if (has_bpar) then
istatus = nf90_get_var (ncid, bparr_id, ftmpr)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, bparr_id)
istatus = nf90_get_var (ncid, bpari_id, ftmpi)
if (istatus /= NF90_NOERR) call netcdf_error (istatus, ncid, bpari_id)
bparnew = cmplx(ftmpr, ftmpi)
end if
end if
if (need_to_broadcast_fields) then
if (has_phi) call broadcast(phinew, proc_to_save_fields)
if (has_apar) call broadcast(aparnew, proc_to_save_fields)
if (has_bpar) call broadcast(bparnew, proc_to_save_fields)
end if
if (has_phi) call zero_array(phi)
if (has_apar) call zero_array(apar)
if (has_bpar) call zero_array(bpar)
if (scale > 0.) then
g = g*scale
phinew = phinew*scale
aparnew = aparnew*scale
bparnew = bparnew*scale
else
fac = - scale/(maxval(abs(phinew)))
g = g*fac
phinew = phinew*fac
aparnew = aparnew*fac
bparnew = bparnew*fac
end if
! RN 2008/05/23: this was commented out. why? HJL 2013/05/15 Because it stops future writing to the file, it's now back in after setting initialized to false
istatus = nf90_close (ncid)
if (istatus /= NF90_NOERR) then
write(error_unit(),*) "nf90_close error: ", nf90_strerror(istatus),' ',iproc
end if
! Deallocate modul level arrays
call deallocate_arrays
# else
write (error_unit(),*) &
'ERROR: gs2_restore_many is called without netcdf'
# endif
end subroutine gs2_restore_many