!### 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 #define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if ! #include "tm5.inc" #include "tmm.inc" ! !----------------------------------------------------------------------------- ! TM5-MP ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: IO_SAVE ! ! !DESCRIPTION: Contains routines to read HDF emissions files and ! the main model state (as replacement for restart file). !\\ !\\ ! !INTERFACE: ! MODULE IO_SAVE ! ! !USES: ! use GO, only : gol, goErr, goPr, goBug use GO, only : TrcFile, InitRc => Init, DoneRc => Done, ReadRc use global_data, only : rcfile IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! PUBLIC :: READHDFMMR ! Read mmix from "old save file" (istart=4 ) PUBLIC :: READ_MMIX ! Read mmix file (mean mixing ratio of transported tracers) (istart=5 ) PUBLIC :: READ_SAVE_FILE_30 ! Read save-file (mixing ratio, opt. slopes, of transported tracers) (istart=30) PUBLIC :: READ_SAVE_FILE ! Read save file (mass of all tracers, grid info) (istart=31) PUBLIC :: WRITE_SAVE_FILE ! Write save-file (to be read with read_save_file) PUBLIC :: READ_HDF4 ! Read a record from hdf file (new way) PUBLIC :: READTM3HDF ! Read a record from hdf file (old way, a la TM3) PUBLIC :: READ_SCATTER_HDF ! Read a record from hdf file, and do the MPI scattering ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'io_save' type(TrcFile) :: rcF integer :: fid ! file id for IF_NOTOK_MDF macro ! ! !INTERFACE: ! INTERFACE READ_HDF4 MODULE PROCEDURE READ_HDF4_2DR END INTERFACE ! ! !REVISION HISTORY: ! 16 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) The READTM3 routines are kept for quick fix when converting old code. ! All new code should use the MDF module to read any data. ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_HDF4_2DR ! ! !DESCRIPTION: read one 2D record of real from an HDF4 file. ! !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_HDF4_2DR(fname, odata, status, dir, vname, vid) ! ! !USES: ! USE MDF, ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close USE global_data, ONLY : inputdir ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: fname ! name of HDF file to read character(len=*), intent(in), optional :: dir ! dir where file is. Default is "inputdir" character(len=*), intent(in), optional :: vname ! record name. Default is to use record #0. integer, intent(in), optional :: vid ! record id #. Default is to use record #0. ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: odata(:,:) integer, intent(out) :: status ! ! !REVISION HISTORY: ! 28 Mar 2012 - P. Le Sager - v0 ! ! !REMARKS: ! (1) This wrapper around MDF routines is only for convenience: a lot of ! calls for the chemistry emissions, and participate to the effort ! of getting rid of file_hdf.f90 routine, in favor of the MDF module. ! (2) By default **FIRST** record is read, unless another variable number ! (vid) and/or variable name (vname) is passed. If both vid and vname ! are passed, vname is ignored. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/read_hdf4_2dr' character(len=512):: Lvname, Ldir integer :: Lvid !--- check input Ldir=inputdir if(present(dir)) Ldir=dir Lvname='' if(present(vname)) Lvname=vname Lvid=1 if (present(vid)) then Lvid=vid Lvname='' end if !--- read CALL MDF_Open( trim(Ldir)//'/'//TRIM(fname), MDF_HDF4, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) if (len(trim(Lvname))/=0) then CALL MDF_Inq_VarID( fid, TRIM(Lvname), Lvid, status ) IF_NOTOK_MDF(status=1) end if CALL MDF_Get_Var( fid, Lvid, odata, status ) IF_NOTOK_MDF(status=1) CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) status=0 END SUBROUTINE READ_HDF4_2DR !EOC ! ====================================================================== !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 ) #ifdef with_hdf4 use file_hdf, only : THdfFile, Init, Done use io_hdf , only : FAIL use io_hdf , only : io_read2d_32g, io_read3d_32g #endif 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 ---------------------------------------- #ifdef with_hdf4 type(THdfFile) :: hdf #endif integer :: status integer :: idx real, allocatable :: field1d(:) ! first dim will be 1 real, allocatable :: field2d(:,:) real, allocatable :: field3d(:,:,:) ! --- begin ------------------------------------------- #ifdef with_hdf4 ! open hdf file: call Init( hdf, trim(filename), 'read', status ) if (status/=0) call escape_tm('ERROR in readtm3hdf') #else call escape_tm('not compiled with hdf4') #endif !FD hdf files indices run from 0 ...n-1 idx = month select case(rank) case(1) ! latitudional field allocate(field1d(jm_emis)) #ifdef with_hdf4 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 write(gol,*) 'readtm3hdf: io_read2d_32g failed:' ; call goPr write(gol,*) 'readtm3hdf: filename : ', trim(filename) ; call goPr write(gol,*) 'readtm3hdf: rank : ', rank ; call goPr write(gol,*) 'readtm3hdf: im_emis : ', im_emis ; call goPr write(gol,*) 'readtm3hdf: jm_emis : ', jm_emis ; call goPr write(gol,*) 'readtm3hdf: lm_emis : ', lm_emis ; call goPr write(gol,*) 'readtm3hdf: month : ', month ; call goPr write(gol,*) 'readtm3hdf: field_name : ', trim(field_name); call goPr if (present(idt)) write(gol,*) 'readtm3hdf: idt : ', idt; call goPr call escape_tm('readtm3hdf: Error reading HDF file 1D') end if #else call escape_tm('not compiled with hdf4') #endif field(1,:,1)=field1d deallocate(field1d) case(2) ! 2D surface field allocate(field2d(im_emis,jm_emis)) #ifdef with_hdf4 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 write(gol,*) 'readtm3hdf: io_read2d_32g failed:' ; call goPr write(gol,*) 'readtm3hdf: filename : ', trim(filename) ; call goPr write(gol,*) 'readtm3hdf: rank : ', rank ; call goPr write(gol,*) 'readtm3hdf: im_emis : ', im_emis ; call goPr write(gol,*) 'readtm3hdf: jm_emis : ', jm_emis ; call goPr write(gol,*) 'readtm3hdf: lm_emis : ', lm_emis ; call goPr write(gol,*) 'readtm3hdf: month : ', month ; call goPr write(gol,*) 'readtm3hdf: field_name : ', trim(field_name) ; call goPr if (present(idt)) write(gol,*) 'readtm3hdf: idt : ', idt ; call goPr call escape_tm('readtm3hdf: Error reading HDF file 2D') endif #else call escape_tm('not compiled with hdf4') #endif field(:,:,1)=field2d deallocate(field2d) case (3) ! 3D field allocate(field3d(im_emis,jm_emis,lm_emis)) #ifdef with_hdf4 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 write(gol,*) 'readtm3hdf: io_read2d_32g failed:' ; call goPr write(gol,*) 'readtm3hdf: filename : ', trim(filename) ; call goPr write(gol,*) 'readtm3hdf: rank : ', rank ; call goPr write(gol,*) 'readtm3hdf: im_emis : ', im_emis ; call goPr write(gol,*) 'readtm3hdf: jm_emis : ', jm_emis ; call goPr write(gol,*) 'readtm3hdf: lm_emis : ', lm_emis ; call goPr write(gol,*) 'readtm3hdf: month : ', month ; call goPr write(gol,*) 'readtm3hdf: field_name : ', trim(field_name) ; call goPr if (present(idt)) write(gol,*) 'readtm3hdf: idt : ', idt ; call goPr call escape_tm('readtm3hdf: Error reading HDF file 3D') endif #else call escape_tm('not compiled with hdf4') #endif field(:,:,:)=field3d deallocate(field3d) case (23) ! 2D field lat-pres allocate(field2d(jm_emis,lm_emis)) #ifdef with_hdf4 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! write(gol,*) '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 write(gol,*) 'readtm3hdf: io_read2d_32g failed:' ; call goPr write(gol,*) 'readtm3hdf: filename : ', trim(filename) ; call goPr write(gol,*) 'readtm3hdf: rank : ', rank ; call goPr write(gol,*) 'readtm3hdf: im_emis : ', im_emis ; call goPr write(gol,*) 'readtm3hdf: jm_emis : ', jm_emis ; call goPr write(gol,*) 'readtm3hdf: lm_emis : ', lm_emis ; call goPr write(gol,*) 'readtm3hdf: month : ', month ; call goPr write(gol,*) 'readtm3hdf: field_name : ', trim(field_name) ; call goPr if (present(idt)) write(gol,*) 'readtm3hdf: idt : ', idt ; call goPr call escape_tm('readtm3hdf: Error reading HDF file 23') end if field(:,:,:)=field3d deallocate(field3d) else field(1,:,:)=field2d end if #else call escape_tm('not compiled with hdf4') #endif 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 #ifdef with_hdf4 ! close: call Done( hdf, status ) if (status/=0) call escape_tm('ERROR in readtm3hdf') #else call escape_tm('not compiled with hdf4') #endif !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 !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_SCATTER_HDF ! ! !DESCRIPTION: simple adaptation of readtm3hdf to use MDF module, and ! scatter data. !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_SCATTER_HDF( filename, rank, region, i1, j1, field, status, & idx, field_name, idt ) ! ! !USES: ! use partools, only : isRoot use tm5_distgrid, only : dgrid, Get_DistGrid, scatter ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: filename ! HDF4 file to read integer, intent(in) :: rank ! 1,2,.. [only 2 implemented for now] integer, intent(in) :: region, i1, j1 ! region and local array start index integer, intent(in), optional :: idx ! variable id in [1..nmax] character(len=*), intent(in), optional :: field_name ! variable name, ignored if idx set integer, intent(in), optional :: idt(6) ! date [not implemented] ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: field(i1:,j1:) ! scattered read data integer, intent(out) :: status ! ! !REVISION HISTORY: ! 23 Jan 2012 - P. Le Sager - for lat-lon mpi decomposition ! ! !REMARKS: ! (1) tested only with rank=2 (see dry_deposition.F90) ! (2) not all features of readtm3hdf are available, but probably also ! not needed... will add as needed. !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/read_scatter_hdf' integer :: im_emis, jm_emis real, allocatable :: field1d(:) real, allocatable :: field2d(:,:) ! --- begin ------------------------------------------- CALL GET_DISTGRID( DGRID(region), I_WRLD=im_emis, J_WRLD=jm_emis ) status=0 READ: if (isRoot) then allocate(field2d(im_emis,jm_emis)) field2D=0. 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 ! write(gol,*) 'readtm3hdf: io_read2d_32g failed:' ! write(gol,*) 'readtm3hdf: filename : ', trim(filename) ! write(gol,*) 'readtm3hdf: rank : ', rank ! write(gol,*) 'readtm3hdf: im_emis : ', im_emis ! write(gol,*) 'readtm3hdf: jm_emis : ', jm_emis ! write(gol,*) 'readtm3hdf: month : ', month ! write(gol,*) 'readtm3hdf: field_name : ', trim(field_name) ! if (present(idt)) write(gol,*) 'readtm3hdf: idt : ', idt ! call escape_tm('readtm3hdf: Error reading HDF file 1D') ! end if ! field2D(1,:)=field1d ! deallocate(field1d) case(2) ! 2D surface field call READ_HDF4(filename, field2d, status, dir='', vname=field_name, vid=idx) IF_NOTOK_RETURN(status=1) case default write (gol,*) rname//' : Please check the rank of the required field,', & ' must be 1,2,3, or 23'; call goPr status=1 end select else allocate(field2d(1,1)) end if READ call SCATTER( dgrid(region), field, field2d, 0, status) IF_NOTOK_RETURN(status=1) deallocate(field2d) END SUBROUTINE READ_SCATTER_HDF !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_SAVE_FILE_30 ! ! !DESCRIPTION: reads initial MIXING RATIO (+ slopes, and second moments if ! they are switched on) of TRANSPORTED tracers from HDF-4 ! "save files". The save file must be specified in rc file with key: ! ! start.30. ! ! No regridding available. ! ! Used for istart=30. !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_SAVE_FILE_30( region, file_name, status ) ! ! !USES: ! #ifdef with_hdf4 use io_hdf #endif use dims, only : im, jm, lm use global_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : ntracet use partools, only : isRoot use tm5_distgrid, only : dgrid, Get_DistGrid, scatter ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region character(len=*), intent(in) :: file_name ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 4 Oct 2011 - P. Le Sager - TM5-MP version ! ! !REMARKS: ! (1) no zoom (must put from_file stuff back in) ! (2) FIXME: not tested... ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/read_save_file_30' 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 #ifdef with_hdf4 integer :: io, sfstart integer :: istat integer :: sfend #endif integer :: i,j,l,n integer :: imr,jmr,lmr, halo real,dimension(:,:,:,:), allocatable :: field real,dimension(:,:,:), allocatable :: msave, msave_local integer :: i1, i2, j1, j2 ! --- begin -------------------------------- m => m_dat(region)%data rm => mass_dat(region)%rm #ifdef slopes rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm #ifdef secmom rxxm => mass_dat(region)%rxxm rxym => mass_dat(region)%rxym rxzm => mass_dat(region)%rxzm ryym => mass_dat(region)%ryym ryzm => mass_dat(region)%ryzm rzzm => mass_dat(region)%rzzm #endif #endif halo = mass_dat(region)%halo ! global array imr = im(region) ; jmr = jm(region) ; lmr = lm(region) if ( isRoot ) then allocate(field(imr, jmr, lmr, ntracet)) field = 0.0 allocate(msave(imr,jmr,lmr)) else allocate(field(1,1,1,1)) allocate(msave(1,1,1)) end if ! local indices call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate(msave_local(i1:i2,j1:j2,lmr)) #ifdef with_hdf4 ! read on root if ( isRoot ) then io = sfstart( file_name, DFACC_READ ) if ( io == FAIL ) then write(gol,*) 'Read_save_file_30 : Could not open file '// trim(file_name) ; call goPr call goErr TRACEBACK ; status=1 ; return end if write(gol,'("Read_save_file_30: opened for READ access: ",a)') trim(file_name); call goPr call io_read3d_32g(io,imr,jmr,lmr,msave,'m',status) IF_NOTOK_RETURN(status=1) call io_read4d_32g(io,imr,jmr,lmr,ntracet, field,'rm',status) ! try again if( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: Read rm, status = ",i6)') status) end if !root ! get local chunk call SCATTER( dgrid(region), rm, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 call SCATTER( dgrid(region), msave_local, msave, 0, status) IF_NOTOK_RETURN(status=1) #ifdef slopes if ( isRoot ) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rxm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) end if ! get local chunk call SCATTER( dgrid(region), rxm, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rym',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rym',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) end if !root ! get local chunk call SCATTER( dgrid(region), rym, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if ( isRoot ) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rzm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rzm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) end if ! get local chunk call SCATTER( dgrid(region), rzm, field, halo, status) IF_NOTOK_RETURN(status=1) ! reset field = 0.0 #endif ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! >>> second moments ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #ifdef secmom if ( isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxxm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rxxm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), rxxm, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxym',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rxym',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), rxym, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rxzm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rxzm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), rxzm, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'ryym',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'ryym',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), ryym, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'ryzm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'ryzm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), ryzm, field, halo, status) IF_NOTOK_RETURN(status=1) ! reset field = 0.0 if(isRoot) then call io_read4d_32g(io,imr,jmr,lmr,ntracet,field,'rzzm',status) if ( status /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,imr,jmr,lmr,field,'rzzm',status) end if IF_NOTOK_RETURN(write (*,'("read_save_file_30: status = ",i6)') status) endif call SCATTER( dgrid(region), rzzm, field, halo, status) IF_NOTOK_RETURN(status=1) field = 0.0 #endif ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! <<< end second moments ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< if ( isRoot ) then istat = sfend(io) if ( istat /= FAIL ) then write(gol,'("read_save_file_30: ... file closed")') ; call goPr write(gol,'("read_save_file_30: ")') ; call goPr else TRACEBACK ; status=1 ; return end if end if ! CONVERSION ! rm/msave ---> mixing ratio on file. mr*m --> something back in kilo's if ( isRoot ) then write(gol,*) 'read_save_file_30: Maximum ratio of the air masses:', & maxval(abs(msave/m(1:imr,1:jmr,1:lmr))) ; call goPr write(gol,*) 'read_save_file_30: Minimum ratio of the air masses:', & minval(abs(msave/m(1:imr,1:jmr,1:lmr))) ; call goPr end if do n = 1, ntracet do l = 1, lmr rm(i1:i2,j1:j2,l,n) = rm (i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) #ifdef slopes rxm(i1:i2,j1:j2,l,n) = rxm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) rym(i1:i2,j1:j2,l,n) = rym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) rzm(i1:i2,j1:j2,l,n) = rzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) #ifdef secmom rxxm(i1:i2,j1:j2,l,n) = rxxm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) rxym(i1:i2,j1:j2,l,n) = rxym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) rxzm(i1:i2,j1:j2,l,n) = rxzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) ryym(i1:i2,j1:j2,l,n) = ryym(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) ryzm(i1:i2,j1:j2,l,n) = ryzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) rzzm(i1:i2,j1:j2,l,n) = rzzm(i1:i2,j1:j2,l,n)*m(i1:i2,j1:j2,l)/msave_local(i1:i2,j1:j2,l) #endif #endif end do end do #else write (gol,'("not compiled with hdf4")'); call goErr TRACEBACK; status=1; return #endif // with_hdf4 ! cleanup deallocate(msave, msave_local) 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 READ_SAVE_FILE_30 !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READHDFMMR ! ! !DESCRIPTION: read tracer mixing ratio from so-called "old save" file and ! uses air mass to convert to mass, in order to initialize the ! model. ! ! Used with option istart=4 !\\ !\\ ! !INTERFACE: ! SUBROUTINE READHDFMMR( region, file_name, status ) ! ! !USES: ! use dims, only : im, jm, lm, nregions, datadir, region_name #ifdef with_hdf4 use io_hdf #endif use global_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : fscale, ntracet use tm5_distgrid, only : dgrid, Get_DistGrid, scatter use partools, only : isRoot ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region character(len=*), intent(in) :: file_name ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 4 Jul 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition (FIXME: not tested) ! ! !REMARKS: ! (1) similar to read_mmix, but all tranported tracers are read at once. ! That is: no check on tracers name, and order of tracer matters, ! which depends on the number of processors in TM5 v3. ! (2) All tracers must be in the file, and the 4D array should be the ! second dataset. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/readhdfmmr' real, dimension(:,:,:,:), pointer :: rm real, dimension(:,:,:), pointer :: m #ifdef with_hdf4 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 #endif 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 ! start m => m_dat(region)%data rm => mass_dat(region)%rm #ifdef slopes mass_dat(region)%rxm = 0. mass_dat(region)%rym = 0. mass_dat(region)%rzm = 0. #ifdef secmom mass_dat(region)%rxxm = 0. mass_dat(region)%rxym = 0. mass_dat(region)%rxzm = 0. mass_dat(region)%ryym = 0. mass_dat(region)%ryzm = 0. mass_dat(region)%rzzm = 0. #endif #endif if (isRoot) then allocate(field(im(region),jm(region),lm(region),ntracet)) else allocate(field(1,1,1,1)) end if #ifdef with_hdf4 ! Read on root if ( isRoot ) then io = sfstart( file_name, DFACC_READ ) if ( io == -1 ) then write(gol,*)'readhdfmmr: Could not open file:'//file_name ; call goErr TRACEBACK; status=1; return end if ind = 1 ! select 1st dataset....(0 = m) istat = sffinfo(io, n_datasets, n_file_attrs) write(gol,*) 'readhdfmmr: Number of datasets and attributes', & n_datasets, n_file_attrs ; call goPr call io_read4d_32g(io,im(region),jm(region),lm(region),ntracet, & field,'rm',ifail,index=ind) if ( ifail /= 0 .and. ntracet == 1 ) then call io_read3d_32g(io,im(region),jm(region),lm(region), & field,'rm',ifail,index=ind) end if if ( ifail /= 0 ) then write(gol,*) 'readhdfmmr: '//trim(file_name) ; call goErr write(gol,*)'readhdfmmr: Failed to read fields'; call goErr TRACEBACK; status=1; return end if write(gol,*) 'readhdfmmr: Read rm, ifail = ', ifail istat = sfend(io) if (istat /= FAIL) then write(gol,*)'readhdfmmr: ... file closed' ; call goErr write(gol,*)' ' ; call goPr else write(gol,*)'readhdfmmr: ERROR in restart from HDF file'; call goErr TRACEBACK; status=1; return end if end if ! root #else write (gol,'("not compiled with hdf4")'); call goErr TRACEBACK; status=1; return #endif // with_hdf4 call SCATTER( dgrid(region), rm, field, mass_dat(region)%halo, status) do n=1,ntracet rm(:,:,:,n) = rm(:,:,:,n)*m(:,:,:)/fscale(n) end do deallocate(field) nullify(m) nullify(rm) END SUBROUTINE READHDFMMR !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_MMIX ! ! !DESCRIPTION: read mean mixing ratio of TRANSPORTED tracers from an "mmix hdf ! file" and uses air mass to convert to mass. Missing tracers ! are set to 0.0 ! ! Used with option istart=5. !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_MMIX( region, file_name, status ) ! ! !USES: ! use dims, only : im, jm, lm, nregions, datadir, region_name, parent #ifdef with_hdf4 use io_hdf #endif use global_data, only : mass_dat use meteodata , only : m_dat use chem_param, only : fscale, ntracet, names use tm5_distgrid, only : dgrid, Get_DistGrid, scatter use ParTools, only : isRoot ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region character(len=*), intent(in) :: file_name ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 4 Jul 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition (FIXME: not tested) ! ! !REMARKS: ! ! (1) This is typically the choice for combining different versions ! or extending the number of tracers. ! (2) The compounds are searched by name. If not in the mmix file ! the field is initialized as zero. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/read_mmix' real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm,rym,rzm #endif real,dimension(:,:,:),pointer :: m #ifdef with_hdf4 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 #endif 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,from_file, my_parent real,dimension(:,:,:,:), allocatable :: field real,dimension(:,:,:), allocatable :: field3d ! start m => m_dat(region)%data rm => mass_dat(region)%rm #ifdef slopes rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm #endif rm = 0.0 #ifdef slopes rxm = 0.0 rym = 0.0 rzm = 0.0 #endif if (isRoot) then allocate(field(im(region),jm(region),lm(region),ntracet)) else allocate(field(1,1,1,1)) end if #ifdef with_hdf4 ! Read on root if ( isRoot ) 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 write(gol,*)'read_mmix : Could not open file and no parent '// & 'available :'//file_name ; call goErr TRACEBACK; status=1; return else ! FIXME : this is not going anywhere from_file = 0 write(gol,*) 'readhdf: Trying to initialise '//trim(region_name(region))// & ' from parent....'; call goPr endif endif if ( from_file == 1 ) then ! read from file: write(gol,*) 'read_mmix: ',file_name,'... opened for READ access.' ; call goPr istat = sffinfo(io, n_datasets, n_file_attrs) write(gol,*) 'read_mmix: Number of datasets and attributes', n_datasets, n_file_attrs call goPr 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 write(gol,*) 'Read_mmix:', names(n), ' not found in dataset --> 0.0 ' ; call goPr field(1:im(region), 1:jm(region), 1:lm(region), n) = tiny(1.0) endif enddo istat = sfend(io) if (istat /= FAIL) then write(gol,*)'read_mmix: ... file closed' ; call goPr write(gol,*)' '; call goPr else write(gol,*)'read_mmix: ERROR in restart from HDF file'; call goErr TRACEBACK; status=1; return end if deallocate(field3d) ! no longer needed endif ! from file endif ! root #else write (gol,'("not compiled with hdf4")'); call goErr TRACEBACK; status=1; return #endif // with_hdf4 call SCATTER( dgrid(region), rm, field, mass_dat(region)%halo, status) ! Convert from mmix to kg tracer do n = 1, ntracet rm(:,:,:,n) = rm(:,:,:,n)*m/fscale(n) end do #ifdef slopes rxm = 0.0 ; rym = 0.0 ; rzm = 0.0 #endif ! Done deallocate(field) nullify(m, rm) #ifdef slopes nullify(rxm, rym, rzm) #endif END SUBROUTINE READ_MMIX !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: READ_SAVE_FILE ! ! !DESCRIPTION: Read mass of both transpoted AND short lived species from so ! -called "save file". No slopes, but regridding to model ! resolution is available, since grid information is also read. ! ! Used with istart=31 option. !\\ !\\ ! !INTERFACE: ! SUBROUTINE READ_SAVE_FILE( region, fname, status, tm4 ) ! ! !USES: ! #ifdef with_hdf4 use file_hdf , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData #endif use grid , only : TllGridInfo, TLevelInfo, Init, Done, Fill3D use dims , only : im, jm, lm use chem_param , only : ntrace, ntracet, fscale, names, ntrace_chem use partools , only : root, ierr, my_real, localComm, isRoot use global_data , only : mass_dat, chem_dat use meteodata , only : global_lli, levi, sp_dat, m_dat use tm5_distgrid, only : dgrid, Get_DistGrid, gather, scatter ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region character(len=*), intent(in) :: fname logical, intent(in), optional :: tm4 ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 4 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) tracers not found in the file are left unchanged (ie to their initialized values) ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Read_save_file' ! 31 (from 60) layer defintition for TM4 save files: real, parameter :: at_tm4_60_31(0:31) = (/ & 0., 7.367743, 467.333588, 1385.912598, 2887.696533, & 4941.778320, 6144.314941, 8802.356445, 10209.500977, 11632.758789, & 14411.124023, 15706.447266, 16899.468750, 17961.357422, 18864.750000, & 19584.330078, 20097.402344, 20384.480469, 20222.205078, 19755.109375, & 19027.695313, 18045.183594, 16819.474609, 15379.805664, 13775.325195, & 12077.446289, 10376.126953, 8765.053711, 6018.019531, 3960.291504, & 2082.273926, 0./) real, parameter :: bt_tm4_60_31(0:31) = (/ & 1.00000000, 0.99401945, 0.96764523, 0.93194032, 0.87965691, & 0.81125343, 0.77159661, 0.68326861, 0.63554746, 0.58616841, & 0.48477158, 0.43396294, 0.38389215, 0.33515489, 0.28832296, & 0.24393314, 0.20247590, 0.16438432, 0.09967469, 0.07353383, & 0.05169041, 0.03412116, 0.02067788, 0.01114291, 0.00508112, & 0.00181516, 0.00046139, 0.00007582, 0.00000000, 0.00000000, & 0.00000000, 0.00000000/) ! 34 (from 91) layer defintition for TM4 save files: real, parameter :: at_tm4_91_34(0:34) = (/ & 0.000000, 6.575628, 336.772369, 1297.656128, 3010.146973, & 5422.802734, 8356.252930, 11543.166992, 14665.645508, 17385.595703, & 19348.775391, 20319.011719, 20348.916016, 19919.796875, 19184.544922, & 18191.029297, 16990.623047, 15638.053711, 14192.009766, 12713.897461, & 11262.484375, 9873.560547, 8564.624023, 7341.469727, 6199.839355, & 4663.776367, 3358.425781, 2292.155518, 1463.163940, 857.945801, & 450.685791, 204.637451, 76.167656, 21.413612, 0.000000/) real, parameter :: bt_tm4_91_34(0:34) = (/ & 1.000000, 0.994204, 0.973466, 0.935157, 0.875518, & 0.795385, 0.698224, 0.589317, 0.475016, 0.362203, & 0.259554, 0.176091, 0.112979, 0.080777, 0.055474, & 0.036227, 0.022189, 0.012508, 0.006322, 0.002765, & 0.001000, 0.000279, 0.000055, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, & 0.000000, 0.000000, 0.000000, 0.000000, 0.000000/) ! 34 (from standard) layer defintition for TM4 save files: real, parameter :: at_tm4_std_34(0:34) = (/ & 0., 7.367743, 210.393890, 855.361755, & 2063.779785, 3850.913330, 6144.314941, 8802.356445, & 11632.758789, 14411.124023, 16899.468750, 17961.357422, & 18864.750000, 19584.330078, 20097.402344, 20384.480469, & 20429.863281, 20222.205078, 19755.109375, 19027.695313, & 18045.183594, 16819.474609, 15379.805664, 13775.325195, & 12077.446289, 10376.126953, 8765.053711, 7306.631348, & 6018.019531, 3960.291504, 1680.640259, 713.218079, & 298.495789, 95.636963, 0./) real, parameter :: bt_tm4_std_34(0:34) = (/ & 1.00000000, 0.99401945, 0.97966272, 0.95182151, & 0.90788388, 0.84737492, 0.77159661, 0.68326861, & 0.58616841, 0.48477158, 0.38389215, 0.33515489, & 0.28832296, 0.24393314, 0.20247590, 0.16438432, & 0.13002251, 0.09967469, 0.07353383, 0.05169041, & 0.03412116, 0.02067788, 0.01114291, 0.00508112, & 0.00181516, 0.00046139, 0.00007582, 0.00000000, & 0.00000000, 0.00000000, 0.00000000, 0.00000000, & 0.00000000, 0.00000000, 0.00000000/) ! --- local ------------------------------- integer :: imr, jmr, lmr real, dimension(:,:,:), pointer :: sp real, dimension(:,:,:), pointer :: m #ifdef with_hdf4 type(THdfFile) :: hdf type(TSds) :: sds_rm, sds_rmc #endif integer :: n, ns integer :: ntrace_in, ntracet_in integer :: n_in, nl character(len=8), allocatable :: names_in(:) integer :: xbeg_in, ybeg_in real :: dx_in, dy_in integer :: im_in, jm_in, lm_in real, allocatable :: at_in(:), bt_in(:) type(TllGridInfo) :: lli_in type(TLevelInfo) :: levi_in real, allocatable :: rm_in(:,:,:,:), sp_gbl(:,:,:) real, allocatable :: rmt(:,:,:,:),rms(:,:,:,:) integer :: l logical :: with_tm4 ! --- begin ------------------------------- ! global grid sizes imr = im(region) jmr = jm(region) lmr = lm(region) ! set pointers sp => sp_dat(region)%data ! surface pressure m => m_dat(region)%data ! air mass ! special input ? with_tm4 = .false. if ( present(tm4) ) with_tm4 = tm4 ! temporary global storage if (isRoot) then allocate( rmt(imr,jmr,lmr,ntracet ) ) allocate( rms(imr,jmr,lmr,ntrace_chem) ) allocate( sp_gbl(imr,jmr,1) ) else allocate( rmt(1,1,1,1) ) allocate( rms(1,1,1,1) ) allocate(sp_gbl(1,1,1)) endif ! init to 0 in case of not found in file rmt=tiny(0.) rms=tiny(0.) call gather( dgrid(region), sp, sp_gbl, sp_dat(region)%halo, status) ! get global surf pressure IF_NOTOK_RETURN(status=1) #ifdef with_hdf4 ! read on root (until we implement the stride in reading) if ( isRoot ) then ! open hdf file call Init( hdf, fname, 'read', status ) IF_NOTOK_RETURN(status=1) ! read list with tracer names in file: call ReadAttribute( hdf, 'ntrace', ntrace_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'ntracet', ntracet_in, status ) IF_NOTOK_RETURN(status=1) allocate( names_in(ntrace_in) ) call ReadAttribute( hdf, 'names', names_in, status ) IF_NOTOK_RETURN(status=1) ! read level definition call ReadAttribute( hdf, 'lm', lm_in, status ) IF_NOTOK_RETURN(status=1) allocate( at_in(0:lm_in) ) allocate( bt_in(0:lm_in) ) if ( with_tm4 ) then select case ( lm_in ) case ( 31 ) at_in = at_tm4_60_31 bt_in = bt_tm4_60_31 case ( 34 ) at_in = at_tm4_std_34 bt_in = bt_tm4_std_34 case default write (gol,'("input from TM4 save file implemented for 60-31 and 91-34levels only ...")'); call goErr TRACEBACK; status=1; return end select else call ReadAttribute( hdf, 'at', at_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'bt', bt_in, status ) IF_NOTOK_RETURN(status=1) end if call Init( levi_in, lm_in, at_in, bt_in, status ) IF_NOTOK_RETURN(status=1) deallocate( at_in ) deallocate( bt_in ) ! read grid definition: call ReadAttribute( hdf, 'im', im_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'jm', jm_in, status ) IF_NOTOK_RETURN(status=1) if ( with_tm4 ) then xbeg_in = -180.0 dx_in = 360.0/im_in ybeg_in = -90.0 dy_in = 180.0/jm_in else call ReadAttribute( hdf, 'xbeg', xbeg_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'dx', dx_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'ybeg', ybeg_in, status ) IF_NOTOK_RETURN(status=1) call ReadAttribute( hdf, 'dy', dy_in, status ) IF_NOTOK_RETURN(status=1) end if call Init( lli_in, xbeg_in+0.5*dx_in, dx_in, im_in, & ybeg_in+0.5*dy_in, dy_in, jm_in, status ) IF_NOTOK_RETURN(status=1) ! temporary storage: allocate( rm_in(im_in,jm_in,lm_in,1) ) ! open tracer masses: call Init( sds_rm, hdf, 'rm', status ) IF_NOTOK_RETURN(status=1) if ( .not. with_tm4 ) then call Init( sds_rmc, hdf, 'rmc', status ) if(status /= 0) write (gol,'("WARNING - No short lived tracer array found in save file...continuing execution")') call goPr end if ! loop over all tracers: do n = 1, ntrace ! search input record: n_in = -1 do nl = 1, ntrace_in if ( names(n) == names_in(nl) ) then n_in = nl exit end if end do ! not found ? if ( n_in < 0 ) then if ( n <= ntracet ) then write (gol,'("WARNING - transported tracer not found in save file:")') else write (gol,'("WARNING - short-lived tracer not found in save file:")') end if call goPr write (gol,'("WARNING - file : ",a)') trim(fname) ; call goPr write (gol,'("WARNING - tracer : ",i3," ",a)') n, trim(names(n)) ; call goPr if ( n <= ntracet ) then write (gol,'(" initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rmt(:,:,:,n)) else ns = n-ntracet write (gol,'(" initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rms(:,:,:,ns)) endif call goPr cycle end if ! transported or short lived ? if ( n <= ntracet ) then ! read 3D tracer field from rm or rmc: if ( n_in <= ntracet_in ) then call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) else call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) end if ! expand polar cell for TM4 files ... if ( with_tm4 ) then do l = 1, lm_in rm_in(:, 1,l,1) = rm_in(1, 1,l,1) / im_in rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in end do end if ! regrid from saved grid to current model grid call Fill3D( global_lli(region), levi, 'n', sp_gbl(1:imr,1:jmr,1), rmt(1:imr,1:jmr,1:lmr,n), & lli_in, levi_in, rm_in(:,:,:,1), 'sum', status ) IF_NOTOK_RETURN(status=1) ! info write (gol,'(" initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rmt(:,:,:,n)) call goPr !Need to gather "m" globally to get VMR ! maxval(rmt(1:imr,1:jmr,1:lmr,n)/m_gbl(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr else ! short lived tracer ! read 3D tracer field: if ( with_tm4 ) then ! tm4 save files have short lived tracers in rm call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) ! expand polar cell ... do l = 1, lm_in rm_in(:, 1,l,1) = rm_in(1, 1,l,1) / im_in rm_in(:,jm_in,l,1) = rm_in(1,jm_in,l,1) / im_in end do else ! read 3D tracer field from rm or rmc: if ( n_in <= ntracet_in ) then call ReadData( sds_rm, rm_in, status, start=(/0,0,0,n_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) else call ReadData( sds_rmc, rm_in, status, start=(/0,0,0,n_in-ntracet_in-1/) ) ! kg IF_NOTOK_RETURN(status=1) end if end if ! regrid from saved grid to current model grid ns = n-ntracet call Fill3D( global_lli(region), levi, 'n', sp_gbl(1:imr,1:jmr,1), rms(1:imr,1:jmr,1:lmr,ns), & lli_in, levi_in, rm_in(:,:,:,1), 'sum', status ) IF_NOTOK_RETURN(status=1) ! info write (gol,'(" initialised tracer ",i2,x,a,"; max mass : ",e10.3)') n, names(n), maxval(rms(1:imr,1:jmr,1:lmr,ns)) call goPr !Need to gather "m" globally to get VMR ! maxval(rms(1:imr,1:jmr,1:lmr,ns)/m(1:imr,1:jmr,1:lmr)*fscale(n)); call goPr end if end do ! cleanup deallocate( names_in ) deallocate( rm_in ) ! done with grid call Done( lli_in, status ) IF_NOTOK_RETURN(status=1) call Done( levi_in, status ) IF_NOTOK_RETURN(status=1) ! close data sets call Done( sds_rm, status ) IF_NOTOK_RETURN(status=1) if ( .not. with_tm4 ) then call Done( sds_rmc, status ) end if ! close hdf file (all pe) call Done( hdf, status ) IF_NOTOK_RETURN(status=1) end if #else write (gol,'("not compiled with hdf4")'); call goErr TRACEBACK; status=1; return #endif ! distribute call scatter( dgrid(region), mass_dat(region)%rm, rmt, mass_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) call scatter( dgrid(region), chem_dat(region)%rm, rms, chem_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ! done deallocate( rmt, rms, sp_gbl ) status = 0 END SUBROUTINE READ_SAVE_FILE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: WRITE_SAVE_FILE ! ! !DESCRIPTION: save all essential model parameters and fields in a HDF ! file. This so-called "save file" is created only if no ! restart is written at the end of the run. It can be read ! (with READ_SAVE_FILE) to initialize a run, by using ! istart=31 option. !\\ !\\ ! !INTERFACE: ! SUBROUTINE WRITE_SAVE_FILE( msg, file_name, status ) ! ! !USES: ! #ifdef with_hdf4 use file_hdf , only : THdfFile, Init, Done, WriteAttribute, TSds, SetDim, WriteData #endif 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 : limits use dims , only : at, bt use dims , only : nsplitsteps, splitorder use global_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 #ifdef with_hdf4 use io_hdf, only : io_write #endif use datetime, only : tstamp use global_data, only : outdir use tm5_distgrid, only : dgrid, Get_DistGrid, GATHER use ParTools, only : isRoot ! ! !INPUT PARAMETERS: ! character(len=*), intent(in) :: file_name ! base filename if not found in rc file character(len=*), intent(in) :: msg ! message attribute ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! return code ! ! !REVISION HISTORY: ! ! 11 Nov 2011 - P. Le Sager - adapted for adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Write_save_file' ! --- 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(:,:,:,:), allocatable :: glbfield #ifdef with_hdf4 type(THdfFile) :: hdf type(TSds) :: sds #endif integer :: region, imr, jmr, lmr, i, n, I1, J1, IH1, halo character(len=64) :: msg_file character(len=256) :: fname character(len=32 ) :: key ! --- begin ----------------------------------- call InitRc( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) do region = 1,nregions #ifdef with_hdf4 if ( isRoot ) 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) ! lib cannot handle kind=8: ! 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) ! lib cannot handle kind=8: ! 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 ! Short cuts !---------------------- rm => mass_dat(region)%rm #ifdef slopes rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm #ifdef secmom rxxm => mass_dat(region)%rxxm rxym => mass_dat(region)%rxym rxzm => mass_dat(region)%rxzm ryym => mass_dat(region)%ryym ryzm => mass_dat(region)%ryzm rzzm => mass_dat(region)%rzzm #endif #endif ! Local & global bounds !----------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, J_STRT=j1, I_STRT_HALO=ih1 ) halo = i1-ih1 imr = im(region) ; jmr = jm(region) ; lmr = lm(region) ! Allocate global array !---------------------- if (isRoot) then allocate( glbfield( imr, jmr, lmr, ntracet) ) glbfield = 0.0 else allocate( glbfield(1,1,1,1) ) end if ! M (Air mass) !---------------------- CALL GATHER( dgrid(region), m_dat(region)%data, glbfield(:,:,:,1), m_dat(region)%halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) then call Init( sds, hdf, 'm', (/imr,jmr,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, glbfield(:,:,:,1), status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RM (transported tracer masses) !----------------------------- CALL GATHER( dgrid(region), rm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) then 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if #ifdef slopes ! RXM !----------------------------- CALL GATHER( dgrid(region), rxm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RYM !----------------------------- CALL GATHER( dgrid(region), rym, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RZM !----------------------------- CALL GATHER( dgrid(region), rzm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if #ifdef secmom ! second moments ! RXXM !----------------------------- CALL GATHER( dgrid(region), rxxm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RXYM !----------------------------- CALL GATHER( dgrid(region), rxym, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RXZM !----------------------------- CALL GATHER( dgrid(region), rxzm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot ) then call Init( sds, hdf, 'rxzm', (/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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RYYM !----------------------------- CALL GATHER( dgrid(region), ryym, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) end if ! RYZM !----------------------------- CALL GATHER( dgrid(region), ryzm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) endif ! RZZM !----------------------------- CALL GATHER( dgrid(region), rzzm, glbfield, halo, status ) IF_ERROR_RETURN(status=1) if ( isRoot) 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,glbfield,status) call Done(sds,status) IF_ERROR_RETURN(status=1) endif #endif #endif ! clear deallocate(glbfield) 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 if any ! ---------------------------- if ( ntrace_chem > 0 ) then ! allocate global array if (isRoot) then allocate(glbfield(imr,jmr,lmr,ntrace_chem)) else allocate(glbfield(1,1,1,1)) end if ! tracer mass rmc => chem_dat(region)%rm CALL GATHER( dgrid(region), rmc, glbfield, 0, status ) IF_NOTOK_RETURN(status=1) if ( isRoot ) then ! write tracers call io_write( hdf%id, imr,'LON'//region_name(region),& jmr,'LAT'//region_name(region), & lmr,'HYBRID',ntrace_chem,'NTRACE_CHEM',glbfield,'rmc') end if ! clear nullify(rmc) deallocate(glbfield) end if ! ntrace_chem > 0 ! Close !------- if ( isRoot ) then call Done( hdf, status ) IF_NOTOK_RETURN(status=1) end if #else write (gol,'("not compiled with hdf4")'); call goErr TRACEBACK; status=1; return #endif end do !regions... ! ------- ! Done ! ------- call DoneRc( rcF, status ) IF_NOTOK_RETURN(status=1) ! ok status = 0 END SUBROUTINE WRITE_SAVE_FILE !EOC END MODULE IO_SAVE