problem with put_att (maximum number of attributes)
Alex Hoffmann
ah519 at cam.ac.uk
Mon Feb 23 05:57:05 CST 2009
Dear Wei-keng,
below is the script section in collective data mode:
-------
! this should have been possible in collective data mode, but did
not seem to work !!!
! call handle_err_nc_p(nfmpi_begin_indep_data(NCIDP),9)
! call
handle_err_nc_p(nfmpi_put_vara_real(NCIDP,array_varIDP,(/istart,jstart/),(/nx_mpi,ny_mpi/),reshape(array,(/nx_mpi,ny_mpi/))),10)
call
handle_err_nc_p(nfmpi_put_vara_real_all(NCIDP,array_varIDP,(/istart,jstart/),(/nx_mpi,ny_mpi/),reshape(array,(/nx_mpi,ny_mpi/))),10)
! call handle_err_nc_p(nfmpi_end_indep_data(NCIDP),11)
call handle_err_nc_p(nfmpi_close(NCIDP),12)
-------
Basically, I replaced the independent call by the collective call for
all processes, and commented out the call into/out of independent data
mode. This produces MPI error messages...
Cheers,
Lex
Wei-keng Liao wrote:
> Hi, Lex,
>
> Have you tried nfmpi_put_vara_real_all() in collective mode?
> Collective mode is the default mode.
>
> If yes, can you show me how you call it?
>
> Wei-keng
>
>
> On Feb 20, 2009, at 1:24 PM, Alex Hoffmann wrote:
>
>> Thank you for the clarifications.
>>
>> I included the mpi setup and netCDF output test code and the
>> corresponding Makefile as attachments. Maybe it can be of use as a
>> simple example.
>>
>> I am currently running my code with the bugs in the output file fixed,
>> and it seems to be working just fine.
>>
>> Thanks again,
>>
>> Lex
>>
>> Wei-keng Liao wrote:
>>> Hi, Lex, Please read below.
>>>
>>> On Feb 20, 2009, at 5:13 AM, Alex Hoffmann wrote:
>>>
>>>> Many thanks, Wei-keng, the issue seems to have been tracked down...
>>>>
>>>> For reasons that are still not completely clear to me, using trim() for
>>>> trimming the coordinate variable names and units (lines 480/81, 485/86)
>>>> seemed to leave a strange character (such as y diacritic or Y accent)
>>>> trailing on e.g. the units (such as degrees_east). Though this did not
>>>> affect the run on its own, any subsequent attribute definition made the
>>>> code crash. I fixed this by removing the trim().
>>>
>>> If you called trim() on a string argument in a function, the tailing
>>> character '\0' will not be passed to the function. So if the string
>>> length
>>> argument has a larger value, the function will include undefined tailing
>>> characters as part of the string. Glad to know you solve the problem.
>>>
>>>
>>>> If I may, I'd like to clarify one further issue with PnetCDF.
>>>> I initially implemented the code (also a small test program) in
>>>> collective data mode. This did not work (the different
>>>> processes/Cartesian sub domains) did not write their data into the
>>>> netCDF file. I then changed the data write statements into independent
>>>> data mode, with which it worked immediately. Am I right to use
>>>> independent data mode for this purpose?
>>>
>>> I am not sure what you meant here.
>>> Can you post your small test program to the list?
>>>
>>> Wei-keng
>>>
>>>>
>>>>
>>>> Thank you very much for your support.
>>>>
>>>> Cheers
>>>>
>>>> Lex
>>>>
>>>> Wei-keng Liao wrote:
>>>>> Hi, Lex,
>>>>>
>>>>> If this error message "NC definations on multiprocesses conflict"
>>>>> is the only one you got before every crash, it is most likely the
>>>>> attribute values are not the same across all processes.
>>>>> Please check that first, so you can rule out that possibility.
>>>>>
>>>>> Also, can you check the attributes in the output file for the case
>>>>> that ran successfully? (using ncmpidump -c filename) This is just
>>>>> to see if the attributes are the ones you expected. Also, please
>>>>> make sure the number of attributes created is less than 4K global
>>>>> or per variable, as Rob indicated.
>>>>>
>>>>> Wei-keng
>>>>>
>>>>>
>>>>> On Feb 19, 2009, at 5:53 PM, Alex Hoffmann wrote:
>>>>>
>>>>>> Dear Wei-keng,
>>>>>>
>>>>>> I understand that the attribute define arguments must be the same
>>>>>> across
>>>>>> all processes. I am also almost certain that the arguments you
>>>>>> cited
>>>>>> are identical across processes but I will perform a check.
>>>>>>
>>>>>> Be this as it may, I get the same arrow message if I uncomment one of
>>>>>> the text attributes specifying the unit or one of the real attributes
>>>>>> specifying the fill value, for instance. There is no possibility
>>>>>> that
>>>>>> these are different across processes as they are introduced straight
>>>>>> into the function call. Any check I performed seemed to indicate
>>>>>> that
>>>>>> the problem arises from the number of arguments I defined, and not
>>>>>> from
>>>>>> the arguments themselves, since they work all fine individually.
>>>>>>
>>>>>> I was actually thinking to write the attributes ad-hoc once the
>>>>>> dataset
>>>>>> is defined using regular netCDF functions. This however seemed to
>>>>>> be a
>>>>>> little elegant solution defeating the purpose of using PnetCDF.
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> Lex
>>>>>>
>>>>>>
>>>>>> Wei-keng Liao a écrit :
>>>>>>> Hi, Lex,
>>>>>>>
>>>>>>> PnetCDF requires the attribute define arguments be the same across
>>>>>>> all processes. This is why you got the error message for defining
>>>>>>> attribute "history" in line 400, if its value "conflicts" among
>>>>>>> processes.
>>>>>>> I think you have got this.
>>>>>>>
>>>>>>> I checked your codes. Starting from line 562, you have a few text
>>>>>>> attributes
>>>>>>> from des_tgas_picture(i), des_trac_picture(i), des_tpas_picture(i),
>>>>>>> and des.
>>>>>>> Similar for out_type == 2. Are you sure their values are all the
>>>>>>> same
>>>>>>> across
>>>>>>> all the processes? May be you can use MPI_Bcast to check if
>>>>>>> anyone has
>>>>>>> different value.
>>>>>>>
>>>>>>> Wei-keng
>>>>>>>
>>>>>>>
>>>>>>> On Feb 19, 2009, at 11:14 AM, Alex Hoffmann wrote:
>>>>>>>
>>>>>>>> Dear Rob
>>>>>>>> Thanks a lot for your answer.
>>>>>>>> Isolating the PnetCDF output module from the code as a stand-alone
>>>>>>>> would
>>>>>>>> require a substantial effort, and I do not encounter this error
>>>>>>>> in my
>>>>>>>> small test programs. I attached the output module for reference.
>>>>>>>>
>>>>>>>> The attribute section that causes trouble is between lines 369-405
>>>>>>>> (global) and 468-613 (variables). I am aware that the attribute
>>>>>>>> history
>>>>>>>> (line 400) may be problematic as every process may have a slightly
>>>>>>>> different value (for the time element taken from the fortran
>>>>>>>> function
>>>>>>>> date_and_time). Yet this is not the only issue.
>>>>>>>>
>>>>>>>> If I prematurely exit at line 492, uncommenting any further
>>>>>>>> attribute
>>>>>>>> definition between 468 and 489 causes a crash with the
>>>>>>>> multiprocesses
>>>>>>>> conflict.
>>>>>>>>
>>>>>>>> Note also that the attribute definitions between 558 and 613 only
>>>>>>>> produce no crash currently because of the print statements on
>>>>>>>> lines 560
>>>>>>>> and 591, which made me think that maybe there is an optimization
>>>>>>>> problem. What is strange is that if these are not included, the
>>>>>>>> crash
>>>>>>>> occurs only when I start writing the first file (movie file,
>>>>>>>> out_type=2), AND if I leave the first block (if out_type==1)
>>>>>>>> uncommented, although this block is not even entered in that
>>>>>>>> particular
>>>>>>>> case, and vice versa. Hope this is more or less clear...
>>>>>>>>
>>>>>>>> Cheers,
>>>>>>>> Lex
>>>>>>>>
>>>>>>>>
>>>>>>>> Robert Latham wrote:
>>>>>>>>> On Thu, Feb 19, 2009 at 10:09:37AM +0000, Alex Hoffmann wrote:
>>>>>>>>>> Hello,
>>>>>>>>>> I am getting the same error message as had already been inquired
>>>>>>>>>> on in
>>>>>>>>>> this list, i.e.:
>>>>>>>>>>
>>>>>>>>>> NC definations on multiprocesses conflict.
>>>>>>>>>> NOTE: Definitions across all process
>>>>>>>>>>
>>>>>>>>>> and this seems to be related somehow randomly to the total
>>>>>>>>>> number of
>>>>>>>>>> attributes (global or for each variable) that I intend to write
>>>>>>>>>> during
>>>>>>>>>> the the definition phase (and, though I am not sure about
>>>>>>>>>> this, to
>>>>>>>>>> the
>>>>>>>>>> number of processes used to run my F90 parallel code).
>>>>>>>>>>
>>>>>>>>>> More precisely, if I uncomment one put_att to replace another,
>>>>>>>>>> the
>>>>>>>>>> program still runs, if I do not comment out another one, the
>>>>>>>>>> program
>>>>>>>>>> crashes with the above error message.
>>>>>>>>>>
>>>>>>>>>> Is there any way of identifying where this problem comes from
>>>>>>>>>> or is
>>>>>>>>>> there a maximum amount of attributes (though I don't write an
>>>>>>>>>> excessive
>>>>>>>>>> amount of them)? There is no problem running it on one
>>>>>>>>>> process(or).
>>>>>>>>>
>>>>>>>>> Maybe you can show us with a small test program how you are
>>>>>>>>> encountering this error? The maximum number of attributes is 4k
>>>>>>>>> (the
>>>>>>>>> constant NC_MAX_ATTRS).
>>>>>>>>>
>>>>>>>>> ==rob
>>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> ______________________________________________________________________
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> Alex HOFFMANN PhD Candidate
>>>>>>>> Centre for Atmospheric Science
>>>>>>>> Department of Geography University of Cambridge
>>>>>>>> Downing Place, Cambridge CB2 3EN, UK
>>>>>>>> e-mail: ah519 at cam.ac.uk tel: +44 (0)1223 766581
>>>>>>>> www.geog.cam.ac.uk/people/hoffmann/
>>>>>>>> ______________________________________________________________________
>>>>>>>>
>>>>>>>>
>>>>>>>> !-*- F90 -*-
>>>>>>>> module atham_pnetcdf
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------------!
>>>>>>>>
>>>>>>>>
>>>>>>>> ! author: Alex
>>>>>>>> HOFFMANN !
>>>>>>>> ! email:
>>>>>>>> ah519 at cam.ac.uk !
>>>>>>>> ! date: February
>>>>>>>> 2009 !
>>>>>>>> ! version:
>>>>>>>> v0.1 !
>>>>>>>> !
>>>>>>>> !
>>>>>>>> ! collection of routines for ATHAM PnetCDF
>>>>>>>> output !
>>>>>>>> !--------------------------------------------------------------------!
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> ! module is to replace existing atham_netcdf.mod (GL) to write
>>>>>>>> ATHAM
>>>>>>>> output as netCDF files
>>>>>>>> ! same as atham_netcdf.mod (AH) but adapted to run with the
>>>>>>>> Parallel
>>>>>>>> NetCDF API,
>>>>>>>> ! developed by Northwestern University and Argonne National
>>>>>>>> Laboratory,
>>>>>>>> ! which is built on top of MPI-IO (tested with MPICH 2.1)
>>>>>>>>
>>>>>>>> use precision, only : kint, kreal, real4, my_real
>>>>>>>>
>>>>>>>> use atham_module, only : nx, ny, nz, ntx, nty, npx, npy
>>>>>>>>
>>>>>>>>
>>>>>>>> use atham_module, only : ntgas_picture, ntrac_picture,
>>>>>>>> ntpas_picture, &
>>>>>>>> ntgas_movie, ntrac_movie, ntpas_movie
>>>>>>>>
>>>>>>>> implicit none
>>>>>>>>
>>>>>>>> private
>>>>>>>> public :: netcdf_data
>>>>>>>>
>>>>>>>> ! 'pnetcdf.inc' file includes cpp directives that were not
>>>>>>>> understood
>>>>>>>> at compile-time
>>>>>>>> ! this is to be changed and hard-wired link to be replaced by the
>>>>>>>> include 'pnetcdf.inc' command
>>>>>>>> !include 'pnetcdf.inc'
>>>>>>>> #include "/usr/local/parallel-netcdf-1.0.1/include/pnetcdf.inc"
>>>>>>>> include 'mpif.h'
>>>>>>>>
>>>>>>>>
>>>>>>>> ! ######### tried saving netCDF variable IDs here so that they must
>>>>>>>> not be retrieved by inquiry
>>>>>>>> ! ######### if this creates conflict, especially between picture
>>>>>>>> and
>>>>>>>> movie file, this must be changed
>>>>>>>> integer(kint), save :: x_dimID, y_dimID,
>>>>>>>> z_dimID, t_dimID
>>>>>>>> integer(kint), save :: x_varID, y_varID,
>>>>>>>> z_varID, t_varID
>>>>>>>> integer(kint), save :: u_varID, v_varID, w_varID
>>>>>>>> integer(kint), save :: pnew_varID, tetnew_varID,
>>>>>>>> tempnew_varID, density_varID
>>>>>>>> integer(kint), save :: turbhor_varID,
>>>>>>>> turbver_varID, turblen_varID
>>>>>>>> integer(kint), dimension(99), save :: tgas_varID, trac_varID,
>>>>>>>> tpas_varID, d_trac_varID
>>>>>>>>
>>>>>>>>
>>>>>>>> integer(MPI_OFFSET_KIND), save :: tcounter
>>>>>>>>
>>>>>>>> real(real4), parameter :: undef = -0.99e33_real4 !
>>>>>>>> make sure precision of _FillValue and output (NF90_FLOAT) is the
>>>>>>>> same
>>>>>>>> !!!
>>>>>>>>
>>>>>>>> contains
>>>>>>>> !=====================================================================
>>>>>>>>
>>>>>>>>
>>>>>>>> subroutine netcdf_data(outfile,out_type,timetot)
>>>>>>>> ! IF: outfile = picture_file or movie_file (path/filename for
>>>>>>>> output created in atham_module)
>>>>>>>> ! out_type = 1 (picture_file) or 2 (movie_file)
>>>>>>>> ! timetot = total simulation time since simulation start
>>>>>>>>
>>>>>>>> use atham_module, only : unew, vnew, wnew, pnew, tetnew, tempnew, &
>>>>>>>> density, turbhor, turbver, turblen
>>>>>>>>
>>>>>>>> use atham_module, only : tgasnew, tracnew, tpasnew, tracdep
>>>>>>>>
>>>>>>>> use atham_module, only : itgas_picture, itrac_picture,
>>>>>>>> itpas_picture, &
>>>>>>>> itgas_movie, itrac_movie, itpas_movie
>>>>>>>>
>>>>>>>> use atham_module, only : mycoor, my_cart, myrank
>>>>>>>>
>>>>>>>> character(*), intent(in ) :: outfile
>>>>>>>> integer(kint), intent(in ) :: out_type
>>>>>>>> real(kreal), intent(in ) :: timetot
>>>>>>>>
>>>>>>>> !------------------------------------------------------------------!
>>>>>>>>
>>>>>>>> ! local
>>>>>>>> variables !
>>>>>>>> !------------------------------------------------------------------!
>>>>>>>>
>>>>>>>> integer(MPI_OFFSET_KIND) :: nx_mpi,ny_mpi, nz_mpi,
>>>>>>>> y_cl_mpi
>>>>>>>> integer(MPI_OFFSET_KIND) :: istart,jstart
>>>>>>>> integer(kint) :: STATUS
>>>>>>>> character(len=50) :: OP ! debugging variable
>>>>>>>> corresponding to the called netCDF OPeration
>>>>>>>> integer(kint) :: NCID
>>>>>>>> integer(kint) :: i
>>>>>>>> logical :: DynMovie_flag,
>>>>>>>> bounds_control
>>>>>>>> integer(kint) :: y_st, y_sp, y_cl
>>>>>>>>
>>>>>>>> ! replace DynMovie precompiler flag by program flag for more
>>>>>>>> convenient code structure
>>>>>>>> DynMovie_flag = .false.
>>>>>>>> #ifdef DynMovie
>>>>>>>> DynMovie_flag = .true.
>>>>>>>> #endif
>>>>>>>>
>>>>>>>> ! re-assign dimension lengths as an MPI_Offset type (necessary for
>>>>>>>> the PnetCDF interface)
>>>>>>>> ! cut-off subdomain borders
>>>>>>>> nx_mpi = nx-2
>>>>>>>> ny_mpi = ny-2
>>>>>>>> nz_mpi = nz
>>>>>>>>
>>>>>>>> ! assign coordinates of first element of each process in the
>>>>>>>> netCDF file
>>>>>>>> istart = mycoor(1)*nx_mpi + 1_MPI_OFFSET_KIND
>>>>>>>> jstart = mycoor(2)*ny_mpi + 1_MPI_OFFSET_KIND
>>>>>>>>
>>>>>>>> bounds_control = .true. ! debugging flag
>>>>>>>> tcounter = tcounter+1_MPI_OFFSET_KIND
>>>>>>>>
>>>>>>>> ! === open netCDF dataset: enter data mode
>>>>>>>> ------------------------!
>>>>>>>> ! if dataset does not yet exist (STATUS = ), create file and
>>>>>>>> enter data mode
>>>>>>>> ! else continue
>>>>>>>> OP = 'open existing netCDF file'
>>>>>>>> STATUS = nfmpi_open(my_cart, trim(outfile), NF_WRITE,
>>>>>>>> MPI_INFO_NULL, NCID)
>>>>>>>> if ( STATUS /= nf_noerr ) then
>>>>>>>> if (myrank==0) then
>>>>>>>> print *, ""
>>>>>>>> print *, "P-NetCDF operation: ",trim(OP)
>>>>>>>> print *, trim(nfmpi_strerror(STATUS))
>>>>>>>> print *, "atham_pnetcdf netcdf_data(): new netCDF file
>>>>>>>> created: ", outfile
>>>>>>>> print *, ""
>>>>>>>> end if
>>>>>>>> call netcdf_define(outfile,out_type,timetot,NCID)
>>>>>>>> tcounter = 1_MPI_OFFSET_KIND ! reset time counter for output if
>>>>>>>> new file is created
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! --- write data values
>>>>>>>> --------------------!
>>>>>>>> ! unlimited record dimension time (only one process writes the
>>>>>>>> whole dimension, i.e. accross all processes)
>>>>>>>> OP = 'write data values - unlimited dimension t'
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_begin_indep_data(NCID))
>>>>>>>>
>>>>>>>> if (myrank==0) then
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID, t_varID,
>>>>>>>> (/tcounter/), (/1_MPI_OFFSET_KIND/), real(timetot,real4)),OP)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_end_indep_data(NCID))
>>>>>>>>
>>>>>>>> ! --- write data values
>>>>>>>> --------------------!
>>>>>>>> ! discriminate between 2D and 3D case for output generation
>>>>>>>> if (nty>4) then
>>>>>>>> y_st = 2 ! y index start
>>>>>>>> y_sp = ny-1 ! y index stop
>>>>>>>> y_cl = ny-2 ! y number of columns
>>>>>>>> else
>>>>>>>> y_st = ny-1 ! 3 in 2D case
>>>>>>>> y_sp = ny-1 ! 3 in 2D case
>>>>>>>> y_cl = 1
>>>>>>>> end if
>>>>>>>> y_cl_mpi = y_cl
>>>>>>>>
>>>>>>>> ! dynamics data
>>>>>>>> OP = 'write data values - data'
>>>>>>>> if (out_type == 1 .OR. DynMovie_flag == .true.) then
>>>>>>>> call write_data(real(unew,real4),u_varID)
>>>>>>>> bounds_control = .false.
>>>>>>>> call write_data(real(vnew, real4),v_varID)
>>>>>>>> call write_data(real(wnew, real4),w_varID)
>>>>>>>> call write_data(real(pnew, real4),pnew_varID)
>>>>>>>> call write_data(real(tetnew, real4),tetnew_varID)
>>>>>>>> call write_data(real(tempnew,real4),tempnew_varID)
>>>>>>>> call write_data(real(density,real4),density_varID)
>>>>>>>> call write_data(real(turbhor,real4),turbhor_varID)
>>>>>>>> call write_data(real(turbver,real4),turbver_varID)
>>>>>>>> call write_data(real(turblen,real4),turblen_varID)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! tracer data
>>>>>>>> if (out_type == 1) then
>>>>>>>> do i=1,ntgas_picture
>>>>>>>> call
>>>>>>>> write_data(real(tgasnew(:,:,:,itgas_picture(i)),real4),tgas_varID(i))
>>>>>>>>
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> call
>>>>>>>> write_data(real(tracnew(:,:,:,itrac_picture(i)),real4),trac_varID(i))
>>>>>>>>
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_picture
>>>>>>>> call
>>>>>>>> write_data(real(tpasnew(:,:,:,itpas_picture(i)),real4),tpas_varID(i))
>>>>>>>>
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> call handle_err(nfmpi_begin_indep_data(NCID))
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID,
>>>>>>>> d_trac_varID(i),(/istart,jstart,tcounter/), &
>>>>>>>> (/nx_mpi,y_cl_mpi,1_MPI_OFFSET_KIND/),
>>>>>>>> real(reshape(tracdep(2:nx-1,y_st:y_sp,itrac_picture(i)),shape=(/nx-2,y_cl,1/)),real4)
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> ),OP)
>>>>>>>> call handle_err(nfmpi_end_indep_data(NCID))
>>>>>>>> enddo
>>>>>>>> else if (out_type == 2) then
>>>>>>>> do i=1,ntgas_movie
>>>>>>>> call
>>>>>>>> write_data(real(tgasnew(:,:,:,itgas_movie(i)),real4),tgas_varID(i))
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_movie
>>>>>>>> call
>>>>>>>> write_data(real(tracnew(:,:,:,itrac_movie(i)),real4),trac_varID(i))
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_movie
>>>>>>>> call
>>>>>>>> write_data(real(tpasnew(:,:,:,itpas_movie(i)),real4),tpas_varID(i))
>>>>>>>> enddo
>>>>>>>> endif
>>>>>>>>
>>>>>>>> ! === close netCDF dataset: exit data mode & save
>>>>>>>> -----------------!
>>>>>>>> OP = 'close netCDF file'
>>>>>>>> call handle_err(nfmpi_close(NCID))
>>>>>>>>
>>>>>>>> if (myrank==0) then
>>>>>>>> print *, ''
>>>>>>>> print *, 'wrote netCDF output to file ', outfile
>>>>>>>> print *, ''
>>>>>>>> end if
>>>>>>>>
>>>>>>>> contains
>>>>>>>> !----------------------------------------------------------------
>>>>>>>> subroutine write_data(array_in,varID)
>>>>>>>> ! IF: array_in = data array (3D) to be written to netCDF
>>>>>>>> ! varID = data variable netCDF ID
>>>>>>>>
>>>>>>>> use atham_module, only : iflgs ! topography mask: topography
>>>>>>>> (iflgs=0) remains undef
>>>>>>>>
>>>>>>>> use atham_module, only : myrank
>>>>>>>>
>>>>>>>> real(real4),dimension(nx,ny,nz), intent (in ) :: array_in
>>>>>>>> integer(kint), intent(in ) :: varID
>>>>>>>>
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>> ! local
>>>>>>>> variables !
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>>
>>>>>>>> real(real4), dimension(nx,ny,nz) :: array_masked
>>>>>>>>
>>>>>>>> ! dimensions & bounds control
>>>>>>>> integer(MPI_OFFSET_KIND) :: dimlenx,
>>>>>>>> dimleny, dimlenz, dimlent
>>>>>>>> integer(kint), dimension(NF_max_var_dims) :: dimIDs
>>>>>>>>
>>>>>>>> ! dimensions & bounds control
>>>>>>>> if (myrank == 0) then
>>>>>>>> if (bounds_control) then
>>>>>>>> print *, ''
>>>>>>>> print *, 'atham var in xyz ',shape(array_in)
>>>>>>>> print *, 'atham var reshaped to xyzt
>>>>>>>> ',shape(reshape(array_in(2:nx-1,y_st:y_sp,1:nz),shape=(/nx-2,y_cl,nz,1/)))
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_inq_vardimid(NCID, u_varID,
>>>>>>>> dimIDs)) ! taken unew as ref
>>>>>>>> call handle_err(nfmpi_inq_dimlen(NCID, dimIDs(1),
>>>>>>>> dimlenx))
>>>>>>>> call handle_err(nfmpi_inq_dimlen(NCID, dimIDs(2),
>>>>>>>> dimleny))
>>>>>>>> call handle_err(nfmpi_inq_dimlen(NCID, dimIDs(3),
>>>>>>>> dimlenz))
>>>>>>>> call handle_err(nfmpi_inq_dimlen(NCID, dimIDs(4),
>>>>>>>> dimlent))
>>>>>>>> print *, 'netCDF var dimensions in xyzt
>>>>>>>> ',dimlenx,dimleny,dimlenz,dimlent
>>>>>>>> end if
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! mask topography
>>>>>>>> where (iflgs == 1)
>>>>>>>> array_masked = array_in
>>>>>>>> elsewhere
>>>>>>>> array_masked = undef
>>>>>>>> endwhere
>>>>>>>>
>>>>>>>> ! mask control
>>>>>>>> !print *, ''
>>>>>>>> !print *, 'topography mask'
>>>>>>>> !print *, iflgs(1:5,3,3)
>>>>>>>> !print *, iflgs(1:5,3,2)
>>>>>>>> !print *, iflgs(1:5,3,1)
>>>>>>>> !print *, ''
>>>>>>>> !print *, 'masked variable array'
>>>>>>>> !print *, array_masked(:,3,3)
>>>>>>>> !print *, array_masked(:,3,2)
>>>>>>>> !print *, array_masked(:,3,1)
>>>>>>>> !print *, ''
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_begin_indep_data(NCID))
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID,
>>>>>>>> varID,(/istart,jstart,1_MPI_OFFSET_KIND,tcounter/), &
>>>>>>>> (/nx_mpi,y_cl_mpi,nz_mpi,1_MPI_OFFSET_KIND/),
>>>>>>>> reshape(array_masked(2:nx-1,y_st:y_sp,1:nz),shape=(/nx-2,y_cl,nz,1/))
>>>>>>>>
>>>>>>>> ),OP)
>>>>>>>> call handle_err(nfmpi_end_indep_data(NCID))
>>>>>>>>
>>>>>>>> end subroutine write_data
>>>>>>>> !
>>>>>>>> ----------------------------------------------------------------
>>>>>>>> end subroutine netcdf_data
>>>>>>>> !=====================================================================
>>>>>>>>
>>>>>>>>
>>>>>>>> subroutine netcdf_define(outfile,out_type,timetot,ncid_def)
>>>>>>>> ! IF: outfile = picture_file or movie_file (path/filename for
>>>>>>>> output created in atham_module)
>>>>>>>> ! out_type = 1 (picture_file) or 2 (movie_file)
>>>>>>>> ! timetot = simulation time since simulation start (from
>>>>>>>> atham/atham_out)
>>>>>>>> ! ncid_def = netCDF ID assigned to file during creation
>>>>>>>>
>>>>>>>> use atham_module, only : var_tgas_picture,
>>>>>>>> var_tgas_movie, &
>>>>>>>> des_tgas_picture,
>>>>>>>> des_tgas_movie, &
>>>>>>>> var_trac_picture,
>>>>>>>> var_trac_movie, &
>>>>>>>> des_trac_picture,
>>>>>>>> des_trac_movie, &
>>>>>>>> var_tpas_picture,
>>>>>>>> var_tpas_movie, &
>>>>>>>> des_tpas_picture, des_tpas_movie
>>>>>>>>
>>>>>>>> use atham_module, only : iindex, spinup, nrep, periodt
>>>>>>>>
>>>>>>>> use atham_module, only : volcano_setup,
>>>>>>>> coignimbrite_setup, &
>>>>>>>> nuclear_setup, convection_setup,
>>>>>>>> procsconfig_setup
>>>>>>>>
>>>>>>>> use atham_module, only : x, y, z, xg, yg
>>>>>>>>
>>>>>>>> use atham_module, only : grid2geo, sim_d_YYYY, sim_d_MM,
>>>>>>>> sim_d_DD, &
>>>>>>>> sim_t_hh, sim_t_mm, sim_t_ss
>>>>>>>>
>>>>>>>> use atham_module, only : my_cart, myrank
>>>>>>>>
>>>>>>>> character(*), intent(in ) :: outfile
>>>>>>>> integer(kint), intent(in ) :: out_type
>>>>>>>> real(kreal), intent(in ) :: timetot
>>>>>>>> integer(kint), intent( out) :: ncid_def
>>>>>>>>
>>>>>>>> !------------------------------------------------------------------!
>>>>>>>>
>>>>>>>> ! local
>>>>>>>> variables !
>>>>>>>> !------------------------------------------------------------------!
>>>>>>>>
>>>>>>>> integer(MPI_OFFSET_KIND) :: ntx_mpi,nty_mpi,nz_mpi
>>>>>>>> character(len=50) :: OP ! debugging variable
>>>>>>>> corresponding to the called netCDF OPeration
>>>>>>>> integer(kint) :: NCID
>>>>>>>> integer(kint) :: i
>>>>>>>> character(len=50) :: var, des
>>>>>>>> real(kreal) :: toffset
>>>>>>>> real(real4) :: coord_scale,
>>>>>>>> coord_x_offset, coord_y_offset, coord_t_offset
>>>>>>>> logical :: DynMovie_flag
>>>>>>>> character(len=3) :: coord_x_name,
>>>>>>>> coord_y_name, coord_z_name
>>>>>>>> character(len=4) :: coord_t_name
>>>>>>>> character(len=9) :: coord_x_lname,
>>>>>>>> coord_y_lname
>>>>>>>> character(len=13) :: coord_x_units,
>>>>>>>> coord_y_units
>>>>>>>> character(len=33) :: coord_t_units
>>>>>>>> character(len=8) :: date
>>>>>>>> character(len=10) :: time
>>>>>>>> character(len=19) :: simdatetime
>>>>>>>> character(len=10) :: simdate
>>>>>>>> character(len=20) :: login
>>>>>>>> character(len=63) :: hist_string
>>>>>>>> integer(kint) :: isetup
>>>>>>>> character(len=12) :: model_setup(6)
>>>>>>>> data model_setup /'volcano ', 'coignimbrite', 'nuclear ',
>>>>>>>> 'convection ', 'procsconfig ','unknown '/
>>>>>>>>
>>>>>>>> ! replace DynMovie precompiler flag by program flag for more
>>>>>>>> convenient code structure
>>>>>>>> DynMovie_flag = .false.
>>>>>>>> #ifdef DynMovie
>>>>>>>> DynMovie_flag = .true.
>>>>>>>> #endif
>>>>>>>>
>>>>>>>> ! re-assign dimension lengths as an MPI_Offset type (necessary for
>>>>>>>> the PnetCDF interface)
>>>>>>>> ! cut-off domain borders
>>>>>>>> ntx_mpi = ntx-2
>>>>>>>> nty_mpi = nty-2
>>>>>>>> nz_mpi = nz
>>>>>>>>
>>>>>>>> ! convert date and time integers into string 'YYYYMMDD hh:mm:ss'
>>>>>>>> write (simdatetime,10)
>>>>>>>> sim_d_YYYY,sim_d_MM,sim_d_DD,sim_t_hh,sim_t_mm,sim_t_ss
>>>>>>>> 10 format (I4,'-',I2,'-',I2,1X,I2,':',I2,':',I2)
>>>>>>>> simdate = simdatetime(1:10)
>>>>>>>>
>>>>>>>> ! switch between cartesian axes/time (default) and geographic
>>>>>>>> axes/date (linear approx)
>>>>>>>> coord_scale = 1._real4
>>>>>>>> coord_x_offset = 0._real4
>>>>>>>> coord_x_name = 'x '
>>>>>>>> coord_x_lname = 'x_coor '
>>>>>>>> coord_x_units = 'm '
>>>>>>>> coord_y_offset = 0._real4
>>>>>>>> coord_y_name = 'y '
>>>>>>>> coord_y_lname = 'y_coor '
>>>>>>>> coord_y_units = 'm '
>>>>>>>> coord_z_name = 'z '
>>>>>>>> coord_t_offset = 0._real4
>>>>>>>> coord_t_name = 'time'
>>>>>>>> coord_t_units = 'sec '
>>>>>>>> if (grid2geo) then
>>>>>>>> call
>>>>>>>> gridconvert(coord_scale,coord_x_offset,coord_x_name,coord_x_lname,coord_x_units,
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> &
>>>>>>>>
>>>>>>>> coord_y_offset,coord_y_name,coord_y_lname,coord_y_units, &
>>>>>>>>
>>>>>>>> coord_z_name,coord_t_offset,coord_t_name,coord_t_units)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! === create netCDF dataset: enter define mode
>>>>>>>> --------------------!
>>>>>>>>
>>>>>>>> OP = 'enter define mode'
>>>>>>>> call handle_err(nfmpi_create(my_cart, trim(outfile), 0,
>>>>>>>> MPI_INFO_NULL,NCID),OP)
>>>>>>>> ncid_def = NCID
>>>>>>>>
>>>>>>>> ! === create netCDF dataset: assign metadata
>>>>>>>> ----------------------!
>>>>>>>>
>>>>>>>> ! --- set global attributes
>>>>>>>> ----------------!
>>>>>>>> OP = 'set global attributes'
>>>>>>>> if (out_type == 1) then
>>>>>>>> toffset = timetot - spinup
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL, 'title',
>>>>>>>> 13_MPI_OFFSET_KIND, 'atham_picture'),OP)
>>>>>>>> else
>>>>>>>> toffset = timetot - spinup - nrep*periodt
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL, 'title',
>>>>>>>> 11_MPI_OFFSET_KIND, 'atham_movie'),OP)
>>>>>>>> end if
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL, 'content',
>>>>>>>> 223_MPI_OFFSET_KIND, &
>>>>>>>> 'This data corresponds to a simulation run of the
>>>>>>>> mesoscale-beta/gamma type non-hydrostatic Active Tracer
>>>>>>>> High-resolution Atmospheric Model ATHAM ' &
>>>>>>>> // '(Oberhuber et al., 1998, Herzog, 1998), for a setup as
>>>>>>>> specified under history'),OP) ! ######### add further
>>>>>>>> description as
>>>>>>>> necessary
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL, 'version',
>>>>>>>> 21_MPI_OFFSET_KIND, 'no current versioning'),OP) ! ######### change
>>>>>>>> once available
>>>>>>>>
>>>>>>>> call date_and_time(date, time)
>>>>>>>> call getlog(login)
>>>>>>>> if (volcano_setup) then
>>>>>>>> isetup=1
>>>>>>>> else if (coignimbrite_setup) then
>>>>>>>> isetup=2
>>>>>>>> else if (nuclear_setup) then
>>>>>>>> isetup=3
>>>>>>>> else if (convection_setup) then
>>>>>>>> isetup=4
>>>>>>>> else if (procsconfig_setup) then
>>>>>>>> isetup=5
>>>>>>>> else
>>>>>>>> isetup=6
>>>>>>>> end if
>>>>>>>> hist_string = trim(date//' '//time//' '//trim(login)//' atham
>>>>>>>> '//model_setup(isetup))
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL, 'history',
>>>>>>>> 63_MPI_OFFSET_KIND, hist_string),OP) ! 1 line/modification: - cmd
>>>>>>>> arg
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL,
>>>>>>>> 'Conventions', 28_MPI_OFFSET_KIND, 'COARDS,
>>>>>>>> //cf-pcmdi.llnl.gov/'),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, NF_GLOBAL,
>>>>>>>> 'sim_onset_date_time', 19_MPI_OFFSET_KIND, simdatetime),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, NF_GLOBAL,
>>>>>>>> 'sim_time_offset', NF_FLOAT, 1_MPI_OFFSET_KIND,
>>>>>>>> real(toffset,real4)),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, NF_GLOBAL,
>>>>>>>> 'x_zoom_location', NF_FLOAT, 1_MPI_OFFSET_KIND,
>>>>>>>> real(x(max(1,iindex)),real4)),OP)
>>>>>>>>
>>>>>>>> ! --- define dimensions
>>>>>>>> --------------------!
>>>>>>>> OP = 'define dimensions'
>>>>>>>> call handle_err(nfmpi_def_dim(NCID, trim(coord_x_name),
>>>>>>>> ntx_mpi, x_dimID),OP)
>>>>>>>> if (nty>4) then
>>>>>>>> call handle_err(nfmpi_def_dim(NCID, trim(coord_y_name),
>>>>>>>> nty_mpi, y_dimID),OP)
>>>>>>>> else
>>>>>>>> call handle_err(nfmpi_def_dim(NCID, trim(coord_y_name),
>>>>>>>> 1_MPI_OFFSET_KIND, y_dimID),OP)
>>>>>>>> end if
>>>>>>>> call handle_err(nfmpi_def_dim(NCID, trim(coord_z_name),
>>>>>>>> nz_mpi, z_dimID),OP)
>>>>>>>> call handle_err(nfmpi_def_dim(NCID, coord_t_name,
>>>>>>>> nfmpi_unlimited, t_dimID),OP)
>>>>>>>>
>>>>>>>> ! --- define variables
>>>>>>>> ---------------------!
>>>>>>>> ! dimensions
>>>>>>>> OP = 'define variables - dimensions'
>>>>>>>> call handle_err(nfmpi_def_var(NCID, trim(coord_x_name),
>>>>>>>> NF_FLOAT, 1, x_dimID, x_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, trim(coord_y_name),
>>>>>>>> NF_FLOAT, 1, y_dimID, y_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, trim(coord_z_name),
>>>>>>>> NF_FLOAT, 1, z_dimID, z_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, coord_t_name,
>>>>>>>> NF_FLOAT, 1, t_dimID, t_varID),OP)
>>>>>>>>
>>>>>>>> ! dynamics data
>>>>>>>> OP = 'define variables - data'
>>>>>>>> if (out_type == 1 .OR. DynMovie_flag == .true.) then
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'u', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), u_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'v', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), v_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'w', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), w_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'pnew', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), pnew_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'tetnew', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), tetnew_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'tempnew', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), tempnew_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'density', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), density_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'turbhor', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), turbhor_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'turbver', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), turbver_varID),OP)
>>>>>>>> call handle_err(nfmpi_def_var(NCID, 'turblen', NF_FLOAT, 4,
>>>>>>>> (/x_dimID,y_dimID,z_dimID,t_dimID/), turblen_varID),OP)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! tracer data
>>>>>>>> if (out_type == 1) then
>>>>>>>> do i=1,ntgas_picture
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_tgas_picture(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> tgas_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_trac_picture(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> trac_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_picture
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_tpas_picture(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> tpas_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> var='d_'//trim(var_trac_picture(i))
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var,
>>>>>>>> NF_FLOAT, 3, (/x_dimID,y_dimID,t_dimID/), d_trac_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> else if (out_type == 2) then
>>>>>>>> do i=1,ntgas_movie
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_tgas_movie(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> tgas_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_movie
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_trac_movie(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> trac_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_movie
>>>>>>>> call handle_err(nfmpi_def_var(NCID, var_tpas_movie(i),
>>>>>>>> NF_FLOAT, 4, (/x_dimID,y_dimID,z_dimID,t_dimID/),
>>>>>>>> tpas_varID(i)),OP)
>>>>>>>> enddo
>>>>>>>> endif
>>>>>>>>
>>>>>>>> ! --- assign attributes
>>>>>>>> --------------------!
>>>>>>>> ! common attributes: long_name, units, valid_range, scale_factor,
>>>>>>>> add_offset, _FillValue, history, conventions
>>>>>>>> ! dimensions
>>>>>>>>
>>>>>>>> OP = 'assign attributes - dimensions'
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, x_varID,
>>>>>>>> 'long_name', 9_MPI_OFFSET_KIND,
>>>>>>>> trim(coord_x_lname)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, x_varID,
>>>>>>>> 'units', 13_MPI_OFFSET_KIND,
>>>>>>>> trim(coord_x_units)),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, x_varID,
>>>>>>>> 'scale_factor', NF_FLOAT, 1_MPI_OFFSET_KIND, coord_scale),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, x_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, coord_x_offset),OP)
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, y_varID,
>>>>>>>> 'long_name', 9_MPI_OFFSET_KIND,
>>>>>>>> trim(coord_y_lname)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, y_varID,
>>>>>>>> 'units', 13_MPI_OFFSET_KIND,
>>>>>>>> trim(coord_y_units)),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, y_varID,
>>>>>>>> 'scale_factor', NF_FLOAT, 1_MPI_OFFSET_KIND, coord_scale),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, y_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, coord_y_offset),OP)
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, z_varID,
>>>>>>>> 'long_name', 6_MPI_OFFSET_KIND, 'height'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, z_varID,
>>>>>>>> 'units', 1_MPI_OFFSET_KIND, 'm'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, z_varID,
>>>>>>>> 'positive', 2_MPI_OFFSET_KIND, 'up'),OP)
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, t_varID,
>>>>>>>> 'long_name', 4_MPI_OFFSET_KIND, coord_t_name),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, t_varID,
>>>>>>>> 'units', 33_MPI_OFFSET_KIND,
>>>>>>>> trim(coord_t_units)),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, t_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, coord_t_offset),OP)
>>>>>>>>
>>>>>>>> ! OP = 'exit define mode NOW for testing
>>>>>>>> ****************************************************' ! ah_testing
>>>>>>>> ********************************
>>>>>>>> ! call
>>>>>>>> handle_err(nfmpi_enddef(NCID),OP)
>>>>>>>> ! ah_testing ********************************
>>>>>>>>
>>>>>>>> ! dynamics data
>>>>>>>> OP = 'assign attributes - data'
>>>>>>>> if (out_type == 1 .OR. DynMovie_flag == .true.) then
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, u_varID,
>>>>>>>> 'long_name', 6_MPI_OFFSET_KIND, 'x_wind'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, u_varID,
>>>>>>>> 'units', 7_MPI_OFFSET_KIND, 'm sec-1'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, u_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, u_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, u_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, v_varID,
>>>>>>>> 'long_name', 6_MPI_OFFSET_KIND, 'y_wind'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, v_varID,
>>>>>>>> 'units', 7_MPI_OFFSET_KIND, 'm sec-1'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, v_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, v_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, v_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, w_varID,
>>>>>>>> 'long_name', 19_MPI_OFFSET_KIND,
>>>>>>>> 'upward_air_velocity'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, w_varID,
>>>>>>>> 'units', 7_MPI_OFFSET_KIND, 'm sec-1'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, w_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, w_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, w_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, pnew_varID,
>>>>>>>> 'long_name', 20_MPI_OFFSET_KIND,
>>>>>>>> 'air_pressure_anomaly'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, pnew_varID,
>>>>>>>> 'units', 2_MPI_OFFSET_KIND, 'Pa'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, pnew_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, pnew_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, pnew_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, tetnew_varID,
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND,
>>>>>>>> 'air_potential_temperature'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, tetnew_varID,
>>>>>>>> 'units', 1_MPI_OFFSET_KIND, 'K'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tetnew_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tetnew_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tetnew_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, tempnew_varID,
>>>>>>>> 'long_name', 15_MPI_OFFSET_KIND, 'air_temperature'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, tempnew_varID,
>>>>>>>> 'units', 7_MPI_OFFSET_KIND, 'Celsius'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tempnew_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tempnew_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, -273.15),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, tempnew_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, density_varID,
>>>>>>>> 'long_name', 11_MPI_OFFSET_KIND, 'air_density'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, density_varID,
>>>>>>>> 'units', 6_MPI_OFFSET_KIND, 'kg m-3'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, density_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, density_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, density_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turbhor_varID,
>>>>>>>> 'long_name', 15_MPI_OFFSET_KIND, 'hor_turb_energy'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turbhor_varID,
>>>>>>>> 'units', 8_MPI_OFFSET_KIND, 'm2 sec-2'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbhor_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbhor_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbhor_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turbver_varID,
>>>>>>>> 'long_name', 15_MPI_OFFSET_KIND, 'ver_turb_energy'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turbver_varID,
>>>>>>>> 'units', 8_MPI_OFFSET_KIND, 'm2 sec-2'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbver_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbver_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turbver_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>>
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turblen_varID,
>>>>>>>> 'long_name', 17_MPI_OFFSET_KIND,
>>>>>>>> 'turb_length_scale'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_text(NCID, turblen_varID,
>>>>>>>> 'units', 1_MPI_OFFSET_KIND, 'm'),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turblen_varID,
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1.),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turblen_varID,
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> ! call handle_err(nfmpi_put_att_real(NCID, turblen_varID,
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! tracer data
>>>>>>>> if (out_type == 1) then
>>>>>>>> print *, 'NetCDF output type (1-picture,2-movie) ',out_type !
>>>>>>>> ##### DO NOT REMOVE!!! optimization problem, trace source
>>>>>>>> do i=1,ntgas_picture
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tgas_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND,
>>>>>>>> des_tgas_picture(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tgas_varID(i),
>>>>>>>> 'units', 9_MPI_OFFSET_KIND, 'gram kg-1'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1000.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, trac_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND,
>>>>>>>> des_trac_picture(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, trac_varID(i),
>>>>>>>> 'units', 9_MPI_OFFSET_KIND, 'gram kg-1'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1000.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_picture
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tpas_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND,
>>>>>>>> des_tpas_picture(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tpas_varID(i),
>>>>>>>> 'units', 1_MPI_OFFSET_KIND, '%'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 100.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_picture
>>>>>>>> des='deposited_'//trim(des_trac_picture(i))
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, d_trac_varID(i),
>>>>>>>> 'long_name', 35_MPI_OFFSET_KIND, des),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, d_trac_varID(i),
>>>>>>>> 'units', 9_MPI_OFFSET_KIND, 'gram m-2'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, d_trac_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1000.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, d_trac_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, d_trac_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> else if (out_type == 2) then
>>>>>>>> print *, 'NetCDF output type (1-picture,2-movie) ',out_type !
>>>>>>>> ##### DO NOT REMOVE!!! optimization problem, trace source
>>>>>>>> do i=1,ntgas_movie
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tgas_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND, des_tgas_movie(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tgas_varID(i),
>>>>>>>> 'units', 9_MPI_OFFSET_KIND, 'gram kg-1'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1000.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tgas_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntrac_movie
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, trac_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND, des_trac_movie(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, trac_varID(i),
>>>>>>>> 'units', 9_MPI_OFFSET_KIND, 'gram kg-1'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 1000.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, trac_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> do i=1,ntpas_movie
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tpas_varID(i),
>>>>>>>> 'long_name', 25_MPI_OFFSET_KIND, des_tpas_movie(i)),OP)
>>>>>>>> call handle_err(nfmpi_put_att_text(NCID, tpas_varID(i),
>>>>>>>> 'units', 1_MPI_OFFSET_KIND, '%'),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> 'scale_factor',NF_FLOAT, 1_MPI_OFFSET_KIND, 100.),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> 'add_offset', NF_FLOAT, 1_MPI_OFFSET_KIND, 0.0),OP)
>>>>>>>> call handle_err(nfmpi_put_att_real(NCID, tpas_varID(i),
>>>>>>>> '_FillValue', NF_FLOAT, 1_MPI_OFFSET_KIND, undef),OP)
>>>>>>>> enddo
>>>>>>>> end if
>>>>>>>>
>>>>>>>> ! === create netCDF dataset: exit define mode, enter data mode
>>>>>>>> ----!
>>>>>>>> OP = 'exit define mode'
>>>>>>>> call handle_err(nfmpi_enddef(NCID),OP)
>>>>>>>>
>>>>>>>> ! --- write data values
>>>>>>>> --------------------!
>>>>>>>> ! dimensions (only one process writes the whole dimension, i.e.
>>>>>>>> accross all processes)
>>>>>>>> OP = 'write data values - dimensions'
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_begin_indep_data(NCID))
>>>>>>>>
>>>>>>>> if (myrank==0) then
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID, x_varID,
>>>>>>>> (/1_MPI_OFFSET_KIND/), (/ntx_mpi/),real(xg(2:ntx-1),real4)),OP)
>>>>>>>> if (nty>4) then
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID, y_varID,
>>>>>>>> (/1_MPI_OFFSET_KIND/), (/nty_mpi/), real(yg(2:nty-1),real4)),OP)
>>>>>>>> else
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID, y_varID,
>>>>>>>> (/1_MPI_OFFSET_KIND/), (/1_MPI_OFFSET_KIND/),
>>>>>>>> real(yg(nty-1),real4)),OP)
>>>>>>>> end if
>>>>>>>> call handle_err(nfmpi_put_vara_real(NCID, z_varID,
>>>>>>>> (/1_MPI_OFFSET_KIND/), (/nz_mpi/), real(z,real4)),OP)
>>>>>>>> end if
>>>>>>>>
>>>>>>>> call handle_err(nfmpi_end_indep_data(NCID))
>>>>>>>>
>>>>>>>> contains
>>>>>>>> !
>>>>>>>> ----------------------------------------------------------------
>>>>>>>> subroutine
>>>>>>>> gridconvert(cd_scale,cd_x_offset,cd_x_name,cd_x_lname,cd_x_units, &
>>>>>>>>
>>>>>>>> cd_y_offset,cd_y_name,cd_y_lname,cd_y_units, &
>>>>>>>>
>>>>>>>> cd_z_name,cd_t_offset,cd_t_name,cd_t_units)
>>>>>>>> ! IF: cd_* = coordinate parameters
>>>>>>>>
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>> ! linear conversion from cartesian grid into lon/lat
>>>>>>>> system !
>>>>>>>> ! only valid as first approxiamtion over small
>>>>>>>> domains !
>>>>>>>> ! to be replaced once a dyn core in a lat/lon system
>>>>>>>> available !
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>>
>>>>>>>> use atham_module, only : deglat_start, deglon_start
>>>>>>>>
>>>>>>>> use phys_constants, only : pi
>>>>>>>>
>>>>>>>> real(real4), intent(inout) :: cd_scale, cd_x_offset,
>>>>>>>> cd_y_offset, cd_t_offset
>>>>>>>> character(len=3), intent(inout) :: cd_x_name, cd_y_name,
>>>>>>>> cd_z_name
>>>>>>>> character(len=4), intent(inout) :: cd_t_name
>>>>>>>> character(len=9), intent(inout) :: cd_x_lname, cd_y_lname
>>>>>>>> character(len=13), intent(inout) :: cd_x_units, cd_y_units
>>>>>>>> character(len=33), intent(inout) :: cd_t_units
>>>>>>>>
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>> ! local
>>>>>>>> variables !
>>>>>>>>
>>>>>>>> !--------------------------------------------------------------!
>>>>>>>>
>>>>>>>> real(kreal) :: r
>>>>>>>> real(kreal), parameter :: a = 6378137._kreal !
>>>>>>>> (WGS 84)
>>>>>>>> real(kreal), parameter :: b = 6356752.3142_kreal !
>>>>>>>> (WGS 84)
>>>>>>>>
>>>>>>>> ! ellipsoid radius at ATHAM domain start
>>>>>>>> r = a*b / sqrt((b*b*(cos(deglat_start*pi/180._kreal))**2
>>>>>>>> )+(a*a*(sin(deglat_start*pi/180._kreal))**2 ))
>>>>>>>>
>>>>>>>> cd_scale = real((1/r) * (180._kreal/pi),real4)
>>>>>>>>
>>>>>>>> cd_x_offset = deglon_start
>>>>>>>> cd_x_name = 'lon'
>>>>>>>> cd_x_lname = 'longitude'
>>>>>>>> cd_x_units = 'degrees_east '
>>>>>>>>
>>>>>>>> cd_y_offset = deglat_start
>>>>>>>> cd_y_name = 'lat'
>>>>>>>> cd_y_lname = 'latitude '
>>>>>>>> cd_y_units = 'degrees_north'
>>>>>>>>
>>>>>>>> cd_z_name = 'lev'
>>>>>>>>
>>>>>>>> coord_t_offset = real((sim_t_hh*3600 + sim_t_mm*60
>>>>>>>> +sim_t_ss),real4)
>>>>>>>> coord_t_name = 'date'
>>>>>>>> coord_t_units = 'seconds since '//simdate//' 00:00:00 [local]'
>>>>>>>> ! COARDS conventions require time zone from UTC to be specified
>>>>>>>> ! not necessary for local area version of ATHAM
>>>>>>>>
>>>>>>>> end subroutine gridconvert
>>>>>>>> !
>>>>>>>> ----------------------------------------------------------------
>>>>>>>>
>>>>>>>> end subroutine netcdf_define
>>>>>>>> !=====================================================================
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> !
>>>>>>>> ======================================================================
>>>>>>>>
>>>>>>>>
>>>>>>>> subroutine handle_err(status,op)
>>>>>>>>
>>>>>>>> use atham_module, only : myrank
>>>>>>>>
>>>>>>>> integer, intent (in ) :: status
>>>>>>>> character(len=*), optional, intent (in ) :: op
>>>>>>>>
>>>>>>>> if ( status /= nf_noerr ) then
>>>>>>>> if (myrank==0) then
>>>>>>>> print *, ""
>>>>>>>> print *, "P-NetCDF operation: ",trim(op)
>>>>>>>> print *, trim(nfmpi_strerror(status))
>>>>>>>> print *, "atham_pnetcdf error: program execution stopped"
>>>>>>>> print *, ""
>>>>>>>>
>>>>>>>> stop
>>>>>>>> end if
>>>>>>>> end if
>>>>>>>>
>>>>>>>> return
>>>>>>>>
>>>>>>>> end subroutine handle_err
>>>>>>>> !
>>>>>>>> ==============================================================================
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> end module atham_pnetcdf
>>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> ______________________________________________________________________
>>>>
>>>> Alex HOFFMANN PhD Candidate
>>>> Centre for Atmospheric Science
>>>> Department of Geography University of Cambridge
>>>> Downing Place, Cambridge CB2 3EN, UK
>>>> e-mail: ah519 at cam.ac.uk tel: +44 (0)1223 766581
>>>> www.geog.cam.ac.uk/people/hoffmann/
>>>> ______________________________________________________________________
>>>>
>>>
>>
>> --
>> ______________________________________________________________________
>>
>> Alex HOFFMANN PhD Candidate
>> Centre for Atmospheric Science
>> Department of Geography University of Cambridge
>> Downing Place, Cambridge CB2 3EN, UK
>> e-mail: ah519 at cam.ac.uk tel: +44 (0)1223 766581
>> www.geog.cam.ac.uk/people/hoffmann/
>> ______________________________________________________________________
>> !-*- F90 -*-
>> program mpi_setup_netCDFout
>>
>> !--------------------------------------------------------------------!
>> ! author: Alex HOFFMANN !
>> ! email: ah519 at cam.ac.uk !
>> ! date: February 2009 !
>> ! version: development !
>> ! !
>> ! simple program to test mpi runs & netCDF output (serial & parallel)!
>> !--------------------------------------------------------------------!
>>
>> ! test code designed to be used with the variables from the mesoscale
>> ! atmospheric model ATHAM
>> ! re-structured as stand-alone
>> ! reads and writes as netCDF an index matrix for a 4x4 matrix in real4
>>
>> ! use precision
>> ! use atham_module, only : nx,ny,ntx,nty,npx,npy
>>
>> implicit none
>>
>> !private
>> !public
>>
>> ! --- mpi interface for Fortran 77 only, use header files-----------
>> include 'mpif.h'
>>
>>
>> ! --- precision ----------------------------------------------------
>> integer, parameter :: real4 =
>> selected_real_kind(6), &
>> real8 =
>> selected_real_kind(12), &
>> kreal = real8
>> integer, parameter :: int4 =
>> selected_int_kind(7), &
>> int8 =
>> selected_int_kind(13), &
>> kint = int4
>> ! --- atham_module -------------------------------------------------
>> integer(kint) :: nx,ny,ntx,nty,npx,npy
>> integer(kint) :: my_ranks,
>> nprocs, &
>> my_cart, myrank, mycoor(2)
>> logical :: mpi_initialized, mpi
>> ! --- other --------------------------------------------------------
>> integer(kint) :: ierror
>> real(real4), allocatable, dimension(:,:) :: array_g, array
>>
>>
>>
>> mpi_initialized=.false.
>> ! ---for initial testing--------------------------------------------
>> nx = 4_kint ! user-specified number of grid points along x
>> ny = 4_kint
>> npx = 2_kint ! number of processes along x
>> npy = 2_kint
>>
>> ntx = nx ! re-adjusted total number of grid points along x
>> nty = ny
>> nx = nx/npx ! number of grid points along x per process
>> ny = ny/npy
>>
>>
>> ! --- MAIN
>> ! --- mpi setup calls ----------------------------------------------
>> call initialize_mpi
>> call setup_mpi
>>
>> ! --- netCDF out calls ---------------------------------------------
>> call nc_out_s
>> call nc_out_p
>>
>> ! --- mpi finalize call --------------------------------------------
>> call mpi_Finalize(ierror)
>> ! --- END MAIN
>>
>> contains
>> !=====================================================================
>> subroutine initialize_mpi()
>>
>> ! use atham_module, only : mpi, mpi_initialized
>> ! use atham_module, only : myrank, nprocs
>>
>> integer(kint) :: ierror, errorclass
>> character(len=100) :: message
>>
>> ! assign variables for serial version
>> ierror=0
>> myrank=0
>> nprocs=1
>> mpi=.false.
>> ! initialize MPI
>> mpi=.true.
>> if (.not. mpi_initialized) call mpi_init(ierror)
>> call MPI_Comm_rank(mpi_comm_world,myrank,ierror)
>> call MPI_Comm_size(mpi_comm_world,nprocs,ierror)
>> mpi=.true.
>> mpi_initialized=.true.
>>
>> if (myrank==0) then
>> print *, 'Myrank 0: ',myrank
>> print *, 'MPI initialized!'
>> print *, 'Number of processes [nprocs]: ',nprocs
>> endif
>>
>> if (ierror>0) then
>> call MPI_Error_string(ierror,message,100)
>> call MPI_Error_class(ierror,errorclass)
>> call exec_stop(message,errorclass)
>> endif
>>
>> end subroutine initialize_mpi
>> !=====================================================================
>> subroutine setup_mpi()
>>
>> ! use atham_module, only : mpi, mpi_initialized
>> ! use atham_module, only : myrank,nprocs,my_cart,mycoor
>>
>> integer(kint) :: ierror, errorclass
>> character(len=100) :: message
>> integer(kint) :: idims(2),shape_array(2)
>> logical :: reorder
>> integer(kint) :: i,j,istart,jstart
>>
>> if (nprocs/=npx*npy) then
>> print *, 'nprocs /= npx*npy - execution aborted'
>> call mpi_abort(mpi_comm_world,ierror)
>> stop
>> endif
>> ! create virtual topology
>> idims(1:2) = (/npx,npy/)
>> reorder = .true.
>>
>> call
>> MPI_Cart_create(mpi_comm_world,2,idims,(/.false.,.false./),reorder,my_cart,ierror)
>>
>> call MPI_Comm_rank(my_cart,myrank,ierror) ! ranks on new communicator
>> call MPI_Cart_coords(my_cart,myrank,2,mycoor,ierror)
>>
>> if(myrank==0) then
>> print *, 'Coords of process ',myrank,' :',mycoor
>> endif
>>
>> if (ierror>0) then
>> call MPI_Error_string(ierror,message,100)
>> call MPI_Error_class(ierror,errorclass)
>> call exec_stop(message,errorclass)
>> endif
>>
>> ! allocate and assign global grid on each process
>> allocate(array_g(ntx,nty))
>>
>> array_g(1,:) = (/ 11.0,12.0,13.0,14.0 /)
>> array_g(2,:) = (/ 21.0,22.0,23.0,24.0 /)
>> array_g(3,:) = (/ 31.0,32.0,33.0,34.0 /)
>> array_g(4,:) = (/ 41.0,42.0,43.0,44.0 /)
>> if (myrank==0) then
>> print *, 'whole 4*4 array on each process:'
>> print *, array_g(1,:)
>> print *, array_g(2,:)
>> print *, array_g(3,:)
>> print *, array_g(4,:)
>> endif
>>
>> ! allocate and assign subgrid portions on various processes
>> allocate(array(nx,ny))
>> istart = mycoor(1)*nx
>> jstart = mycoor(2)*ny
>>
>> do i=1,nx
>> do j=1,ny
>> array(i,j) = array_g(istart+i,jstart+j)
>> enddo
>> enddo
>>
>> ! output results from one process
>> if (myrank==0) then
>> print *, ''
>> print *, 'subgrid array on process with coords ',mycoor,':'
>> shape_array = shape(array)
>> print *, 'shape: ',shape_array
>> do i=1,shape_array(1)
>> print *, array(i,:)
>> enddo
>> print *, ''
>> endif
>>
>> end subroutine setup_mpi
>> !=====================================================================
>> subroutine nc_out_s()
>>
>> use netcdf
>>
>> ! use atham_module, only : makename
>> ! use atham_module, only : mycoor
>> ! use atham_module, only : nx,ny
>>
>> integer(kint) :: NCID
>> character(len=10) :: name
>> integer(kint) :: x_dimID,y_dimID,x_varID,y_varID,array_varID
>>
>>
>> call makename('out',mycoor,name)
>> name = trim(name)//'.nc'
>>
>> call handle_err_nc_s(nf90_create(path = './'//trim(name), cmode =
>> 0, ncid = NCID))
>> call handle_err_nc_s(nf90_def_dim(ncid = NCID, name = 'x', len
>> = nx, dimid = x_dimID))
>> call handle_err_nc_s(nf90_def_dim(ncid = NCID, name = 'y', len
>> = ny, dimid = y_dimID))
>> call handle_err_nc_s(nf90_def_var(ncid = NCID, name = 'x',
>> xtype = NF90_FLOAT, dimids = x_dimID, varid = x_varID))
>> call handle_err_nc_s(nf90_def_var(ncid = NCID, name = 'y',
>> xtype = NF90_FLOAT, dimids = y_dimID, varid = y_varID))
>> call handle_err_nc_s(nf90_def_var(ncid = NCID, name = 'array',
>> xtype = NF90_FLOAT, dimids = (/x_dimID,y_dimID/), varid = array_varID))
>> call handle_err_nc_s(nf90_enddef(ncid = NCID))
>> call handle_err_nc_s(nf90_put_var(ncid = NCID, varid = array_varID,
>> values = array))
>> call handle_err_nc_s(nf90_close(ncid = NCID))
>>
>> end subroutine nc_out_s
>> !=====================================================================
>> subroutine nc_out_p
>>
>> ! use atham_module, only : mycoor,my_cart,myrank
>> ! use atham_module, only : nx,ny
>>
>> ! precompiler problem: include file had to be hardwired into code
>> ! include 'pnetcdf.inc'
>> #include "/usr/local/parallel-netcdf-1.0.1/include/pnetcdf.inc"
>>
>>
>> integer(kint) :: NCIDP
>> character(len=10) :: name
>> integer(kint) ::
>> x_dimIDP,y_dimIDP,x_varIDP,y_varIDP,array_varIDP
>> integer(kint) :: dimIDPs(2)
>> integer(MPI_OFFSET_KIND) :: dimlenx,dimleny
>> integer(MPI_OFFSET_KIND) :: istart,jstart
>> integer(MPI_OFFSET_KIND) :: ntx_mpi,nty_mpi,nx_mpi,ny_mpi
>>
>> ntx_mpi = ntx
>> nty_mpi = nty
>> nx_mpi = nx
>> ny_mpi = ny
>>
>> name = 'out_glb.nc'
>> istart = mycoor(1)*nx_mpi + 1
>> jstart = mycoor(2)*ny_mpi + 1
>>
>>
>> call
>> handle_err_nc_p(nfmpi_create(my_cart,'./'//trim(name),0,MPI_INFO_NULL,NCIDP),0)
>>
>>
>> call handle_err_nc_p(nfmpi_def_dim(NCIDP,'x',ntx_mpi,x_dimIDP),1)
>>
>> call handle_err_nc_p(nfmpi_def_dim(NCIDP,'y',nty_mpi,y_dimIDP),2)
>>
>> call
>> handle_err_nc_p(nfmpi_def_var(NCIDP,'x',NF_FLOAT,1,x_dimIDP,x_varIDP),3)
>>
>> call
>> handle_err_nc_p(nfmpi_def_var(NCIDP,'y',NF_FLOAT,1,y_dimIDP,y_varIDP),4)
>>
>> call
>> handle_err_nc_p(nfmpi_def_var(NCIDP,'array',NF_FLOAT,2,(/x_dimIDP,y_dimIDP/),array_varIDP),5)
>>
>>
>> call handle_err_nc_p(nfmpi_enddef(NCIDP),6)
>>
>> if (myrank==0) then
>> print *, ''
>> print *, 'array in nc_out_p subroutine'
>> print *, array
>> print *, ''
>> print *, '(/nx,ny/) ',(/nx,ny/)
>> print *, '(/istart,jstart/) ',(/istart,jstart/)
>> endif
>>
>> if (myrank==0) then
>> print *, ''
>> call handle_err_nc_p(nfmpi_inq_dimlen(NCIDP,x_dimIDP,dimlenx),7)
>> call handle_err_nc_p(nfmpi_inq_dimlen(NCIDP,y_dimIDP,dimleny),8)
>> print *, 'pnetCDF var dimensions in xy ',dimlenx,dimleny
>> print *, 'ntx_mpi ',ntx_mpi
>> print *, 'ntx ',ntx
>> print *, 'nx_mpi ',nx_mpi
>> print *, 'nx ',nx
>> endif
>>
>> ! this should have been possible in collective data mode, but did
>> not seem to work !!!
>> call handle_err_nc_p(nfmpi_begin_indep_data(NCIDP),9)
>>
>> call
>> handle_err_nc_p(nfmpi_put_vara_real(NCIDP,array_varIDP,(/istart,jstart/),(/nx_mpi,ny_mpi/),reshape(array,(/nx_mpi,ny_mpi/))),10)
>>
>>
>> call handle_err_nc_p(nfmpi_end_indep_data(NCIDP),11)
>>
>> call handle_err_nc_p(nfmpi_close(NCIDP),12)
>>
>> end subroutine nc_out_p
>> !=====================================================================
>> subroutine handle_err_nc_s(status)
>>
>> use netcdf
>>
>> integer, intent (in ) :: status
>>
>> if ( status /= nf90_noerr ) then
>> print *, ""
>> print *, trim(nf90_strerror(status))
>> print *, "atham_netcdf error: program execution stopped"
>> print *, ""
>> stop
>> end if
>> return
>> end subroutine handle_err_nc_s
>> !=====================================================================
>> subroutine handle_err_nc_p(status,callID)
>>
>> ! use atham_module, only : myrank
>>
>> #include "/usr/local/parallel-netcdf-1.0.1/include/pnetcdf.inc"
>>
>> integer, intent (in ) :: status
>> integer, intent (in ) :: callID
>>
>> if ( status /= nf_noerr ) then
>> if (myrank==0) then
>> print *, ""
>> print *, trim(nfmpi_strerror(status))
>> print *, "atham_pnetcdf error: program execution stopped
>> ",callID
>> print *, ""
>> ! stop
>> end if
>> endif
>> return
>> end subroutine handle_err_nc_p
>> !=====================================================================
>> subroutine exec_stop(MESSAGE,ERRORCLASS)
>>
>> character(len=*), intent(in ) :: MESSAGE
>> integer(kint), intent(in ) :: ERRORCLASS
>>
>> integer(kint) :: ierror
>>
>> print *, 'ERROR: MPI execution stopped'
>> print *, 'error class: ', ERRORCLASS
>> print *, 'message: ', trim(MESSAGE)
>> call mpi_abort(mpi_comm_world,ierror)
>> stop
>>
>> end subroutine exec_stop
>> !=====================================================================
>> subroutine makename(prefix,icoor,name)
>> !------------------------------------------------------------------!
>> ! create character name=prefix_(icoor(1)+1)_(icoor(2)+1) !
>> !------------------------------------------------------------------!
>> character(len=*), intent(in) :: prefix
>> integer(kint), dimension(2), intent(in) :: icoor
>> character(len=*), intent(out) :: name
>> !------------------------------------------------------------------!
>> ! local variables !
>> !------------------------------------------------------------------!
>> character(len=25) :: numx,numy
>> integer(kint) :: lenx,leny,lenf
>>
>> write(numx,'(i5)') icoor(1)+1
>> numx=adjustl(numx)
>> lenx=len_trim(numx)
>>
>> write(numy,'(i5)') icoor(2)+1
>> numy=adjustl(numy)
>> leny=len_trim(numy)
>>
>> name=adjustl(prefix)
>> lenf=len_trim(name)
>>
>> name=name(1:lenf)//'_'//numx(1:lenx)//'_'//numy(1:leny)
>>
>> end subroutine makename
>> !=====================================================================
>>
>> end program mpi_setup_netCDFout
>> # MAKEFILE used for development module testing
>> # Alex Hoffmann
>>
>> # set netCDF path
>> NETCDF = /usr/local/netcdf-3.6.2-ifort
>> PNETCDF= /usr/local/parallel-netcdf-1.0.1
>>
>> # set compiler
>> F90 = mpif90 # mpif90 refers to a shell script that includes links
>> to the libraries
>>
>>
>> # set compiler flags
>> OFLAGS = -g
>> # set include flags
>> IFLAGS = -I$(NETCDF)/include/ -I$(PNETCDF)/include/
>> # set pre-processing flags
>> #PFLAGS = -fpp
>> # concatenate flags
>> FLAGS = $(OFLAGS) $(IFLAGS) $(PFLAGS)
>> # set library flags
>> LFLAGS = -L$(NETCDF)/lib -lnetcdf -L$(PNETCDF)/lib -lpnetcdf
>>
>> RM = rm -f
>>
>> #MAIN =
>> MAIN = mpi_setup_netCDFout
>> OBJECTS =
>>
>>
>> all: $(MAIN)
>>
>>
>> $(MAIN): $(MAIN).F90 $(OBJECTS)
>> $(F90) $(FLAGS) -o $@ $^ $(LFLAGS)
>>
>> clean:
>> $(RM) *.o *.mod
>
--
______________________________________________________________________
Alex HOFFMANN PhD Candidate
Centre for Atmospheric Science
Department of Geography University of Cambridge
Downing Place, Cambridge CB2 3EN, UK
e-mail: ah519 at cam.ac.uk tel: +44 (0)1223 766581
www.geog.cam.ac.uk/people/hoffmann/
______________________________________________________________________
More information about the parallel-netcdf
mailing list