#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_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; if (isRoot) call MDF_CLose(fid,status); status=1; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if #include "tm5.inc" !---------------------------------------------------------------------------- ! TM5 ! !---------------------------------------------------------------------------- !BOP ! ! !MODULE: RESTART ! ! !DESCRIPTION: Write and read restart files. !\\ !\\ ! !INTERFACE: ! MODULE RESTART ! ! !USES: ! use GO , only : gol, goPr, goErr use dims , only : nregions implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Restart_Init ! read restart keys in rc file public :: Restart_Done ! nothing yet public :: Restart_Save ! wrapper around Restart_Write public :: Restart_Write ! write a restart file public :: Restart_Read ! read a restart file public :: rs_write ! model must write restart ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'Restart' character(len=256) :: rs_write_dir logical :: rs_write logical :: rs_write_extra integer :: rs_write_extra_dhour, rs_write_extra_hour integer :: fid ! file id for IF_NOTOK_MDF macro ! ! !REVISION HISTORY: ! 8 Apr 2011 - P. Le Sager - Close MDF file if error occurs. This is ! needed for mpi_abort not to hang. See TM5_MPI_Abort in ! partools, and remarks below. Made IF_NOTOK_MDF macro for ! that purpose. ! 28 Apr 2011 - P. Le Sager - Read method : handle restart file with extra ! tracers. ! 10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) when an error occurs when accessing MDF files, you should first close ! the file before returning. The IF_NOTOK_MDF macro takes care of that. ! The only thing you need is to call it like that : ! ! IF_NOTOK_MDF(fid=xxxx) ! ! where you replace xxxx with the integer id (file handler) of the file ! you are accessing. Note that this does not solve all problems (but ! probably most of them): it is still possible that MDF_Close hangs... ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_INIT ! ! !DESCRIPTION: read settings from rcfile !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_INIT( status ) ! ! !USES: ! use GO , only : TrcFile, Init, Done, ReadRc use global_data, only : rcfile use global_data, only : outdir use meteodata , only : lli ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = 'Restart_Init' type(TrcFile) :: rcF ! ---- begin call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) ! write restart files at all ? call ReadRc( rcF, 'restart.write', rs_write, status, default=.false. ) IF_ERROR_RETURN(status=1) ! further settings ... if ( rs_write ) then ! output directory: call ReadRc( rcF, 'restart.write.dir', rs_write_dir, status, default=outdir ) IF_ERROR_RETURN(status=1) ! extra restart files ? call ReadRc( rcF, 'restart.write.extra', rs_write_extra, status, default=.false. ) IF_ERROR_RETURN(status=1) if ( rs_write_extra ) then call ReadRc( rcF, 'restart.write.extra.hour', rs_write_extra_hour, status, default=0 ) IF_ERROR_RETURN(status=1) call ReadRc( rcF, 'restart.write.extra.dhour', rs_write_extra_dhour, status, default=24 ) IF_ERROR_RETURN(status=1) end if end if ! write restart files call Done( rcF, status ) IF_NOTOK_RETURN(status=1) status = 0 END SUBROUTINE RESTART_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_DONE ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_DONE( status ) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = 'Restart_Done' ! --- begin -------------------------------- ! nothing to be done ... ! ok status = 0 END SUBROUTINE RESTART_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_SAVE ! ! !DESCRIPTION: !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_SAVE( status, extra, isfirst ) ! ! !USES: ! use dims, only : idate ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! logical, intent(in), optional :: extra logical, intent(in), optional :: isfirst ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = 'Restart_Save' logical :: is_extra real :: t1, t2 ! --- begin -------------------------------- ! options ... is_extra = .false. if ( present(extra) ) is_extra = extra ! write restart files at all ? if ( rs_write ) then ! end or extra ? if ( is_extra ) then ! save extra restart files ? if ( rs_write_extra ) then ! every hour+n*dhour only : if ( modulo( idate(4) - rs_write_extra_hour, rs_write_extra_dhour ) == 0 .and. & all( idate(5:6) == 0 ) ) then ! write restart file for this time: call Restart_Write( status, isfirst=isfirst ) IF_NOTOK_RETURN(status=1) end if ! for this hour end if ! extra restart files ? else ! write restart file : call cpu_time(t1) call Restart_Write( status, isfirst=isfirst ) IF_NOTOK_RETURN(status=1) call cpu_time(t2) write (gol,*) " time to write restart [s]: ", t2-t1 ; call goPr end if ! not extra end if ! write at all ! ok status = 0 END SUBROUTINE RESTART_SAVE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_FILENAME ! ! !DESCRIPTION: Build restart filename from inputs. !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_FILENAME( region, fname, status, key, dir, isfirst ) ! ! !USES: ! use dims , only : idate use global_data, only : outdir use meteodata , only : lli ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region logical, intent(in), optional :: isfirst character(len=*), intent(in), optional :: dir character(len=*), intent(in), optional :: key ! ! !OUTPUT PARAMETERS: ! character(len=*), intent(out) :: fname integer, intent(out) :: status ! ! !REVISION HISTORY: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = 'Restart_FileName' character(len=256) :: adir character(len=32) :: akey ! --- begin -------------------------------- ! destination directory: adir = trim(outdir) if ( present(dir) ) adir = trim(dir) ! extra key, for example '_x' to denote that ! a restart file was dumped after process 'x': akey = '' if ( present(key) ) akey = trim(key) ! if this is the initial time, add an extra key to avoid ! that the restart file for this hour from the previous ! run is overwritten: if ( present(isfirst) ) then if ( isfirst ) akey = trim(akey)//'_initial' end if ! write filename: write (fname,'(a,"/TM5_restart_",i4.4,2i2.2,"_",2i2.2,"_",a,a,".nc")') & trim(adir), idate(1:5), trim(lli(region)%name), trim(akey) ! ok status = 0 END SUBROUTINE RESTART_FILENAME !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_WRITE ! ! !DESCRIPTION: write restart !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_WRITE( status, key, region, isfirst ) ! ! !USES: ! use GO , only : Get use dims , only : nregions, at, bt use chem_param , only : ntracet, ntrace_chem, ntrace, names use partools , only : isRoot use tm5_distgrid, only : dgrid, Get_DistGrid, gather use global_data , only : mass_dat, chem_dat use meteodata , only : global_lli, levi use meteodata , only : sp_dat, phlb_dat, m_dat use MDF , only : MDF_Create, MDF_EndDef, MDF_Close use MDF , only : MDF_Def_Dim, MDF_Def_Var use MDF , only : MDF_Put_Att, MDF_Put_Var use MDF , only : MDF_REPLACE, MDF_NETCDF4 use MDF , only : MDF_FLOAT, MDF_DOUBLE, MDF_CHAR ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! character(len=*), intent(in), optional :: key integer, intent(in), optional :: region logical, intent(in), optional :: isfirst ! ! !REVISION HISTORY: ! 8 Apr 2011 - P. Le Sager - use IF_NOTOK_MDF macro ! 10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = 'Restart_Write' integer :: imr, jmr, lmr, n character(len=256) :: fname integer :: ftype integer :: ncid integer :: dimid_lon, dimid_lat, dimid_lev, dimid_hlev integer :: dimid_lon_sfc, dimid_lat_sfc integer :: dimid_trace, dimid_trace_transp, dimid_trace_chem integer :: dimid_name integer :: varid, varid_at, varid_bt integer :: varid_sp, varid_ph, varid_m integer :: varid_names, varid_rm #ifdef slopes integer :: varid_rxm, varid_rym, varid_rzm #endif integer :: varid_rmc integer :: rtype real, allocatable :: arr4d(:,:,:,:), arr3d(:,:,:) ! --- begin -------------------------------- write (gol,'("write restart file(s) ...")'); call goPr ! loop over regions: REG: do n = 1, nregions ! only selected region ? if ( present(region) ) then if ( n /= region ) cycle end if ! entire region grid size imr = global_lli(n)%nlon jmr = global_lli(n)%nlat lmr = levi%nlev ! allocate 3D and 4D global arrays for gathering data if (isRoot) then allocate( arr4d(imr,jmr,lmr,ntracet) ) allocate( arr3d(imr,jmr,lmr+1) ) else allocate( arr4d(1,1,1,1) ) allocate( arr3d(1,1,1) ) endif ! name of restart file call Restart_FileName( n, fname, status, key=key, dir=rs_write_dir, isfirst=isfirst ) IF_NOTOK_RETURN(status=1) write (gol,'(" destination : ",a)') trim(fname); call goPr if (isRoot) then !------------------ ! OPEN NETCDF FILE !------------------ ! overwrite existing files (clobber) call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, ncid, status ) IF_NOTOK_RETURN(status=1) !------------------ ! DEFINE DIMENSIONS !------------------ call MDF_Def_Dim( ncid, 'lon', imr, dimid_lon, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Def_Dim( ncid, 'lat', jmr, dimid_lat, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Def_Dim( ncid, 'lev', lmr, dimid_lev, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Def_Dim( ncid, 'hlev', lmr+1, dimid_hlev, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Def_Dim( ncid, 'trace_transp', ntracet, dimid_trace_transp, status ) IF_NOTOK_MDF(fid=ncid) if ( ntrace_chem > 0 ) then call MDF_Def_Dim( ncid, 'trace_chem', ntrace_chem, dimid_trace_chem, status ) IF_NOTOK_MDF(fid=ncid) else dimid_trace_chem = -1 end if call MDF_Def_Dim( ncid, 'trace', ntrace, dimid_trace, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Def_Dim( ncid, 'name', len(names(1)), dimid_name, status ) IF_NOTOK_MDF(fid=ncid) !------------------ ! DEFINE VARIABLES !------------------ select case ( kind(m_dat(n)%data) ) case ( 4 ) ; rtype = MDF_FLOAT case ( 8 ) ; rtype = MDF_DOUBLE case default write (gol,'("unsupported real kind : ",i6)') kind(m_dat(n)%data) TRACEBACK; status=1; return end select ! surface pressure call MDF_Def_Var( ncid, 'sp', rtype, (/dimid_lon,dimid_lat/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'surface pressure', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'Pa', status ) IF_NOTOK_MDF(fid=ncid) varid_sp = varid ! at, bt coefficients for hybrid grid call MDF_Def_Var( ncid, 'at', rtype, (/dimid_hlev/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'hybrid grid a_t coefficient', status ) IF_NOTOK_MDF(fid=ncid) varid_at = varid call MDF_Def_Var( ncid, 'bt', rtype, (/dimid_hlev/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'hybrid grid b_t coefficient', status ) IF_NOTOK_MDF(fid=ncid) varid_bt = varid ! half level pressure call MDF_Def_Var( ncid, 'ph', rtype, & (/dimid_lon,dimid_lat,dimid_hlev/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'half level pressure', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'Pa', status ) IF_NOTOK_MDF(fid=ncid) varid_ph = varid ! air mass call MDF_Def_Var( ncid, 'm', rtype, & (/dimid_lon,dimid_lat,dimid_lev/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'air mass', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg', status ) IF_NOTOK_MDF(fid=ncid) varid_m = varid !! accumulated surface fluxes !! !call MDF_Def_Var( ncid, 'slhf', rtype, (/dimid_lon_sfc,dimid_lat_sfc/), varid, status ) !IF_NOTOK_MDF(fid=ncid) !call MDF_Put_Att( ncid, varid, 'long_name', 'surface latent heat flux', status ) !IF_NOTOK_MDF(fid=ncid) !call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status ) !IF_NOTOK_MDF(fid=ncid) !varid_slhf = varid !! !call MDF_Def_Var( ncid, 'sshf', rtype, (/dimid_lon_sfc,dimid_lat_sfc/), varid, status ) !IF_NOTOK_MDF(fid=ncid) !call MDF_Put_Att( ncid, varid, 'long_name', 'surface sensible heat flux', status ) !IF_NOTOK_MDF(fid=ncid) !call MDF_Put_Att( ncid, varid, 'unit', 'W/m2', status ) !IF_NOTOK_MDF(fid=ncid) !varid_sshf = varid ! tracer names call MDF_Def_Var( ncid, 'names', MDF_CHAR, (/dimid_name,dimid_trace/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'tracer names', status ) IF_NOTOK_MDF(fid=ncid) varid_names = varid ! tracer mass call MDF_Def_Var( ncid, 'rm', rtype, & (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'transported tracer mass', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg', status ) IF_NOTOK_MDF(fid=ncid) varid_rm = varid ! tracer mass slopes: #ifdef slopes call MDF_Def_Var( ncid, 'rxm', rtype, & (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in x direction', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status ) IF_NOTOK_MDF(fid=ncid) varid_rxm = varid call MDF_Def_Var( ncid, 'rym', rtype, & (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in y direction', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status ) IF_NOTOK_MDF(fid=ncid) varid_rym = varid call MDF_Def_Var( ncid, 'rzm', rtype, & (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_transp/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'tracer mass slope in z direction', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg/(half cell)', status ) IF_NOTOK_MDF(fid=ncid) varid_rzm = varid #endif ! non-transported tracers: if ( ntrace_chem > 0 ) then call MDF_Def_Var( ncid, 'rmc', rtype, & (/dimid_lon,dimid_lat,dimid_lev,dimid_trace_chem/), varid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'long_name', 'non-transported tracer mass', status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Att( ncid, varid, 'unit', 'kg', status ) IF_NOTOK_MDF(fid=ncid) varid_rmc = varid end if !------------------ ! END DEFINITION MODE !------------------ call MDF_EndDef( ncid, status ) IF_NOTOK_MDF(fid=ncid) endif !------------------ ! WRITE VARIABLES !------------------ ! surface pressure call gather( dgrid(n), sp_dat(n)%data, arr3d(:,:,1:1), sp_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_sp, arr3d(:,:,1), status ) IF_NOTOK_MDF(fid=ncid) ! half level pressure call gather( dgrid(n), phlb_dat(n)%data, arr3d, phlb_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_ph, arr3d, status) IF_NOTOK_MDF(fid=ncid) ! at, bt coefficients if (isRoot) then call MDF_Put_Var( ncid, varid_at, at(1:lmr+1), status ) IF_NOTOK_MDF(fid=ncid) call MDF_Put_Var( ncid, varid_bt, bt(1:lmr+1), status ) IF_NOTOK_MDF(fid=ncid) end if ! air mass call gather( dgrid(n), m_dat(n)%data, arr4d(:,:,:,1), m_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_m, arr4d(:,:,:,1), status) IF_NOTOK_MDF(fid=ncid) !! surface latent heat flux; global surface field ! !call MDF_Put_Var( ncid, varid_slhf, slhf_dat(iglbsfc)%data(1:n360,1:n180,1), status ) !IF_NOTOK_MDF(fid=ncid) ! !! surface sensible heat flux; global surface field ! !call MDF_Put_Var( ncid, varid_sshf, sshf_dat(iglbsfc)%data(1:n360,1:n180,1), status ) !IF_NOTOK_MDF(fid=ncid) ! tracer names if (isRoot) call MDF_Put_Var( ncid, varid_names, names, status ) IF_NOTOK_MDF(fid=ncid) ! write transported tracers call gather( dgrid(n), mass_dat(n)%rm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rm, arr4d, status) IF_NOTOK_MDF(fid=ncid) #ifdef slopes call gather( dgrid(n), mass_dat(n)%rxm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rxm, arr4d, status) IF_NOTOK_MDF(fid=ncid) call gather( dgrid(n), mass_dat(n)%rym, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rym, arr4d, status) IF_NOTOK_MDF(fid=ncid) call gather( dgrid(n), mass_dat(n)%rzm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rzm, arr4d, status) IF_NOTOK_MDF(fid=ncid) #endif ! write transported tracers call gather( dgrid(n), mass_dat(n)%rm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rm, arr4d, status) IF_NOTOK_MDF(fid=ncid) #ifdef slopes call gather( dgrid(n), mass_dat(n)%rxm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rxm, arr4d, status) IF_NOTOK_MDF(fid=ncid) call gather( dgrid(n), mass_dat(n)%rym, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rym, arr4d, status) IF_NOTOK_MDF(fid=ncid) call gather( dgrid(n), mass_dat(n)%rzm, arr4d, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rzm, arr4d, status) IF_NOTOK_MDF(fid=ncid) #endif ! write non-transported tracers if (ntrace_chem > 0) then call gather( dgrid(n), chem_dat(n)%rm, arr4d(:,:,:,1:ntrace_chem), chem_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if (isRoot) call MDF_Put_Var( ncid, varid_rmc, arr4d(:,:,:,1:ntrace_chem), status) IF_NOTOK_MDF(fid=ncid) end if ! Done if (isRoot) call MDF_Close( ncid, status ) IF_NOTOK_RETURN(status=1) deallocate(arr4d, arr3d) end do REG status = 0 END SUBROUTINE RESTART_WRITE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: RESTART_READ ! ! !DESCRIPTION: Read restart file. Case of istart=33 (can read any of the ! available variables) or 32 (can read only tracer mass). !\\ !\\ ! !INTERFACE: ! SUBROUTINE RESTART_READ( status, region, & surface_pressure, pressure, air_mass, surface_fluxes, & tracer_mass, tendencies, megan_history, nox_pulsing ) ! ! !USES: ! use GO, only : TrcFile, Init, Done, ReadRc use GO, only : goMatchValue use dims, only : nregions, im, jm, istart, idate, idatei use grid, only : TllGridInfo, TLevelInfo, Init, Done, Fill3D use chem_param, only : ntracet, ntrace_chem, ntrace use chem_param, only : names, tracer_name_len use partools, only : isRoot, par_broadcast use tm5_distgrid, only : dgrid, gather, scatter use global_data, only : rcfile, mass_dat, chem_dat use meteodata, only : levi, global_lli, sp_dat, phlb_dat, m_dat !use meteodata, only : slhf_dat, sshf_dat use MDF, only : MDF_Open, MDF_Close, MDF_Inquire_Dimension use MDF, only : MDF_Inq_VarID, MDF_Inquire_Variable, MDF_Inq_DimID use MDf, only : MDF_Var_Par_Access, MDF_INDEPENDENT, MDF_COLLECTIVE use MDF, only : MDF_Get_Att, MDF_Get_Var use MDF, only : MDF_READ, MDF_NETCDF4 ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! integer, intent(in), optional :: region logical, intent(in), optional :: surface_pressure, pressure, air_mass, surface_fluxes logical, intent(in), optional :: tracer_mass, tendencies, megan_history, nox_pulsing ! ! !REVISION HISTORY: ! 8 Apr 2011 - P. Le Sager - use IF_NOTOK_MDF macro ! 28 Apr 2011 - P. Le Sager - Check on tracer availability in restart file. ! - Allows for more tracers in restart file than needed ! 10 May 2011 - P. Le Sager - Added deallocate statement to work with zoom regions ! 10 Oct 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! - logically, if we need to remap, then meteo is not read from restart ! (but from met field and used for remapping): in other words, only ! tracer mass is read, and istart should be 32. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Restart_Read' character(len=tracer_name_len), allocatable :: values_names(:) character(len=256) :: rs_read_dir, fname type(TrcFile) :: rcF logical :: exist logical :: do_sp, do_ph, do_m, do_sflux, do_rm, do_megan, do_pulse integer :: imr, jmr, lmr, imr_restart, jmr_restart, lmr_restart integer :: n, region1, region2 integer :: ncid integer :: varid_sp, varid_ph, varid_m, varid_rm, varid_rmc, varid_names !integer :: varid_slhf, varid_sshf integer :: itr, itr_file integer :: ntracet_restart, dimid integer :: shp(2) #ifdef slopes integer :: varid_rxm, varid_rym, varid_rzm #endif ! global work arrays to read data real, allocatable :: tmp3d(:,:,:) real, allocatable :: rmt(:,:,:,:),rms(:,:,:,:), rmx(:,:,:,:),rmy(:,:,:,:), rmz(:,:,:,:) ! for remapping: logical :: need_vremap, need_hremap, need_remap integer :: varid_at, varid_bt real :: dx, dy real, allocatable :: sp_gbl(:,:,:) real, allocatable :: at_restart(:), bt_restart(:) real, allocatable :: src_glb(:,:,:) type(TllGridInfo) :: lli_restart type(TLevelInfo) :: levi_restart ! --- begin -------------------------------- if ( istart /= 33 .and. istart /= 32 ) then write (gol,'(" skip read restart; istart not 33 or 32 but ",i2)') istart; call goPr status=0; return endif if ( any( idate /= idatei ) ) then write (gol,'(" skip read restart; idate not idatei but ",i4,5i2.2)') idate; call goPr status=0; return endif ! input directory: call Init( rcF, rcfile, status ) IF_NOTOK_RETURN(status=1) call ReadRc( rcF, 'restart.read.dir', rs_read_dir, status ) IF_NOTOK_RETURN(status=1) call Done( rcF, status ) IF_NOTOK_RETURN(status=1) ! region range: if ( present(region) ) then region1 = region region2 = region else region1 = 1 region2 = nregions end if ! data sets: do_rm = .false. ; if ( present(tracer_mass ) ) do_rm = tracer_mass do_m = .false. ; if ( present(air_mass ).and.(istart==33) ) do_m = air_mass do_sp = .false. ; if ( present(surface_pressure ).and.(istart==33) ) do_sp = surface_pressure do_ph = .false. ; if ( present(pressure ).and.(istart==33) ) do_ph = pressure do_sflux = .false. ; if ( present(surface_fluxes ).and.(istart==33) ) do_sflux = surface_fluxes do_megan = .false. ; if ( present(megan_history ).and.(istart==33) ) do_megan = megan_history do_pulse = .false. ; if ( present(nox_pulsing ).and.(istart==33) ) do_pulse = nox_pulsing ! sorry .. if ( do_sflux ) then write (gol,'("no surface fluxes in restart files until somebody")') ; call goErr write (gol,'("has a good idea on what should be storred:")') ; call goErr write (gol,'(" o global surface field (1x1 ?)")') ; call goErr write (gol,'(" o zoom regions")') ; call goErr write (gol,'(" o both")') ; call goErr TRACEBACK; status=1; return end if ! do we need anything? if(.not.(do_rm.or.do_m.or.do_sp.or.do_ph.or.do_sflux.or.do_megan.or.do_pulse))then status=0; return endif REG: do n = region1, region2 imr = global_lli(n)%nlon jmr = global_lli(n)%nlat lmr = levi%nlev ! name of restart file call Restart_FileName( n, fname, status, dir=trim(rs_read_dir) ) IF_NOTOK_RETURN(status=1) write (gol,'(" read restart file: ",a)') trim(fname); call goPr inquire( file=fname, exist=exist ) if ( .not. exist ) then write (gol,'("restart file not found : ",a)') trim(fname); call goErr TRACEBACK; status=1; return end if ! ** open netcdf file if (isRoot) then call MDF_Open( trim(fname), MDF_NETCDF4, MDF_READ, ncid, status ) IF_NOTOK_RETURN(status=1) ! ** check for dimension compatibility call MDF_Inq_DimID( ncid, 'lev', dimid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inquire_Dimension( ncid, dimid, status, length=lmr_restart ) IF_NOTOK_MDF(fid=ncid) call MDF_Inq_DimID( ncid, 'lat', dimid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inquire_Dimension( ncid, dimid, status, length=jmr_restart ) IF_NOTOK_MDF(fid=ncid) call MDF_Inq_DimID( ncid, 'lon', dimid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inquire_Dimension( ncid, dimid, status, length=imr_restart ) IF_NOTOK_MDF(fid=ncid) need_vremap = (lmr /= lmr_restart) need_hremap = (jmr /= jmr_restart) .or. (imr /= imr_restart) need_remap = need_hremap .or. need_vremap endif call par_broadcast( need_remap, status) IF_NOTOK_RETURN(status=1) if ((istart /= 32).and.need_remap) then status=1 write(gol,*)' you must use istart=32 for using a restart file at different resolution' call goErr TRACEBACK; return endif ! work arrays if (isRoot) then allocate( rmt(imr,jmr,lmr,ntracet) ) allocate( rmx(imr,jmr,lmr,ntracet) ) allocate( rmy(imr,jmr,lmr,ntracet) ) allocate( rmz(imr,jmr,lmr,ntracet) ) if ( ntrace_chem > 0 ) allocate( rms(imr,jmr,lmr,ntracet+1:ntracet+ntrace_chem) ) allocate( tmp3d(imr,jmr,lmr+1 ) ) else allocate( rmt(1,1,1,1) ) allocate( rmx(1,1,1,1) ) allocate( rmy(1,1,1,1) ) allocate( rmz(1,1,1,1) ) if ( ntrace_chem > 0 ) allocate( rms(1,1,1,1) ) allocate( tmp3d(1,1,1) ) endif ! prepare for remap if (need_remap .and. do_rm) then write (gol,'(" remap tracer from restart file")') ; call goPr if (isRoot) then allocate( sp_gbl(imr,jmr,1) ) allocate( src_glb(imr_restart,jmr_restart,lmr_restart)) else allocate(sp_gbl(1,1,1)) allocate(src_glb(1,1,1)) endif call gather( dgrid(n), sp_dat(n)%data, sp_gbl, sp_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) ! init to 0 in case of data not found in file rmt=0. rms=0. ! init lli_restart/levi_restart allocate(at_restart(lmr_restart+1)) allocate(bt_restart(lmr_restart+1)) ! call MDF_Inq_VarID( ncid, 'at', varid_at, status ) IF_NOTOK_MDF(fid=ncid) ! call MDF_Get_Var( ncid, varid_at, at_restart(1:(lmr_restart+1)), status ) IF_NOTOK_MDF(fid=ncid) ! call MDF_Inq_VarID( ncid, 'bt', varid_bt, status ) IF_NOTOK_MDF(fid=ncid) ! call MDF_Get_Var( ncid, varid_bt, bt_restart(1:(lmr_restart+1)), status ) IF_NOTOK_MDF(fid=ncid) ! call Init( levi_restart, lmr_restart, at_restart, bt_restart, status ) IF_NOTOK_RETURN(status=1) ! deallocate(at_restart,bt_restart) ! dx=360./imr_restart dy=180./jmr_restart call Init( lli_restart, -180.+0.5*dx, dx, imr_restart, & -90.+0.5*dy, dy, jmr_restart, status ) IF_NOTOK_RETURN(status=1) endif ! ** get variables id if (isRoot) then ! surface pressure if ( do_sp ) call MDF_Inq_VarID( ncid, 'sp', varid_sp, status ) IF_NOTOK_MDF(fid=ncid) ! half level pressure if ( do_ph ) call MDF_Inq_VarID( ncid, 'ph', varid_ph, status ) IF_NOTOK_MDF(fid=ncid) ! air mass if ( do_m ) call MDF_Inq_VarID( ncid, 'm', varid_m, status ) IF_NOTOK_MDF(fid=ncid) !! surface fluxes !if ( do_sflux ) then !end if ! tracer mass if ( do_rm ) then call MDF_Inq_VarID( ncid, 'names', varid_names, status ) if ( status /= 0 ) then write (gol,'("could not find variable `names` in restart file;")'); call goErr write (gol,'(" using an old restart file to initialize the model ?")'); call goErr status=1 end if IF_NOTOK_MDF(fid=ncid) ! get dimension of "names" call MDF_Inquire_Variable( ncid, varid_names, status, shp=shp ) IF_NOTOK_MDF(fid=ncid) ! get number of transported tracer in restart file call MDF_Inq_DimID( ncid, 'trace_transp', dimid, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inquire_Dimension( ncid, dimid, status, length=ntracet_restart ) IF_NOTOK_MDF(fid=ncid) ! tracers mass id call MDF_Inq_VarID( ncid, 'rm', varid_rm, status ) IF_NOTOK_MDF(fid=ncid) #ifdef slopes call MDF_Inq_VarID( ncid, 'rxm', varid_rxm, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inq_VarID( ncid, 'rym', varid_rym, status ) IF_NOTOK_MDF(fid=ncid) call MDF_Inq_VarID( ncid, 'rzm', varid_rzm, status ) IF_NOTOK_MDF(fid=ncid) #endif ! read non-transported tracers if any if ( ntrace_chem > 0 ) then call MDF_Inq_VarID( ncid, 'rmc', varid_rmc, status ) IF_NOTOK_MDF(fid=ncid) end if end if end if ! *** READ VARIABLES *** if ( do_sp ) then write (gol,'(" restore surface pressure ...")'); call goPr if (isRoot) call MDF_Get_Var( ncid, varid_sp, tmp3d(:,:,1), status ) IF_NOTOK_MDF(fid=ncid) call scatter( dgrid(n), sp_dat(n)%data, tmp3d(:,:,1:1), sp_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) end if if ( do_ph ) then write (gol,'(" restore half level pressure ...")'); call goPr if (isRoot) call MDF_Get_Var( ncid, varid_ph, tmp3d, status) IF_NOTOK_MDF(fid=ncid) call scatter( dgrid(n), phlb_dat(n)%data, tmp3d, phlb_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) end if if ( do_m ) then write (gol,'(" restore air mass ...")'); call goPr if (isRoot) call MDF_Get_Var( ncid, varid_m, tmp3d(:,:,1:lmr), status) IF_NOTOK_MDF(fid=ncid) call scatter( dgrid(n), m_dat(n)%data, tmp3d(:,:,1:lmr), m_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) end if ! tracer mass READRM: if ( do_rm ) then write (gol,'(" restore tracer mass ...")'); call goPr ! read list with tracer names in rcfile allocate( values_names(shp(2)) ) call MDF_Get_Var( ncid, varid_names, values_names, status ) IF_NOTOK_MDF(fid=ncid) ! loop over all model tracers do itr = 1, ntrace if (isRoot) then ! search in list: call goMatchValue( names(itr), values_names, itr_file, status ) if ( status < 0 ) then write(gol,'("*WARNING* Requested tracer `",a,"` not FOUND in restart file!")') trim(names(itr)) if (istart /= 32) then call goErr IF_NOTOK_MDF(fid=ncid) else status=0 call goPr if ( itr <= ntracet ) then rmt(:,:,:,itr) = 1.e-30 write(gol,'("*WARNING* Requested TRANSPORTED tracer `",a,"` has been SET to a default value of 1.e30")') trim(names(itr)) else rms(:,:,:,itr) = 1.e-30 write(gol,'("*WARNING* Requested SHORT-LIVED tracer `",a,"` has been SET to a default value of 1.e30")') trim(names(itr)) endif call goPr endif else ! transported or short lived tracer ? if ( itr <= ntracet ) then if ( itr_file > ntracet_restart ) then write (gol,'("tracer `",a,"` is transported but seems to be not-transported in restart file")') trim(names(itr)); call goErr status=1 IF_NOTOK_MDF(fid=ncid) end if if (need_remap) then call MDF_Get_Var( ncid, varid_rm, src_glb, status, start=(/1,1,1,itr_file/)) IF_NOTOK_MDF(fid=ncid) call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rmt(:,:,:,itr), & lli_restart, levi_restart, src_glb, 'sum', status ) IF_NOTOK_RETURN(status=1) else call MDF_Get_Var( ncid, varid_rm, rmt(:,:,:,itr), status, start=(/1,1,1,itr_file/)) IF_NOTOK_MDF(fid=ncid) endif #ifdef slopes ! read slopes if (.not. need_remap) then if (isRoot) call MDF_Get_Var( ncid, varid_rxm, rmx(:,:,:,itr), status, start=(/1,1,1,itr_file/)) IF_NOTOK_MDF(fid=ncid) if (isRoot) call MDF_Get_Var( ncid, varid_rym, rmy(:,:,:,itr), status, start=(/1,1,1,itr_file/)) IF_NOTOK_MDF(fid=ncid) if (isRoot) call MDF_Get_Var( ncid, varid_rzm, rmz(:,:,:,itr), status, start=(/1,1,1,itr_file/)) IF_NOTOK_MDF(fid=ncid) ! Scale methane concentration slopes by a factor specified in the rc file if ( (factor_ch4 /= 1.) .and. (itr == ich4) ) then mass_dat(n)%rxm(:,:,:,itr)= mass_dat(n)%rxm(:,:,:,itr) * factor_ch4 mass_dat(n)%rym(:,:,:,itr)= mass_dat(n)%rym(:,:,:,itr) * factor_ch4 mass_dat(n)%rzm(:,:,:,itr)= mass_dat(n)%rzm(:,:,:,itr) * factor_ch4 endif endif #endif else ! short lived tracer: if ( itr_file <= ntracet_restart ) then write (gol,'("tracer `",a,"` is not-transported but seems to be transported in restart file")') trim(names(itr)); call goErr status=1 IF_NOTOK_MDF(fid=ncid) end if itr_file = itr_file - ntracet_restart if (need_remap) then call MDF_Get_Var( ncid, varid_rmc, src_glb, status, start=(/1,1,1,itr_file/) ) IF_NOTOK_MDF(fid=ncid) call Fill3D( global_lli(n), levi, 'n', sp_gbl(:,:,1), rms(:,:,:,itr), & lli_restart, levi_restart, src_glb, 'sum', status ) IF_NOTOK_RETURN(status=1) else call MDF_Get_Var( ncid, varid_rmc, rms(:,:,:,itr), status, start=(/1,1,1,itr_file/) ) IF_NOTOK_MDF(fid=ncid) endif end if ! transported or short-lived endif ! in the file endif ! root end do ! tracers ! distribute call scatter( dgrid(n), mass_dat(n)%rm, rmt, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) if ( ntrace_chem > 0 ) then call scatter( dgrid(n), chem_dat(n)%rm, rms, chem_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) endif #ifdef slopes if (.not. need_remap) then call scatter( dgrid(n), mass_dat(n)%rxm, rmx, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call scatter( dgrid(n), mass_dat(n)%rym, rmy, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) call scatter( dgrid(n), mass_dat(n)%rzm, rmz, mass_dat(n)%halo, status) IF_NOTOK_RETURN(status=1) else ! Ensure that slopes are initialized to "unset" values of 0.0. Wouter says that ! we could remap levels for rxm et al., but 0.0 will also work. The noise ! induced from remapping the rm array is almost certainly bigger than any issues ! from having this "default=0.0" slopes information. -ARJ 1 Jan 12 mass_dat(n)%rxm = 0.0 mass_dat(n)%rym = 0.0 mass_dat(n)%rzm = 0.0 endif #endif ! free mem for next region deallocate( values_names) if (need_remap) deallocate(sp_gbl,src_glb) ENDIF READRM ! clean deallocate(rmt) if ( ntrace_chem > 0 ) deallocate(rms) #ifdef slopes deallocate(rmx, rmy, rmz) #endif deallocate( tmp3d ) if (isRoot) call MDF_Close( ncid, status ) IF_NOTOK_RETURN(status=1) ENDDO REG status = 0 END SUBROUTINE RESTART_READ !EOC END MODULE RESTART