!################################################################# ! ! contains routines to read emissions and ! to read and write the main model state from/to file ! !### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module io_save use GO, only : gol, goErr, goPr use GO, only : TrcFile, InitRc => Init, DoneRc => Done, ReadRc use global_data, only : rcfile implicit none ! --- in/out -------------------------- private public :: readhdf, readhdfmmr, savehdf, readhdfmmix public :: readtm3hdf, readtm3hdf_4d ! --- const ---------------------------------- character(len=*), parameter :: mname = 'io_save' type(TrcFile) :: rcF contains ! ====================================================================== !WP! this subroutine reads HDF emissions files on ! 1x1 degree resolution and returns the specified field. !WP! input for this subroutine are: !WP! !WP! filename --> name of hdf file to be opened or read !MK! rank --> rank of the dataset (1,2,3) !WP! im_emis --> x-DIMENSION of requested field !WP! jm_emis --> y-DIMENSION of requested field !WP! lm_emis --> z-DIMENSION of requested field !WP! month --> the month to read (0 for january or simply the first field) !WP! field --> the field (im_emis,jm_emis,lm_emis) to put the values in !WP! field_name --> the name of the field to match the one in the HDF file !WP! !WP! The routine checks the dimensions and reports the succes/failure. subroutine readtm3hdf( filename, rank, im_emis, jm_emis, lm_emis, & month, field, field_name, idt ) use file_hdf, only : THdfFile, Init, Done use io_hdf , only : FAIL use io_hdf , only : io_read2d_32g, io_read3d_32g use toolbox , only : escape_tm ! --- in/out -------------------------------------- character(len=*),intent(in) :: filename character(len=*),intent(in) :: field_name integer,intent(in) :: rank integer,intent(in) :: month integer,intent(in) :: im_emis, jm_emis, lm_emis real,dimension(im_emis,jm_emis,lm_emis),intent(out) :: field integer,dimension(6),optional :: idt ! --- local ---------------------------------------- type(THdfFile) :: hdf integer :: status integer :: idx real, allocatable :: field1d(:) ! first dim will be 1 real, allocatable :: field2d(:,:) real, allocatable :: field3d(:,:,:) ! --- begin ------------------------------------------- ! open hdf file: call Init( hdf, trim(filename), 'read', status ) if (status/=0) call escape_tm('ERROR in readtm3hdf') !FD hdf files indices run from 0 ...n-1 idx = month select case(rank) case(1) ! latitudional field allocate(field1d(jm_emis)) if(present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,idate=idt) if(.not.present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field1d,field_name, status,index=idx) if ( status /= 0 ) then print *, 'readtm3hdf: io_read2d_32g failed:' print *, 'readtm3hdf: filename : ', trim(filename) print *, 'readtm3hdf: rank : ', rank print *, 'readtm3hdf: im_emis : ', im_emis print *, 'readtm3hdf: jm_emis : ', jm_emis print *, 'readtm3hdf: lm_emis : ', lm_emis print *, 'readtm3hdf: month : ', month print *, 'readtm3hdf: field_name : ', trim(field_name) if (present(idt)) print *, 'readtm3hdf: idt : ', idt call escape_tm('readtm3hdf: Error reading HDF file 1D') end if field(1,:,1)=field1d deallocate(field1d) case(2) ! 2D surface field allocate(field2d(im_emis,jm_emis)) if(present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field2d,field_name, status,idate=idt) if(.not.present(idt)) call io_read2d_32g(hdf%id,im_emis,jm_emis,field2d,field_name, status,index=idx) if ( status /= 0 ) then print *, 'readtm3hdf: io_read2d_32g failed:' print *, 'readtm3hdf: filename : ', trim(filename) print *, 'readtm3hdf: rank : ', rank print *, 'readtm3hdf: im_emis : ', im_emis print *, 'readtm3hdf: jm_emis : ', jm_emis print *, 'readtm3hdf: lm_emis : ', lm_emis print *, 'readtm3hdf: month : ', month print *, 'readtm3hdf: field_name : ', trim(field_name) if (present(idt)) print *, 'readtm3hdf: idt : ', idt call escape_tm('readtm3hdf: Error reading HDF file 2D') endif field(:,:,1)=field2d deallocate(field2d) case (3) ! 3D field allocate(field3d(im_emis,jm_emis,lm_emis)) if(present(idt)) call io_read3d_32g(hdf%id,im_emis,jm_emis,lm_emis,field3d,field_name, status,idate=idt) if(.not.present(idt)) call io_read3d_32g(hdf%id,im_emis,jm_emis,lm_emis,field3d,field_name, status,index=idx) if ( status /= 0 ) then print *, 'readtm3hdf: io_read2d_32g failed:' print *, 'readtm3hdf: filename : ', trim(filename) print *, 'readtm3hdf: rank : ', rank print *, 'readtm3hdf: im_emis : ', im_emis print *, 'readtm3hdf: jm_emis : ', jm_emis print *, 'readtm3hdf: lm_emis : ', lm_emis print *, 'readtm3hdf: month : ', month print *, 'readtm3hdf: field_name : ', trim(field_name) if (present(idt)) print *, 'readtm3hdf: idt : ', idt call escape_tm('readtm3hdf: Error reading HDF file 3D') endif field(:,:,:)=field3d deallocate(field3d) case (23) ! 2D field lat-pres allocate(field2d(jm_emis,lm_emis)) if(present(idt)) call io_read2d_32g(hdf%id,jm_emis,lm_emis,field2d,field_name, status,idate=idt) if(.not.present(idt)) call io_read2d_32g(hdf%id,jm_emis,lm_emis,field2d,field_name, status,index=idx) if ( status /= 0 ) then !might be due to (1,jm,lm) file! print *, 'readtm3hdf: try shape : ', 1, jm_emis, lm_emis allocate(field3d(1,jm_emis,lm_emis)) if ( present(idt) ) then call io_read3d_32g(hdf%id,1,jm_emis,lm_emis,field3d,field_name, status,idate=idt) else call io_read3d_32g(hdf%id,1,jm_emis,lm_emis,field3d,field_name, status,index=idx) end if if ( status /= 0 ) then print *, 'readtm3hdf: io_read2d_32g failed:' print *, 'readtm3hdf: filename : ', trim(filename) print *, 'readtm3hdf: rank : ', rank print *, 'readtm3hdf: im_emis : ', im_emis print *, 'readtm3hdf: jm_emis : ', jm_emis print *, 'readtm3hdf: lm_emis : ', lm_emis print *, 'readtm3hdf: month : ', month print *, 'readtm3hdf: field_name : ', trim(field_name) if (present(idt)) print *, 'readtm3hdf: idt : ', idt call escape_tm('readtm3hdf: Error reading HDF file 23') end if field(:,:,:)=field3d deallocate(field3d) else field(1,:,:)=field2d end if deallocate(field2d) case default write (gol,*) 'readtm3hdf: Please check the rank of the required field,', & ' must be 1,2,3, or 23'; call goPr end select ! close: call Done( hdf, status ) if (status/=0) call escape_tm('ERROR in readtm3hdf') !if ( idx >= 0 ) then ! write (gol,*) 'readtm3hdf: ',trim(filename),' ; ds#= ',idx+1; call goPr ! write (gol,*) 'readtm3hdf: total month sum after reading the full field', & ! sum(field)/1e9 ,' (x1e9) '; call goPr !end if end subroutine readtm3hdf ! *** subroutine readtm3hdf_4d( filename, rank, im_emis, jm_emis, lm_emis, km_emis, & month, field, field_name, idt ) use file_hdf, only : THdfFile, Init, Done use io_hdf , only : FAIL use io_hdf , only : io_read4d_32g use toolbox , only : escape_tm ! --- in/out -------------------------------------- character(len=*),intent(in) :: filename character(len=*),intent(in) :: field_name integer,intent(in) :: rank integer,intent(in) :: month integer,intent(in) :: im_emis, jm_emis, lm_emis,km_emis real,dimension(im_emis,jm_emis,lm_emis),intent(out) :: field integer,dimension(6),optional :: idt ! --- local ---------------------------------------- type(THdfFile) :: hdf integer :: status integer :: idx real, allocatable :: field1d(:) ! first dim will be 1 real, allocatable :: field2d(:,:) real, allocatable :: field3d(:,:,:) real, allocatable :: field4d(:,:,:,:) ! --- begin ------------------------------------------- ! open hdf file: call Init( hdf, trim(filename), 'read', status ) if (status/=0) call escape_tm('ERROR in readtm3hdf_4d') !FD hdf files indices run from 0 ...n-1 idx = month select case(rank) case (4) ! 4D field allocate(field4d(im_emis,jm_emis,lm_emis,km_emis)) if(present(idt)) call io_read4d_32g(hdf%id,im_emis,jm_emis,lm_emis,km_emis,field4d,field_name, status,idate=idt) if(.not.present(idt)) call io_read4d_32g(hdf%id,im_emis,jm_emis,lm_emis,km_emis,field4d,field_name, status,index=idx) if ( status /= 0 ) then print *, 'readtm3hdf_4d: io_read4d_32g failed:' print *, 'readtm3hdf_4d: filename : ', trim(filename) print *, 'readtm3hdf_4d: rank : ', rank print *, 'readtm3hdf_4d: im_emis : ', im_emis print *, 'readtm3hdf_4d: jm_emis : ', jm_emis print *, 'readtm3hdf_4d: lm_emis : ', lm_emis print *, 'readtm3hdf_4d: month : ', month print *, 'readtm3hdf_4d: field_name : ', trim(field_name) if (present(idt)) print *, 'readtm3hdf_4d: idt : ', idt call escape_tm('readtm3hdf_4d: Error reading HDF file 3D') endif field(:,:,:)=field4d(:,:,:,month) deallocate(field4d) case default write (gol,*) 'readtm3hdf_4d: Please check the rank of the required field,', & ' must be 4'; call goPr end select ! close: call Done( hdf, status ) if (status/=0) call escape_tm('ERROR in readtm3hdf_4d') end subroutine readtm3hdf_4d ! *** ! ! data are expected in kg ! (typical save.hdf files from previous run...) ! only long-lived tracer are read, since the fields ! for chemistry are not yet allocated ! subroutine ReadHDF( region, file_name ) use io_hdf use dims, only : region_name, im, jm, lm, nregions, & datadir, parent use global_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : ntrace, ntracet use toolbox, only : escape_tm use ParTools , only : myid, root_t use ParTools , only : ntracet_ar, ntracetloc #ifdef MPI use mpi_comm, only : scatter_after_read_t use mpi_const, only : com_trac, ierr, my_real, mpi_integer #endif ! --- in/out ------------------------------ integer, intent(in) :: region character(len=*), intent(in) :: file_name ! --- local ------------------------------- real,dimension(:,:,:,:),pointer :: rm real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #ifdef secmom real,dimension(:,:,:,:),pointer :: rxxm,rxym,rxzm,ryym,ryzm,rzzm #endif real,dimension(:,:,:),pointer :: m integer :: io, sfstart, ifail, sds_id integer :: n_datasets, n_file_attrs integer :: istat, attributes, num_type integer :: sffinfo, sfselect, sfginfo integer :: sfend, sfrnatt, sfrcatt, sffattr integer,dimension(6) :: idate_save integer,dimension(nregions) :: im_file,jm_file,lm_file integer :: ntrace_file, i,j,l,n character(len=80) :: msg_file integer :: ind, my_parent,imr,jmr,lmr real,dimension(:,:,:,:), allocatable :: field, field2 real,dimension(:,:,:), allocatable :: msave integer :: from_file ! signals reading from file ! --- begin -------------------------------- m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #ifdef secmom rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t #endif #endif imr = im(region) ; jmr = jm(region) ; lmr = lm(region) allocate(field(-1:imr+2,-1:jmr+2,lmr,ntracet)) field = 0.0 allocate(msave(imr,jmr,lmr)) if ( myid == root_t ) then from_file = 1 allocate(field2(imr,jmr,lmr,ntracet)) ind = len_trim(datadir) !io = sfstart(datadir(1:ind)//file_name//trim(region_name(region)), DFACC_READ) io = sfstart( file_name, DFACC_READ ) if ( io == FAIL ) then my_parent = parent(region) if ( my_parent == 0 ) then !call escape_tm('readhdf: Could not open file and no parent '// & ! 'available :'//datadir(1:ind)//file_name//trim(region_name(region))) call escape_tm('readhdf: Could not open file and no parent '// & 'available :'//trim(file_name) ) else from_file = 0 write (*,'("readhdf: Trying to initialise ",a," from parent....")') trim(region_name(region)) end if end if if ( from_file == 1 ) then write (*,'("readhdf:")') write (*,'("readhdf: opened for READ access: ",a)') trim(file_name) call io_read3d_32g(io,imr,jmr,lmr,msave,'m',ifail) if ( ifail /= 0 ) then call escape_tm('readhdf: Failed to read m from saveold.hdf in readhdf') end if call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rm',ifail) if( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rm',ifail) end if if ( ifail /= 0 ) & call escape_tm('readhdf: Failed to read in saveold.hdf in readhdf') write (*,'("readhdf: Read rm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! end if !from_file end if !root_t #ifdef MPI ! broadcast from_file ! call mpi_bcast(from_file, 1, mpi_integer , root_t, com_trac, ierr) #endif if ( from_file == 1 ) then #ifdef MPI call mpi_bcast(msave, imr*jmr*lmr, my_real, root_t, com_trac, ierr) call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rm,root_t) #else rm = field #endif else call init_child(region) ! get info from parent.... end if field = 0.0 if ( from_file == 1 ) then #ifdef slopes if ( myid == root_t ) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rxm',ifail) end if if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rxm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rxm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! end if !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxm,root_t) #else rxm = field #endif field = 0.0 if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rym',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rym',ifail) end if if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rym in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rym, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! end if !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rym,root_t) #else rym = field #endif field = 0.0 if ( myid == root_t ) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rzm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rzm',ifail) end if if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rzm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rzm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! deallocate(field2) end if !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rzm,root_t) #else rzm = field #endif #endif ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! >>> second moments ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #ifdef secmom allocate(field2(imr,jmr,lmr,ntracet)) if ( myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxxm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rxxm',ifail) endif if ( ifail /= 0 ) then call escape_tm('readhdf: Failed to read in rxxm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rxxm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxxm,root_t) #else rxxm = field #endif field = 0.0 if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxym',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rxym',ifail) endif if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rxym in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rxym, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxym,root_t) #else rxym = field #endif field = 0.0 if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rxzm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rxzm',ifail) endif if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rxzm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rxzm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rxzm,root_t) #else rxzm = field #endif if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'ryym',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'ryym',ifail) endif if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in ryym in saveold.hdf in readhdf') end if write (*,'("readhdf: Read ryym, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,ryym,root_t) #else ryym = field #endif field = 0.0 if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'ryzm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'ryzm',ifail) endif if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in ryzm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read ryzm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,ryzm,root_t) #else ryzm = field #endif field = 0.0 if(myid == root_t) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field2,'rzzm',ifail) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field2,'rzzm',ifail) endif if ( ifail /= 0 ) then call escape_tm( & 'readhdf: Failed to read in rzzm in saveold.hdf in readhdf') end if write (*,'("readhdf: Read rzzm, ifail = ",i6)') ifail field(1:imr,1:jmr,1:lmr,:)=field2 ! only transported tracers! deallocate(field2) endif !root_t #ifdef MPI call scatter_after_read_t(field, imr,jmr,lmr,2,2,0,rzzm,root_t) #else rzzm = field #endif #endif ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! <<< end second moments ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< if ( myid == root_t ) then istat = sfend(io) if ( istat /= FAIL ) then write (*,'("readhdf: ... file closed")') write (*,'("readhdf: ")') else call escape_tm('readhdf: ERROR in restart from HDF file') end if end if ! scale in input by a possible factor. ! rm/msave ---> mixing ratio on file. mr*m --> something back in kilo's if ( myid == root_t ) then print *, 'readhdf: Maximum ratio of the air masses:', & maxval(abs(msave/m(1:imr,1:jmr,1:lmr))) print *, 'readhdf: Minimum ratio of the air masses:', & minval(abs(msave/m(1:imr,1:jmr,1:lmr))) end if do n=1,ntracetloc do l = 1,lmr do j = 1,jmr do i = 1,imr rm(i,j,l,n) = rm (i,j,l,n)*m(i,j,l)/msave(i,j,l) #ifdef slopes rxm(i,j,l,n) = rxm(i,j,l,n)*m(i,j,l)/msave(i,j,l) rym(i,j,l,n) = rym(i,j,l,n)*m(i,j,l)/msave(i,j,l) rzm(i,j,l,n) = rzm(i,j,l,n)*m(i,j,l)/msave(i,j,l) #ifdef secmom rxxm(i,j,l,n) = rxxm(i,j,l,n)*m(i,j,l)/msave(i,j,l) rxym(i,j,l,n) = rxym(i,j,l,n)*m(i,j,l)/msave(i,j,l) rxzm(i,j,l,n) = rxzm(i,j,l,n)*m(i,j,l)/msave(i,j,l) ryym(i,j,l,n) = ryym(i,j,l,n)*m(i,j,l)/msave(i,j,l) ryzm(i,j,l,n) = ryzm(i,j,l,n)*m(i,j,l)/msave(i,j,l) rzzm(i,j,l,n) = rzzm(i,j,l,n)*m(i,j,l)/msave(i,j,l) #endif #endif end do end do end do end do end if ! from_file == 1 deallocate(msave) deallocate(field) nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #ifdef secmom nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) #endif #endif end subroutine ReadHDF ! *** subroutine init_child(region) ! ! Parent is used to initilise rm, rxm,rym,rzm of region. ! Assumed parallel over tracers... ! use dims, only : lm, ibeg, iend, jbeg, jend, xref, yref, parent use tracer_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : ntrace, ntracet use toolbox, only : escape_tm use ParTools, only : myid, ntracetloc, root_t implicit none ! in/out integer, intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:),pointer :: m real,dimension(:,:,:,:),pointer :: rmp !parent real,dimension(:,:,:),pointer :: mp real :: mm_parent integer :: my_parent integer :: n,l,jp,ip,ic,jc, xref_, yref_, i,j ! start m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif my_parent = parent(region) if ( my_parent == 0 ) then call escape_tm( 'init_child: Region 1 cannot be initialised') end if mp => m_dat(my_parent)%data rmp => mass_dat(my_parent)%rm_t xref_ = xref(region)/xref(my_parent) yref_ = yref(region)/yref(my_parent) do n=1,ntracetloc do l=1,lm(region) do jp=jbeg(region), jend(region) jc = (jp-jbeg(region))*yref_ do ip=ibeg(region), iend(region) mm_parent = rmp(ip,jp,l,n)/mp(ip,jp,l) ic = (ip-ibeg(region))*xref_ do j=1,yref_ do i=1,xref_ rm (ic+i,jc+j,l,n) = mm_parent*m(ic+i,jc+j,l) #ifdef slopes rxm(ic+i,jc+j,l,n) = 0.0 rym(ic+i,jc+j,l,n) = 0.0 rzm(ic+i,jc+j,l,n) = 0.0 #endif end do end do end do end do end do end do if ( myid == root_t ) print *, 'readhdf: Initialisation successful ' nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif nullify(mp) nullify(rmp) end subroutine init_child subroutine init_child_edge(region) ! ! ONLY the edges! ! Parent is used to initilise rm, rxm,rym,rzm of region. ! Assumed parallel over tracers... ! use dims, only : lm, ibeg, iend, jbeg, jend, xref, yref, parent, im, jm use global_data, only : region_dat use tracer_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : ntrace, ntracet use toolbox, only : escape_tm use ParTools, only : myid, ntracetloc, root_t implicit none ! in/out integer, intent(in) :: region ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:),pointer :: m real,dimension(:,:,:,:),pointer :: rmp !parent real,dimension(:,:,:),pointer :: mp integer,dimension(:,:),pointer :: edge real :: mm_parent integer :: my_parent integer :: n,l,jp,ip,ic,jc, xref_, yref_, i,j ! start my_parent = parent(region) if ( my_parent == 0 ) return m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif edge => region_dat(region)%edge mp => m_dat(my_parent)%data rmp => mass_dat(my_parent)%rm_t xref_ = xref(region)/xref(my_parent) yref_ = yref(region)/yref(my_parent) do n=1,ntracetloc do l=1,lm(region) do j=1,jm(region) jp = jbeg(region) + (j-1)/yref_ do i=1,im(region) if (edge(i,j) /= -1 ) cycle ip = ibeg(region) + (i-1)/xref_ mm_parent = rmp(ip,jp,l,n)/mp(ip,jp,l) rm (i,j,l,n) = mm_parent*m(i,j,l) #ifdef slopes rxm(i,j,l,n) = 0.0 rym(i,j,l,n) = 0.0 rzm(i,j,l,n) = 0.0 #endif end do end do end do end do if ( myid == root_t ) print *, 'readhdf: Edges Initialisation successful ' nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif nullify(mp) nullify(rmp) end subroutine init_child_edge subroutine readhdfmmr(region,file_name) ! ! data are expected in mixing ratio... ! use dims, only : im, jm, lm, nregions, datadir, region_name use io_hdf use tracer_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : fscale, ntracet use toolbox, only : escape_tm #ifdef MPI use mpi_comm, only : scatter_after_read_t use mpi_const, only : root, myid, ntracet_ar, ntracetloc #endif implicit none ! in/out integer, intent(in) :: region character(len=*), intent(in) :: file_name ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #ifdef secmom real,dimension(:,:,:,:),pointer :: rxxm,rxym,rxzm,ryym,ryzm,rzzm #endif #endif real,dimension(:,:,:),pointer :: m integer :: io, sfstart, ifail, sds_id integer :: n_datasets, n_file_attrs integer :: istat, attributes, num_type integer :: sffinfo, sfselect, sfginfo integer :: sfend, sfrnatt, sfrcatt, sffattr integer,dimension(6) :: idate_save integer,dimension(nregions) :: im_file,jm_file,lm_file integer :: ntrace_file character(len=80) :: msg_file integer :: ind,n,i,j,l,offsetn real,dimension(:,:,:,:), allocatable :: field,field2 ! start m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #ifdef secmom rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t #endif #endif allocate(field(-1:im(region)+2,-1:jm(region)+2,1:lm(region),1:ntracet)) field=0.0 #ifdef MPI if ( myid == root ) then #endif ! to hold all ntrace values of rm allocate(field2(im(region),jm(region),lm(region),ntracet)) io = sfstart( file_name, DFACC_READ ) if ( io == -1 ) call escape_tm('readhdfmmr: Could not open file:'//file_name) ind = 1 ! select 1st dataset....(0 = m) istat = sffinfo(io, n_datasets, n_file_attrs) print*, 'readhdfmmr: Number of datasets and attributes', & n_datasets, n_file_attrs call io_read4d_32g(io,im(region),jm(region),lm(region),ntracet, & field2,'rm',ifail,index=ind) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,im(region),jm(region),lm(region), & field2,'rm',ifail,index=ind) end if if ( ifail /= 0 ) then print *, 'readhdfmmr: '//trim(file_name) call escape_tm('readhdfmmr: Failed to read fields') end if print*, 'readhdfmmr: Read rm, ifail = ', ifail istat = sfend(io) if (istat /= FAIL) then print*,'readhdfmmr: ... file closed' print*,' ' else call escape_tm('readhdfmmr: ERROR in restart from HDF file') end if ! only transported tracers! field(1:im(region),1:jm(region),1:lm(region),:)=field2 deallocate(field2) ! no longer needed #ifdef MPI end if ! root call scatter_after_read_t(field,im(region),jm(region),lm(region),& 2,2,0,rm,root) #else rm = field #endif #ifdef MPI offsetn = sum(ntracet_ar(0:myid-1)) do n=1,ntracetloc #else offsetn = 0 do n=1,ntracet #endif do l = 1,lm(region) do j = 1,jm(region) do i = 1,im(region) rm(i,j,l,n) = rm(i,j,l,n)*m(i,j,l)/fscale(offsetn+n) #ifdef slopes rxm(i,j,l,n) = 0.0 rym(i,j,l,n) = 0.0 rzm(i,j,l,n) = 0.0 #ifdef secmom rxxm(i,j,l,n) = 0.0 rxym(i,j,l,n) = 0.0 rxzm(i,j,l,n) = 0.0 ryym(i,j,l,n) = 0.0 ryzm(i,j,l,n) = 0.0 rzzm(i,j,l,n) = 0.0 #endif #endif end do end do end do end do deallocate(field) nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #ifdef secmom nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) #endif #endif !WP! checked rm with the values in savehdf after scatter ! and gather on 4 processors=>OK end subroutine readhdfmmr subroutine readhdfmmix(region,file_name) ! ! data are expected in mixing ratio from an mmix file ! use dims, only : im, jm, lm, nregions, datadir, region_name, parent use io_hdf use tracer_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : fscale, ntracet, names use toolbox, only : escape_tm #ifdef MPI use mpi_comm, only : scatter_after_read_t use mpi_const, only : ntracet_ar, ntracetloc use mpi_const, only : com_trac, mpi_integer, root_t, ierr #endif use ParTools, only : root, myid, root_t implicit none ! in/out integer, intent(in) :: region character(len=*), intent(in) :: file_name ! local real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:),pointer :: m integer :: io, sfstart, ifail, sds_id integer :: n_datasets, n_file_attrs integer :: istat, attributes, num_type integer :: sffinfo, sfselect, sfginfo integer :: sfend, sfrnatt, sfrcatt, sffattr integer,dimension(6) :: idate_save integer,dimension(nregions) :: im_file,jm_file,lm_file integer :: ntrace_file character(len=80) :: msg_file integer :: ind,n,i,j,l,offsetn, from_file, my_parent real,dimension(:,:,:,:), allocatable :: field real,dimension(:,:,:), allocatable :: field3d ! start m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #endif rm = 0.0 #ifdef slopes rxm = 0.0 rym = 0.0 rzm = 0.0 #endif allocate(field(-1:im(region)+2,-1:jm(region)+2,1:lm(region),1:ntracet)) field = 0.0 !09/2004 cmk: to initialise the edges! if ( myid == root_t ) then from_file = 1 allocate(field3d(im(region),jm(region),lm(region))) io = sfstart(file_name, DFACC_READ) if ( io == FAIL ) then my_parent = parent(region) if ( my_parent == 0 ) then call escape_tm('readhdfmmix : Could not open file and no parent '// & 'available :'//file_name ) else from_file = 0 print *, 'readhdf: Trying to initialise '//trim(region_name(region))// & ' from parent....' endif endif if ( from_file == 1 ) then ! read from file: print*,' ' print*,'readhdfmmix: ',file_name,'... opened for READ access.' istat = sffinfo(io, n_datasets, n_file_attrs) print*, 'readhdfmmix: Number of datasets and attributes', & n_datasets, n_file_attrs do n = 1, ntracet call io_read3d_32(io, im(region), jm(region), lm(region), field3d, names(n), ifail) if (ifail == 0) then field(1:im(region), 1:jm(region), 1:lm(region), n) = field3d else print *, 'Readhdfmmix:', names(n), ' not found in dataset --> 0.0 ' field(1:im(region), 1:jm(region), 1:lm(region), n) = 0.0 endif enddo istat = sfend(io) if (istat /= FAIL) then print*,'readhdfmmix: ... file closed' print*,' ' else call escape_tm('readhdfmmix: ERROR in restart from HDF file') end if deallocate(field3d) ! no longer needed endif ! from file endif ! root_t #ifdef MPI ! broadcast from_file ! call mpi_bcast(from_file, 1, mpi_integer , root_t, com_trac, ierr) #endif if ( from_file == 1 ) then #ifdef MPI call scatter_after_read_t(field,im(region),jm(region),lm(region),& 2,2,0,rm,root) #else rm = field #endif #ifdef MPI offsetn = sum(ntracet_ar(0:myid-1)) do n=1,ntracetloc #else offsetn = 0 do n=1,ntracet #endif do l = 1,lm(region) do j = 1,jm(region) do i = 1,im(region) rm(i,j,l,n) = rm(i,j,l,n)*m(i,j,l)/fscale(offsetn+n) ! from mmix to kg tracer #ifdef slopes rxm(i,j,l,n) = 0.0 rym(i,j,l,n) = 0.0 rzm(i,j,l,n) = 0.0 #endif end do end do end do end do deallocate(field) call init_child_edge(region) ! get edge info from parent.... else ! initialise from parent: call init_child(region) ! get info from parent.... end if nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #endif !WP! checked rm with the values in savehdf after scatter ! and gather on 4 processors=>OK end subroutine readhdfmmix !------------------------------------------------------------------ ! save all essential model parameters and fields on unit kdisk ! ! v 8.4 ! modified for HDF output !------------------------------------------------------------------ subroutine SaveHDF( msg, file_name, status ) use file_hdf , only : THdfFile, Init, Done, WriteAttribute, TSds, SetDim, WriteData use dims , only : splitsave use dims , only : nregions, region_name use dims , only : im, jm, lm use dims , only : dx, xref, xbeg, xend, ibeg, iend use dims , only : dy, yref, ybeg, yend, jbeg, jend use dims , only : dz, zref, zbeg, zend, lbeg, lend use dims , only : areag use dims , only : czeta, czetak use dims , only : itau, itaui, itaue, itau0, itaut use dims , only : idate, idatei, idatee, idate0, idatet use dims , only : icalendo, iyear0, julday0 use dims , only : newyr, newmonth, newday, newsrun use dims , only : tref use dims , only : ndyn, nconv, ndiag, nchem, nsrce use dims , only : nread, ninst, ncheck, ndiff, ndiagp1, ndiagp2, nstep use dims , only : nwrite use dims , only : istart use dims , only : cpu0, cpu1 use dims , only : xlabel use dims , only : cdebug, kmain, kdebug use dims , only : limits use dims , only : at, bt use dims , only : nsplitsteps, splitorder use tracer_data, only : mass_dat, chem_dat use meteodata , only : m_dat use chem_param , only : ntrace, ntrace_chem, ntracet, ra, names, fscale use chem_param , only : nstd, istd use io_hdf, only : io_write use datetime, only : tstamp use toolbox, only : escape_tm use partools , only : myid, root_t, root_k use partools , only : Par_Barrier, Par_Gather_Tracer_t use partools , only : Par_Gather_Tracer_k use partools , only : Par_Barrier use partools , only : which_par, previous_par use global_data, only : outdir ! --- in/out ---------------------------------- character(len=*), intent(in) :: file_name character(len=*), intent(in) :: msg integer, intent(out) :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/SaveHDF' ! --- local ---------------------------------- real, dimension(:,:,:,:), pointer :: rm #ifdef slopes real, dimension(:,:,:,:), pointer :: rxm,rym,rzm #ifdef secmom real, dimension(:,:,:,:), pointer :: rxxm,rxym,rxzm,ryym,ryzm,rzzm #endif #endif real, dimension(:,:,:,:), pointer :: rmc real, dimension(:,:,:), pointer :: m real, dimension(:,:,:,:), allocatable :: fieldglob4d type(THdfFile) :: hdf type(TSds) :: sds integer :: region, imr,jmr, lmr,i character(len=64) :: msg_file character(len=256) :: fname character(len=32 ) :: key ! --- begin ----------------------------------- if ( .not. splitsave ) call escape_tm('savehdf: splitsave should be true ') call InitRc( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) do region = 1,nregions if ( myid == root_t ) then ! try if filename for istart.3 was specified in rcfile: write (key,'("start.3.",a)') trim(region_name(region)) call ReadRc( rcF, key, fname, status, default='file_name_empty' ) if(status == -1) then ! no specific istart.3 file name found, use default name write (fname,'(a,"_",i4.4,3i2.2,"_",a,".hdf")') trim(file_name), idate(1:4), trim(region_name(region)) else write (gol,*) 'Using save file names from rc-file start.3.* values'; call goPr endif call Init( hdf, trim(fname), 'create', status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'itau' , itau , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nregions' , nregions , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'region_name', trim(region_name(region)), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'im' , im(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'jm' , jm(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'lm' , lm(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'dx' , dx/xref(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'dy' , dy/yref(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'dz' , dz/zref(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'xbeg' , xbeg(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'xend' , xend(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ybeg' , ybeg(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'yend' , yend(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'zbeg' , zbeg(region) , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'zend' , zend(region) , status ) IF_NOTOK_RETURN(status=1) if ( region /= 1 ) then call WriteAttribute( hdf, 'ibeg', ibeg(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'iend', iend(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'jbeg', jbeg(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'jend', jend(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'lbeg', lbeg(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'lend', lend(region), status ) IF_NOTOK_RETURN(status=1) end if call WriteAttribute( hdf, 'xref' , xref(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'yref' , yref(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'zref' , zref(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'tref' , tref(region), status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ntracet' , ntracet , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ntrace ' , ntrace , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nstd' , nstd , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'idate' , idate , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'istart' , istart , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ndyn' , ndyn , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nconv' , nconv , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ndiag' , ndiag , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nchem' , nchem , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nsrce' , nsrce , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nread' , nread , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nwrite' , nwrite , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ninst' , ninst , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ncheck' , ncheck , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ndiff' , ndiff , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'itaui' , itaui , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'itaue' , itaue , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'itaut' , itaut , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'itau0' , itau0 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'idatei' , idatei , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'idatee' , idatee , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'idatet' , idatet , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'idate0' , idate0 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'icalendo' , icalendo , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'iyear0' , iyear0 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'julday0' , julday0 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ndiagp1' , ndiagp1 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ndiagp2' , ndiagp2 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nstep' , nstep , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'cpu0' , cpu0 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'cpu1' , cpu1 , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'ra' , ra , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'fscale' , fscale , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'names' , names , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'areag' , areag , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'czeta' , czeta , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'czetak' , czetak , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'xlabel' , xlabel , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'istd' , istd , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'newyr' , newyr , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'newmonth' , newmonth , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'newday' , newday , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'newsrun' , newsrun , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'cdebug' , cdebug , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'limits' , limits , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'at' , at , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'bt' , bt , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'nsplitsteps', nsplitsteps , status ) IF_NOTOK_RETURN(status=1) call WriteAttribute( hdf, 'splitorder' , splitorder , status ) IF_NOTOK_RETURN(status=1) msg_file = msg call WriteAttribute( hdf, 'msg', msg_file, status ) IF_NOTOK_RETURN(status=1) end if ! root all PEs from here: call Par_Barrier m => m_dat(region)%data rm => mass_dat(region)%rm_t #ifdef slopes rxm => mass_dat(region)%rxm_t rym => mass_dat(region)%rym_t rzm => mass_dat(region)%rzm_t #ifdef secmom rxxm => mass_dat(region)%rxxm_t rxym => mass_dat(region)%rxym_t rxzm => mass_dat(region)%rxzm_t ryym => mass_dat(region)%ryym_t ryzm => mass_dat(region)%ryzm_t rzzm => mass_dat(region)%rzzm_t #endif #endif which_par=previous_par(region) if ( which_par /= 'tracer' ) & call escape_tm('savehdf: Data should be parallel over tracers') imr = im(region) ; jmr = jm(region) ; lmr = lm(region) ! allocate global array on all PEs allocate(fieldglob4d(-1:imr+2,-1:jmr+2,lmr, ntracet)) fieldglob4d = 0.0 ! false signals: do not broadcast call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rm,.false.) call Par_Barrier if ( myid == root_t ) then call Init( sds, hdf, 'm', shape(m(1:imr,1:jmr,1:lmr)), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call WriteData(sds,m(1:imr,1:jmr,1:lmr),status) call Done(sds,status) IF_ERROR_RETURN(status=1) call Init( sds, hdf, 'rm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier #ifdef slopes ! false signals: do not broadcast call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.) if ( myid == root_t ) then call Init( sds, hdf, 'rxm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier ! false signals: do not broadcast call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rym,.false.) if ( myid == root_t ) then call Init( sds, hdf, 'rym', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier ! false signals: do not broadcast call Par_Gather_Tracer_t(fieldglob4d,imr,jmr,lmr,2,2,0,ntracet,rzm ,.false.) if ( myid == root_t ) then call Init( sds, hdf, 'rzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if #ifdef secmom ! second moments call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxxm ,.false.) ! false signals: do not broadcast if ( myid == root_t ) then call Init( sds, hdf, 'rxxm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxym ,.false.) ! false signals: do not broadcast if ( myid == root_t ) then call Init( sds, hdf, 'rxym', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rxzm ,.false.) ! false signals: do not broadcast if ( myid == root_t ) then call Init( sds, hdf, 'rxym', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, ryym ,.false.) ! false signals: do not broadcast if ( myid == root_t) then call Init( sds, hdf, 'ryym', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, ryzm ,.false.) ! false signals: do not broadcast if ( myid == root_t) then call Init( sds, hdf, 'ryzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) endif call Par_Barrier call Par_Gather_Tracer_t(fieldglob4d,imr,jmr, lmr, 2,2,0,ntracet, rzzm ,.false.) ! false signals: do not broadcast if ( myid == root_t) then call Init( sds, hdf, 'rzzm', (/imr,jmr,lmr,ntracet/), 'real(4)', status ) call SetDim(sds,0,'LON'//region_name(region),'deg',(/(xbeg(region)+(i+0.5)*dx/xref(region),i=0,imr-1)/),status) call SetDim(sds,1,'LAT'//region_name(region),'deg',(/(ybeg(region)+(i+0.5)*dy/yref(region),i=0,jmr-1)/),status) call SetDim(sds,2,'HYBRID','number',(/(i,i=1,lmr)/),status) call SetDim(sds,3,'ntracet','None',(/(i,i=1,ntracet)/),status) call WriteData(sds,fieldglob4d(1:imr,1:jmr,1:lmr,1:ntracet),status) call Done(sds,status) IF_ERROR_RETURN(status=1) endif #endif #endif ! clear deallocate(fieldglob4d) nullify(m) nullify(rm) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #ifdef secmom nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) #endif #endif ! ! short lived tracers ! ! only if there are non-transported tracers if ( ntrace_chem > 0 ) then ! ahoc ... if ( root_t /= root_k ) then write (gol,'("roots for t and k not equal: ",2i4)') root_t, root_k; call goErr write (gol,'("in ",a)') rname; status=1; call goErr endif ! allocate global array on all PEs allocate(fieldglob4d(imr,jmr,lmr,ntrace_chem)) ! tracers parallel over levels: rmc => chem_dat(region)%rm_k ! collect on root_k call Par_Gather_Tracer_k(fieldglob4d,imr,jmr,lmr,0,0,0,ntrace_chem,rmc,.false.,status) IF_NOTOK_RETURN(status=1) if ( myid == root_k ) then ! write tracers: call io_write( hdf%id, imr,'LON'//region_name(region),& jmr,'LAT'//region_name(region), & lmr,'HYBRID',ntrace_chem,'NTRACE_CHEM',fieldglob4d,'rmc') end if call Par_Barrier ! clear nullify(rmc) deallocate(fieldglob4d) end if ! ntrace_chem > 0 ! ! close ! if ( myid == root_t ) then call Done( hdf, status ) IF_NOTOK_RETURN(status=1) end if call Par_Barrier end do !regions... call DoneRc( rcF, status ) IF_NOTOK_RETURN(status=1) call tstamp(kmain,itau,msg) if (cdebug) then call tstamp(kdebug,itau,'savehdf') end if ! ok status = 0 end subroutine SaveHDF end module io_save