program prep_chimere ! Extracts the regional domain boundary conditions ! from coarse-mesh previous CHIMERE run ! ! 201701: G.Siour, L.Menut: Boundary conditions are now written as panels (top, west, east etc..) use netcdf use chimere_consts use chimere_params use chimere_common, only : read_nml_geometry, read_nml_aerosol, read_nml_chemistry use interp_tools use ncvar implicit none #define NCERR(lnum) if(ncstat/=NF90_NOERR) call nc_err(ncstat,lnum,__FILE__) !----------------------------------------------------------------------------------------- ! Input arguments integer :: id0 ! start date of the simulation integer :: nhourrun ! Model run duration integer :: nspecies,nbins_read integer :: iglobfile ! Current iteration ! number of species in the input fine (included species written per size bins) integer :: ninputspecies ! Whether to use global top concentrations when same vertical grid integer :: nest_topconc ! horizontal coordinates file for nest domain character(len=256) :: fnigeog character(len=256) :: fni character(len=256) :: fniVCO character(len=256) :: fnaerosol character(len=256) :: fnoutspec integer :: noutspec character(len=256) :: fno character(len=256) :: fnini character(len=256) :: fncb namelist /args/ id0,nhourrun,fnigeog,fni,fncb,iglobfile,nest_topconc,fniVCO,fnaerosol,fnoutspec,noutspec,fno,fnini ! Local constants character(len=*),parameter :: fniARGS = 'prep_chimere.nml' ! Constant netCDF attributes character(len=*),parameter :: title='CHIMERE SUITE' character(len=*),parameter :: subtitle='Fine boundary concentrations file' character(len=*),parameter :: generator='Generated by prep_CHIMERE' character(len=*),parameter :: conventions='' ! Local variables ! arrays real(kind=iprec),allocatable,dimension(:) :: wx_b,wy_b real(kind=iprec),allocatable,dimension(:) :: wx_l,wy_l integer,allocatable,dimension(:) :: ix_b,iy_b integer,allocatable,dimension(:) :: ix_l,iy_l real(kind=iprec),allocatable,dimension(:) :: a,b,ac,bc real(kind=iprec),allocatable,dimension(:,:,:) :: hlay real(kind=iprec),allocatable,dimension(:,:,:) :: airm real(kind=iprec),allocatable,dimension(:,:) :: psfc real(kind=iprec),allocatable,dimension(:,:,:) :: xin1 real(kind=iprec),allocatable,dimension(:,:) :: xout_bulk real(kind=iprec),allocatable,dimension(:,:) :: xout_lat real(kind=iprec),allocatable,dimension(:,:) :: buf_top real(kind=iprec),allocatable,dimension(:,:,:) :: buf_lat real(kind=iprec),allocatable,dimension(:,:,:) :: xin real(kind=iprec),allocatable,dimension(:,:) :: c_lon,c_lat real(kind=iprec),allocatable,dimension(:) :: xb,yb real(kind=iprec),allocatable,dimension(:,:,:) :: pc,pn real(kind=iprec),allocatable,dimension(:,:) :: ps real(kind=iprec),allocatable,dimension(:) :: c_cutoff,n_cutoff real(kind=iprec),allocatable,dimension(:,:,:,:) :: corners_3d ! Coarse mesh top boundary concentration from the BOUND.nc file real,allocatable,dimension(:,:) :: coarse_topconc ! dimensions integer :: c_nzonal,c_nmerid,c_nlevels,f_nlevels,c_times,dl integer :: c_vcoord integer :: n_ncutoff,c_ncutoff integer :: c_ncutoff_dimid integer :: c_cutoff_varid ! indexes integer :: nl,nx,ny,ilcoarse,ilnest,ihz,izo,ime,ispec,ivar,nbd,ihb integer :: ixx,iyy,lstix,lstiy,ibin,i, j, k, iaero, indic_aero real(kind=iprec) :: lxx,lyy ! misc integer :: dim1, dim2 integer,allocatable,dimension(:) :: ibvarid integer :: nvars integer :: ierr integer :: nv,ideb integer :: ihr integer :: ifn_v integer :: ifn_args integer :: ifn_aerosol integer :: ifnbc !file unit for coarse boundary condition list file character(len=dlen) :: datebuf character(len=dlen) :: c_start_date,c_end_date ! time limits of coarse-mesh file real(kind=iprec) :: c_deltax,c_deltay logical :: coarse_is_recti, found logical :: same_vert_grid real(kind=iprec) :: wx,wy logical :: exist_psfc ! (MM5 format) character(len=dlen) :: c_cur_date ! current time integer :: c_id0,c_ide ! time limits of coarse-mesh file ! (CHIMERE format) integer :: delta_time_idx ! time index difference between coarse- ! and fine-mesh files real(kind=iprec) :: cf,aread,bread,deltap real(kind=iprec),allocatable,dimension(:,:,:) :: bulkbuf real(kind=iprec),allocatable,dimension(:,:,:) :: buf3d character(len=128) :: hname,usrname,systime character(len=256) :: cwd character(len=splen),allocatable,dimension(:) :: specname,inpspecname,aero_name integer,allocatable,dimension(:) :: inpbin,tmpinpbin character(len=splen),allocatable,dimension(:) :: qqqname !name of the variables in the output file character(len=splen),allocatable,dimension(:) :: inpname !name of the variables in the input file character(len=splen),allocatable,dimension(:) :: varname character(len=splen) :: tmpname,tmpinpname character(len=2) :: strbin integer:: ifle character(len=1024),allocatable,dimension(:) :: bndname integer :: nbn ! Number of coarse mesh boundary files (aer, gas,dust) ! netCDF stuff integer :: c_id,f_id ! IDs of coarse-mesh/fine-mesh files integer,allocatable,dimension(:) ::cb_id ! IDs of coarse-mesh boundary conditions !IDs of species variable is coarse mesh boundary files integer,allocatable,dimension(:) ::cb_species_varid !number of species in each coarse/mesh boundary file (gas, aer, dust...) integer,allocatable,dimension(:) ::cb_speclen !IDs of species variable dimension in coarse/mesh boundary file integer,allocatable,dimension(:) ::cb_species_dimid integer,allocatable,dimension(:) ::cb_topconc_varid !IDs of top conc variable for each file integer :: ncstat ! return code for netCDF functions integer :: c_date_dimid,f_date_dimid ! Date string dimension IDs integer :: f_spstr_dimid ! Species string dimension ID integer :: f_species_dimid ! species dimension IDs integer :: c_time_dimid,f_time_dimid ! Time dimension IDs integer :: c_we_dimid ,f_we_dimid integer :: c_sn_dimid ,f_sn_dimid integer :: c_bt_dimid ,f_bt_dimid integer :: f_hb_dimid integer :: c_vc_dimid ! VCOORD level dimension IDs integer :: c_times_varid,f_times_varid ! Times variable IDs integer :: f_species_varid ! species variable IDs integer :: c_lon_varid, f_lon_varid ! longitude variable ID integer :: c_lat_varid, f_lat_varid ! latitude variable ID integer :: c_hlay_varid ! layer altitude variable ID integer :: c_airm_varid ! air density variable ID integer :: c_psfc_varid ! surface pressure variable ID integer :: c_avcoord_varid,c_bvcoord_varid integer :: c_varid ! current variable ID in coarse-mesh file character(len=1024) :: c_history integer :: c_hlen ! length of history attribute in ! coarse-mesh file integer :: ini_ncid ! ID of output file integer :: ini_time_dimid integer :: ini_date_dimid integer :: ini_splen_dimid integer :: ini_species_dimid integer :: ini_we_dimid integer :: ini_sn_dimid integer :: ini_bt_dimid integer :: ini_times_varid integer :: ini_lon_varid integer :: ini_lat_varid integer :: ini_species_varid integer,dimension(:),allocatable :: ini_spec_varid character(len=256) :: long_name character(len=16) :: units logical :: accfull ! conv: conversion factor (from ug m-3 to ppb or 1) ! molar: scaled molar mass to convert from ug m-3 to ppb real(kind=iprec),allocatable,dimension(:,:,:) :: conv real(kind=iprec),allocatable,dimension(:) :: molar ! functions character(len=dlen),external :: numeric2mm5date type(ncfile) :: nc_geog !------------------------------------------------------------------------------------------ ! Remark : ! The interpolation formula used later on is equivalent to ! xout=(1-wx)*xin(i) + wx*xin(i+1) ! it expects wx=0 for two coincident points ! Thus the "real" weight is not wx, but (1-wx) ! but for historical reasons, we keep calling wx what should be called (1-wx) ! Conversion factor for hydrostatic relation cf = g/(Rcte*7.2868e16) ! Read input arguments from namelist call opfi(ifn_args,fniARGS,'f','o') read(ifn_args,nml=args,iostat=ierr) if(ierr/=0) then print *,'prep_chimere : problem reading arguments namelist' call exit1('Exiting') endif close (ifn_args) ! Read chimere namelist: geometry and aerosol call read_nml_geometry call read_nml_aerosol call read_nml_chemistry call get_system(usrname,hname,systime,cwd) ! open coarse-mesh (input) file ncstat=nf90_open(fni,NF90_NOWRITE,c_id) NCERR(__LINE__) ! read attributes from coarse-mesh file call read_attributes ! How many variables in coarse file ? ncstat=nf90_inquire(c_id, nVariables=nvars) NCERR(__LINE__) ! read dimensions from coarse-mesh file call read_dims if (nest_topconc == 1) then ! Reading coarse run boundary conditions' list of files to get the top boundaries call opfi(ifnbc,fncb,'f','o') read(ifnbc,*) nbn allocate(bndname(nbn)) allocate(cb_id(nbn)) allocate(cb_species_varid(nbn)) allocate(cb_speclen(nbn)) allocate(cb_species_dimid(nbn)) allocate(cb_topconc_varid(nbn)) ! Skip records of previous iterations do i = 1, (iglobfile-1)*nbn read(ifnbc,*) end do ! Read current iteration's file names do i = 1,nbn read(ifnbc,'(A)',iostat=ierr) bndname(i) enddo close(ifnbc) end if ! Once we know dimensions, we can allocate storage allocate(wx_b(nhoriz)) allocate(wy_b(nhoriz)) allocate(wx_l(nhbound)) allocate(wy_l(nhbound)) allocate(c_lon(c_nzonal,c_nmerid)) allocate(c_lat(c_nzonal,c_nmerid)) allocate(xb(nhbound)) allocate(yb(nhbound)) allocate(ix_b(nhoriz)) allocate(iy_b(nhoriz)) allocate(ix_l(nhbound)) allocate(iy_l(nhbound)) allocate(ac(c_nlevels)) allocate(bc(c_nlevels)) allocate(xin(c_nzonal,c_nmerid,c_nlevels)) allocate(hlay(c_nzonal,c_nmerid,c_nlevels)) allocate(airm(c_nzonal,c_nmerid,c_nlevels)) allocate(psfc(c_nzonal,c_nmerid)) allocate(pc(c_nzonal,c_nmerid,c_nlevels)) allocate(ps(c_nzonal,c_nmerid)) allocate(c_cutoff(c_ncutoff)) allocate(conv(c_nzonal,c_nmerid,c_nlevels)) allocate(corners_3d(c_nzonal,c_nmerid,4,3)) if (nest_topconc == 1) allocate(coarse_topconc(c_nzonal,c_nmerid)) nbins_read=c_ncutoff-1 ! due to binnes aerosol species in input, the number of "real" variables ! including pxxDUST kind of variables is typically GREATER than nvars. ! An upper bound is given by nvars * nbins_read only. ! Read species names allocate(varname(noutspec)) allocate(qqqname(noutspec)) allocate(inpname(noutspec)) allocate(tmpinpbin(noutspec)) ninputspecies=0 ispec = 0 ! Check if there are vars with concentration per size bin (made with $accur = full) if(aero.ge.1) then accfull = .false. do ivar=1,nvars ncstat = nf90_inquire_variable(c_id, ivar, name = varname(ivar)) ncstat = nf90_get_att(c_id, ivar, 'long_name', long_name) if (index(long_name,'concentration per size bin') /= 0) then if (index(long_name,'Particle concentration per size bin')==0) then ! ignore number concentration accfull = .true. exit endif endif enddo if (.not.accfull) then print *, '***For NESTED runs you must have concentrations per size bin in your output file! ***' print *, ' BC file ', trim(fni), ' : ', nvars, ' variables' call exit1('The [$accur] flag in chimere.par must be "full" for the coarse run!') endif endif allocate(aero_name(nkc)) if(bins.eq.1) call read_aerosol do ivar=1,nvars do j=1,256 long_name(j:j) = ' ' enddo ncstat = nf90_inquire_variable(c_id, ivar, name = varname(ivar)) ncstat = nf90_get_att(c_id, ivar, 'long_name', long_name) ! Treating aerosol species written per bins if (index(long_name,'concentration per size bin') /= 0) then if (index(long_name,'Particle concentration per size bin')==0) then ! ignore number concentration indic_aero=0 do iaero=1,nkc if(varname(ivar)(2:len_trim(varname(ivar))).eq.trim(aero_name(iaero))) then indic_aero=1 endif enddo if(indic_aero==0) cycle ninputspecies=ninputspecies + 1 ! All species count for 1 even if they have multiple bins tmpname= varname(ivar) tmpname=tmpname(2:)//' ' ! remove initial 'b' from species name do i=1,nbins_read ispec=ispec + 1 call i2str(i,strbin) qqqname(ispec)='p'//strbin//tmpname inpname(ispec)=varname(ivar) tmpinpbin (ispec)=i enddo endif endif ! Treating other species if (index(long_name,'Concentration') /= 0) then if (index(long_name,'Particle concentration per size bin')==0) then ! ignore number concentration ninputspecies=ninputspecies + 1 ! All species count for 1 even if they have multiple bins ispec = ispec+1 qqqname(ispec) = varname(ivar) inpname(ispec)= varname(ivar) tmpinpbin(ispec)=0 endif endif nspecies = ispec enddo allocate(ibvarid(5*nspecies)) allocate(specname(nspecies)) allocate(molar(nspecies)) allocate(inpspecname(nspecies)) allocate(inpbin(nspecies)) specname = qqqname(1:nspecies) inpspecname = inpname(1:nspecies) inpbin=tmpinpbin(1:nspecies) if (nest_topconc == 1) then do i = 1,nbn print *, "Opening coarse run's boundary condition file", trim(bndname(i)) ncstat=nf90_open(bndname(i),NF90_NOWRITE,cb_id(i)) NCERR(__LINE__) enddo else print *, "top conc is taken from the coarse run's output: nest_topconc =", nest_topconc end if ! read molar mass call read_molar_mass do i=1,nspecies print*,'i', i, trim(specname(i)), molar(i) enddo ! Read auxilliary (non-concentration) variables ncstat=nf90_inq_varid(c_id,'lon',c_lon_varid) NCERR(__LINE__) ncstat=nf90_inq_varid(c_id,'lat',c_lat_varid) NCERR(__LINE__) ncstat=nf90_inq_varid(c_id,'a_vcoord',c_avcoord_varid) NCERR(__LINE__) ncstat=nf90_inq_varid(c_id,'b_vcoord',c_bvcoord_varid) NCERR(__LINE__) ncstat=nf90_inq_varid(c_id,'cut_off_diameters',c_cutoff_varid) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_lon_varid,c_lon) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_lat_varid,c_lat) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_avcoord_varid,ac) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_bvcoord_varid,bc) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_cutoff_varid,c_cutoff) NCERR(__LINE__) if(bins.eq.1) then if(n_ncutoff.ne.c_ncutoff) then call exit1('***ERROR in prep_chimere.f90 : INCONSISTENCY n_ncutoff =/= c_ncutoff') else do i=1,n_ncutoff if(abs(n_cutoff(i)-c_cutoff(i))/c_cutoff(i).gt.1e-4) then call exit1('***ERROR in prep_chimere.f90 : INCONSISTENCY n_ncutoff =/= c_ncutoff') endif enddo endif endif ! Reading nest grid coordinate call nc_get_struct(fnigeog, nc_geog) ! Read output grid call nc_get_coord(nc_geog) ! Reading nest vertical grid call read_nest_vertical_grid same_vert_grid = ((c_nlevels==f_nlevels) .and. all(a==ac) .and. all(b==bc)) ! the number of vertical levels from VCOORD of fine mesh is f_nlevels. It should be equal to nverti if (f_nlevels /= nverti) then print*, 'prep_chimere : Vertical levels inconsistency:', f_nlevels, nverti call exit1('Exiting') endif allocate(xin1(c_nzonal,c_nmerid,f_nlevels)) allocate(xout_bulk(nhoriz,f_nlevels)) allocate(xout_lat(nhbound,f_nlevels)) allocate(buf_top(nspecies,nhoriz)) allocate(buf_lat(nspecies,nhbound,f_nlevels)) allocate(bulkbuf(nhoriz,f_nlevels,nspecies)) allocate(buf3d(nc_geog%nx,nc_geog%ny,f_nlevels)) allocate(ini_spec_varid(nspecies)) ! check if time span for fine-mesh run is included in coarse-mesh file call check_time ! Defining boundary coordinates call define_bound_coords coarse_is_recti= test_rectilon(c_lon,c_nzonal,c_nmerid) .and. & & test_rectilat(c_lat,c_nzonal,c_nmerid) if (coarse_is_recti) then print *,'**** prep_chimere: coarse domain is rectilinear, fast !' c_deltax = (c_lon(c_nzonal, c_nmerid/2) - c_lon(1, c_nmerid/2))/(c_nzonal-1) c_deltay = (c_lat(c_nzonal/2, c_nmerid) - c_lat(c_nzonal/2, 1))/(c_nmerid-1) else print *,'**** prep_chimere: coarse domain is NOT rectilinear, be patient !' ! Calculate the corners call calculate_corners(c_lon,c_lat,corners_3d,c_nzonal,c_nmerid) endif ! bulk ixx=1 iyy=1 lstix=1 lstiy=1 do ime=1,nc_geog%ny do izo=1,nc_geog%nx ihz=nc_geog%nx*(ime-1)+izo found=.false. if (coarse_is_recti) then lxx = 1 + (nc_geog%lon(izo,ime)-c_lon(1, c_nmerid/2))/c_deltax lyy = 1 + (nc_geog%lat(izo,ime)-c_lat(c_nzonal/2, 1))/c_deltay found = (lxx>=1).and.(lxx<=c_nzonal).and.(lyy>=1).and.(lyy<=c_nmerid) ixx = floor(lxx) iyy = floor(lyy) wx = real(lxx-ixx) wy = real(lyy-iyy) else call findcell(nc_geog%lon(izo,ime),nc_geog%lat(izo,ime),ixx,iyy, & & c_nzonal,c_nmerid,corners_3d,lstix,lstiy,found) if (found) call wcal(nc_geog%lon(izo,ime),nc_geog%lat(izo,ime), & & ixx,iyy,c_nzonal,c_nmerid,corners_3d,wx,wy) endif if (.not. found) then print *,'prep_chimere: fine mesh point outside coarse mesh !' print*,'izo,ime',izo,ime print*,'coordinates of point in fine mesh',nc_geog%lon(izo,ime),nc_geog%lat(izo,ime) print *,'Should never happen ! Bye.' call exit1('Exiting') else ix_b(ihz)=ixx iy_b(ihz)=iyy wx_b(ihz)=wx wy_b(ihz)=wy endif enddo enddo ! boundary ixx=1 iyy=1 lstix=1 lstiy=1 do ihb=1,nhbound found=.false. if (coarse_is_recti) then lxx = 1 + (xb(ihb)-c_lon(1, c_nmerid/2))/c_deltax lyy = 1 + (yb(ihb)-c_lat(c_nzonal/2, 1 ))/c_deltay found = (lxx>=1).and.(lxx<=c_nzonal).and.(lyy>=1).and.(lyy<=c_nmerid) ixx = floor(lxx) iyy = floor(lyy) wx = real(lxx-ixx) wy = real(lyy-iyy) else call findcell(xb(ihb),yb(ihb),ixx,iyy,c_nzonal,c_nmerid,corners_3d, & lstix,lstiy,found) if (found) call wcal(xb(ihb),yb(ihb),ixx,iyy,c_nzonal,c_nmerid,corners_3d,wx,wy) endif if (.not. found) then print *,'prep_chimere: fine mesh boundary point outside coarse mesh !' print*,'coordinates of point in fine mesh',xb(ihb),yb(ihb) print *,'Should never happen ! Bye.' call exit1('Exiting') else ix_l(ihb)=ixx iy_l(ihb)=iyy wx_l(ihb)=wx wy_l(ihb)=wy endif enddo ! Preparing the netCDF output and ini files call prepare_out_file if (len_trim(fnini) > 0) call prepare_ini_file ! find IDs of hlay and airm variables ncstat=nf90_inq_varid(c_id,'hlay',c_hlay_varid) NCERR(__LINE__) ncstat=nf90_inq_varid(c_id,'airm',c_airm_varid) NCERR(__LINE__) ! Check existence of psfc variable ncstat=nf90_inq_varid(c_id,'psfc',c_psfc_varid) exist_psfc = (ncstat == NF90_NOERR) ! write non time-varying data ! lon ncstat=nf90_put_var(f_id,f_lon_varid,nc_geog%lon,start=(/1,1/),count=(/nc_geog%nx,nc_geog%ny/)) NCERR(__LINE__) ! lat ncstat=nf90_put_var(f_id,f_lat_varid,nc_geog%lat,start=(/1,1/),count=(/nc_geog%nx,nc_geog%ny/)) NCERR(__LINE__) ! species. ! This list contains the names of the species without the extension of location (top, north etc.) ! This list contains the species as O3, NO2, p01DUST, p02DUST etc. ncstat=nf90_put_var(f_id,f_species_varid,specname) NCERR(__LINE__) ! Loop over time do ihr=0,nhourrun ncstat=nf90_get_var(c_id,c_times_varid,c_cur_date, & (/1,ihr+1+delta_time_idx/),(/dlen, 1/)) NCERR(__LINE__) if(ihr.eq.0) datebuf=c_cur_date ! read hlay, airm and psfc ncstat=nf90_get_var(c_id,c_hlay_varid,hlay,(/1, 1, 1, ihr+1+delta_time_idx/), & (/c_nzonal, c_nmerid, c_nlevels, 1/)) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_airm_varid,airm,(/1, 1, 1, ihr+1+delta_time_idx/), & (/c_nzonal, c_nmerid, c_nlevels, 1/)) NCERR(__LINE__) if (exist_psfc) then ncstat=nf90_get_var(c_id,c_psfc_varid,psfc,(/1, 1, ihr+1+delta_time_idx/), & (/c_nzonal, c_nmerid, 1/)) NCERR(__LINE__) endif do ny=1,c_nmerid do nx=1,c_nzonal if (exist_psfc) then ps(nx,ny) = psfc(nx,ny) else deltap = cf*airm(nx,ny,1)*hlay(nx,ny,1) ps(nx,ny) = (ac(1)*1e5 + real(deltap))/(1.-bc(1)) end if do ilcoarse=1,c_nlevels pc(nx,ny,ilcoarse) = ac(ilcoarse)*1e5 + bc(ilcoarse)*ps(nx,ny) enddo do ilnest=1,f_nlevels pn(nx,ny,ilnest) = a(ilnest)*1e5 + b(ilnest)*ps(nx,ny) enddo if(pn(nx,ny,f_nlevels).lt.pc(nx,ny,c_nlevels)) then print *,'*** ERROR: Nest domain TOP higher than coarse' call exit1('Exiting') endif enddo enddo ! writing Times ncstat=nf90_put_var(f_id,f_times_varid,c_cur_date,start=(/1,ihr+1/),count=(/dlen, 1/)) NCERR(__LINE__) ! Reading, interpolating, writing concentrations ! we follow the order of "specname" to write concentrations ! to the array. this ensures consistency with the order of ! "species" variable get the id for this variable do ispec=1,nspecies if(ihr.eq.0)print*,'reading and interpolating for ', trim(specname(ispec)), & & ' from ', trim(inpspecname(ispec)), ' bin ',inpbin(ispec) tmpname = specname(ispec) tmpinpname = inpspecname(ispec) ibin = inpbin(ispec) ncstat = nf90_inq_varid(c_id, tmpinpname, c_varid) if(ibin.eq.0)then ncstat=nf90_get_var(c_id,c_varid,xin(:,:,:),(/1, 1, 1, ihr+1+delta_time_idx/), & (/c_nzonal, c_nmerid, c_nlevels, 1/)) else ncstat=nf90_get_var(c_id,c_varid,xin(:,:,:),(/1, 1, 1,ibin , ihr+1+delta_time_idx/), & (/c_nzonal, c_nmerid, c_nlevels, 1, 1/)) endif NCERR(__LINE__) ! units : ppb or ug m-3? ! If ug/m3, we need to convert from [ug m-3] to [ppb]: ! initial and boundary conditions are read in PPB in inichimere ! conc factor: [](ppb)=[](ug/m3)/conv ! airm in [molecules/cm3] must be converted in [ug/m3] ! -> /avogadro*29 to have [g/cm3] ! -> *1.e12 from [g/cm3] to [ug/m3] ! Molar mass in g, conc in ug/m3 ! -> / 1.e9 from fraction to ppb ncstat = nf90_get_att(c_id, c_varid, 'units', units) if (units(1:5) == 'ug/m3') then conv(:,:,:) = airm(:,:,:) * 1.e12 /1.e9 * real(molar(ispec)) / avogadro if (ihr == 0) then print*,'Units for '//trim(tmpname)//' are '//trim(units)// & & '. Converting to ppb with a factor of ', sum(conv(:,:,1))/size(conv(:,:,1)),'...' if (molar(ispec) < 1.e-20) then print *,'* Error The molar mass specified corresponds to "ppb" units?', molar(ispec) call exit1('Check '//trim(fni)//' and '//trim(fnoutspec)) endif endif xin(:,:,:) = xin(:,:,:) / conv(:,:,:) else if (units(1:7) == 'ppb vol') then conv = 1. else print *, '**** Unknown units:'//trim(units)//'!' call exit1('Bye!') endif if (nest_topconc == 1) then ! Read the top_conc from the BOUND nc files for the species and time step coarse_topconc(:,:)=1.d-20 ! Find the current species in appropriate "super-coarse" file tmpinpname = trim(specname(ispec))//'_top' ifle = 0 do i=1,nbn ifle = ifle + 1 ncstat = nf90_inq_varid(cb_id(ifle), trim(tmpinpname), cb_topconc_varid(ifle)) if (ncstat == 0) exit end do if (ncstat == 0) then ! Read coarse_topconc ncstat=nf90_get_var(cb_id(ifle), cb_topconc_varid(ifle), coarse_topconc, & & start=(/1, 1, ihr+1/), count=(/c_nzonal, c_nmerid, 1/) ) NCERR(__LINE__) end if end if ! Vertical interpolation of concentrations if (same_vert_grid) then xin1=xin ! Replace the top boundary of the coarse-mesh run from the top boundary of the global climatology if (nest_topconc == 1) xin1(:,:,c_nlevels) = coarse_topconc else do ny=1,c_nmerid do nx=1,c_nzonal call lintp(pc(nx,ny,:),xin(nx,ny,:), & & c_nlevels,ps(nx,ny),pn(nx,ny,:),xin1(nx,ny,:),f_nlevels) enddo enddo endif ! Horizontal interpolation for top/bulk if(ihr.eq.0) then ideb=1 else ideb=f_nlevels endif call horint(xin1, ideb, f_nlevels, nhoriz, wx_b, wy_b, ix_b, iy_b, xout_bulk) if(ihr.eq.0) bulkbuf(:,:,ispec)=xout_bulk(:,:) ! Horizontal interpolation for lateral boundary call horint(xin1, 1, f_nlevels, nhbound, wx_l, wy_l, ix_l, iy_l, xout_lat) buf_top(ispec,:) = xout_bulk(:,f_nlevels) buf_lat(ispec,:,:) = xout_lat enddo ! ispec = 1,nspecies ! Write results for the BOUND file j=0 do k=1,5 ! Top if(k.eq.1)then do ispec=1,nspecies j=j+1 ncstat=nf90_put_var(f_id,ibvarid(j),buf_top(ispec,:),start=(/1, 1, ihr+1/), & & count=(/nc_geog%nx,nc_geog%ny,1/)) NCERR(__LINE__) enddo endif ! West if(k.eq.2)then do ispec=1,nspecies j=j+1 dim1=1 dim2=nc_geog%ny ncstat=nf90_put_var(f_id,ibvarid(j),buf_lat(ispec,dim1:dim2,:),& & start=(/1, 1, ihr+1/),count=(/nc_geog%ny,f_nlevels,1/)) NCERR(__LINE__) enddo endif ! East if(k.eq.3)then do ispec=1,nspecies j=j+1 dim1=nc_geog%ny+1 dim2=2*nc_geog%ny ncstat=nf90_put_var(f_id,ibvarid(j),buf_lat(ispec,dim1:dim2,:),& & start=(/1, 1, ihr+1/),count=(/nc_geog%ny,f_nlevels,1/)) NCERR(__LINE__) enddo endif ! South if(k.eq.4)then do ispec=1,nspecies j=j+1 dim1=2*nc_geog%ny+1 dim2=2*nc_geog%ny+nc_geog%nx ncstat=nf90_put_var(f_id,ibvarid(j),buf_lat(ispec,dim1:dim2,:),& & start=(/1, 1, ihr+1/),count=(/nc_geog%nx,f_nlevels,1/)) NCERR(__LINE__) enddo endif ! North if(k.eq.5)then do ispec=1,nspecies j=j+1 dim1=2*nc_geog%ny+nc_geog%nx+1 dim2=2*nc_geog%ny+2*nc_geog%nx ncstat=nf90_put_var(f_id,ibvarid(j),buf_lat(ispec,dim1:dim2,:),& & start=(/1, 1, ihr+1/),count=(/nc_geog%nx,f_nlevels,1/)) NCERR(__LINE__) enddo endif enddo ! of k enddo ! ihr = 0,nhourrun ncstat=nf90_close(f_id) ncstat=nf90_close(c_id) ! Save ini data only if fini is not empty if (len_trim(fnini) > 0) then ! INI file ! write data ! LON, LAT ncstat=nf90_put_var(ini_ncid,ini_lon_varid,nc_geog%lon) NCERR(__LINE__) ncstat=nf90_put_var(ini_ncid,ini_lat_varid,nc_geog%lat) NCERR(__LINE__) ! SPECIES ncstat=nf90_put_var(ini_ncid,ini_species_varid,specname) NCERR(__LINE__) ! TIMES ncstat=nf90_put_var(ini_ncid,ini_times_varid,datebuf,(/1,1/),(/dlen,1/)) NCERR(__LINE__) ! CONCENTRATIONS do ispec=1,nspecies do nv=1,nverti ! should be equal to f_nlevels, otherwise segfault with bulkbuf ihz=1 do ime=1,nc_geog%ny do izo=1,nc_geog%nx buf3d(izo,ime,nv)=bulkbuf(ihz,nv,ispec) ihz=ihz+1 enddo enddo enddo ncstat=nf90_put_var(ini_ncid,ini_spec_varid(ispec),buf3d, & (/1,1,1,1/),(/nc_geog%nx,nc_geog%ny,nverti,1/)) NCERR(__LINE__) enddo ncstat=nf90_close(ini_ncid) NCERR(__LINE__) end if ! len_trim(fnini) call nc_close(nc_geog) if (allocated(varname)) deallocate(varname) if (allocated(qqqname)) deallocate(qqqname) if (allocated(specname)) deallocate(specname) if (allocated(molar)) deallocate(molar) if (allocated(wx_b)) deallocate(wx_b) if (allocated(wy_b)) deallocate(wy_b) if (allocated(wx_l)) deallocate(wx_l) if (allocated(wy_l)) deallocate(wy_l) if (allocated(c_lon)) deallocate(c_lon) if (allocated(c_lat)) deallocate(c_lat) if (allocated(xb)) deallocate(xb) if (allocated(yb)) deallocate(yb) if (allocated(ix_b)) deallocate(ix_b) if (allocated(iy_b)) deallocate(iy_b) if (allocated(ix_l)) deallocate(ix_l) if (allocated(iy_l)) deallocate(iy_l) if (allocated(ac)) deallocate(ac) if (allocated(bc)) deallocate(bc) if (allocated(xin)) deallocate(xin) if (allocated(hlay)) deallocate(hlay) if (allocated(airm)) deallocate(airm) if (allocated(psfc)) deallocate(psfc) if (allocated(pc)) deallocate(pc) if (allocated(ps)) deallocate(ps) if (allocated(c_cutoff)) deallocate(c_cutoff) if (allocated(conv)) deallocate(conv) if (allocated(buf_top)) deallocate(buf_top) if (allocated(buf_lat)) deallocate(buf_lat) if (allocated(xout_bulk)) deallocate(xout_bulk) if (allocated(xout_lat)) deallocate(xout_lat) if (nest_topconc == 1) then if (allocated(coarse_topconc)) deallocate(coarse_topconc) if (allocated(bndname)) deallocate(bndname) if (allocated(cb_id)) deallocate(cb_id) if (allocated(cb_species_varid)) deallocate(cb_species_varid) if (allocated(cb_speclen)) deallocate(cb_speclen) if (allocated(cb_species_dimid)) deallocate(cb_species_dimid) if (allocated(cb_topconc_varid)) deallocate(cb_topconc_varid) end if contains !======================================================================================= subroutine read_attributes ncstat=nf90_inquire_attribute(c_id,NF90_GLOBAL,'history',len=c_hlen) NCERR(__LINE__) ncstat=nf90_get_att(c_id,NF90_GLOBAL,'history',c_history) NCERR(__LINE__) c_history=c_history(1:c_hlen) end subroutine read_attributes !---------------------------------------------------------------------------------------- subroutine read_dims ! dimensions IDs ncstat=nf90_inq_dimid(c_id,'Time', c_time_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'DateStrLen', c_date_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'west_east', c_we_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'south_north', c_sn_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'bottom_top', c_bt_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'vcoord_dim', c_vc_dimid) NCERR(__LINE__) ncstat=nf90_inq_dimid(c_id,'ninterv', c_ncutoff_dimid) NCERR(__LINE__) ! lengths ncstat=nf90_inquire_dimension(c_id,c_time_dimid, len=c_times) NCERR(__LINE__) ncstat=nf90_inquire_dimension(c_id,c_date_dimid, len=dl) NCERR(__LINE__) if (dl/=dlen) call exit1('prep_chimere : date format error in coarse-mesh file') ncstat=nf90_inquire_dimension(c_id,c_we_dimid, len=c_nzonal) NCERR(__LINE__) ncstat=nf90_inquire_dimension(c_id,c_sn_dimid, len=c_nmerid) NCERR(__LINE__) ncstat=nf90_inquire_dimension(c_id,c_bt_dimid, len=c_nlevels) NCERR(__LINE__) ncstat=nf90_inquire_dimension(c_id,c_vc_dimid, len=c_vcoord) NCERR(__LINE__) ncstat=nf90_inquire_dimension(c_id,c_ncutoff_dimid, len=c_ncutoff) NCERR(__LINE__) if(c_vcoord.ne.c_nlevels) call exit1('* ERROR in prep_chimere.f90 : nlevels coarse file too low') end subroutine read_dims !---------------------------------------------------------------------------------------- subroutine read_molar_mass ! Assign molar masses corresponding to array specname(nspecies) character(len=16) :: iread real(kind=iprec),dimension(noutspec) :: mmasses integer :: ns, ifnoutspec, io, ispec character(len=16),dimension(noutspec) :: outspecs real(kind=iprec) :: xmass ! read molar masses from OUTPUT_SPECIES call opfi(ifnoutspec,fnoutspec,'f','o') do ispec = 1, noutspec read(ifnoutspec,*,iostat=io) iread, xmass if (io > 0) then call exit1('*** Error reading '//trim(fnoutspec)) else if (io < 0) then print *, 'N output species read: ', ispec exit else outspecs(ispec) = iread mmasses(ispec) = xmass end if end do close(ifnoutspec) ! arrange according to specname(nspecies) molar = -999. do ns = 1, nspecies do ispec = 1, noutspec if (trim(outspecs(ispec)) == trim(specname(ns)) ) then molar(ns) = mmasses(ispec) exit endif enddo if (molar(ns) < -1.) then print *,'**** prep_chimere::read_molar_mass(): molar mass for species ' & &,ns,' : "', trim(specname(ns)), '" could not be assigned!' print *,' Species in '//trim(fni)//':' print *, (trim(specname(ispec))//' ',ispec=1,nspecies) print *,' Species in '//trim(fnoutspec)//':' print *, (trim(outspecs(ispec))//' ',ispec=1,noutspec) call exit1('Bye!') end if end do end subroutine read_molar_mass subroutine read_aerosol integer :: i real :: bidon call opfi(ifn_aerosol,fnaerosol,'f','o') read(ifn_aerosol,*) read(ifn_aerosol,*) n_ncutoff n_ncutoff=n_ncutoff+1 allocate(n_cutoff(n_ncutoff)) read(ifn_aerosol,*) n_cutoff(1:n_ncutoff) read(ifn_aerosol,*) read(ifn_aerosol,*) do i=1,nkc print*,i,nkc read(ifn_aerosol,*) aero_name(i),bidon,bidon,bidon,bidon,bidon enddo close(ifn_aerosol) end subroutine read_aerosol !---------------------------------------------------------------------------------------- subroutine define_bound_coords nbd=0 do ime=1,nc_geog%ny nbd=nbd+1 xb(nbd) = 2*nc_geog%lon(1,ime) - nc_geog%lon(2,ime) yb(nbd) = 2*nc_geog%lat(1,ime) - nc_geog%lat(2,ime) enddo do ime=1,nc_geog%ny nbd=nbd+1 xb(nbd) = 2*nc_geog%lon(nc_geog%nx,ime) - nc_geog%lon(nc_geog%nx-1,ime) yb(nbd) = 2*nc_geog%lat(nc_geog%nx,ime) - nc_geog%lat(nc_geog%nx-1,ime) enddo do izo=1,nc_geog%nx nbd=nbd+1 xb(nbd) = 2*nc_geog%lon(izo,1) - nc_geog%lon(izo,2) yb(nbd) = 2*nc_geog%lat(izo,1) - nc_geog%lat(izo,2) enddo do izo=1,nc_geog%nx nbd=nbd+1 xb(nbd) = 2*nc_geog%lon(izo,nc_geog%ny) - nc_geog%lon(izo,nc_geog%ny-1) yb(nbd) = 2*nc_geog%lat(izo,nc_geog%ny) - nc_geog%lat(izo,nc_geog%ny-1) enddo end subroutine define_bound_coords !---------------------------------------------------------------------------------------- subroutine read_nest_vertical_grid call opfi(ifn_v, fniVCO, 'f','o') ! number of lines in file do nl=1,100000 read(ifn_v,*,end=4755,iostat=ierr)aread,bread if (ierr/=0) then print *, 'prep_chimere : problem reading ', fniVCO call exit1('Exiting') end if enddo 4755 continue f_nlevels = nl - 1 allocate(a(f_nlevels)) allocate(b(f_nlevels)) allocate(pn(c_nzonal,c_nmerid,f_nlevels)) ! real read rewind(ifn_v) do nl=1,f_nlevels read(ifn_v,*)aread,bread a(nl) = real(aread) b(nl) = real(bread) enddo close(ifn_v) end subroutine read_nest_vertical_grid !---------------------------------------------------------------------------------------- subroutine check_time ! checks if time span for fine-mesh run is included in coarse-mesh file ncstat=nf90_inq_varid(c_id,'Times',c_times_varid) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_times_varid,c_start_date,(/1,1/),(/dlen,1/)) NCERR(__LINE__) ncstat=nf90_get_var(c_id,c_times_varid,c_end_date,(/1,c_times/),(/dlen,1/)) NCERR(__LINE__) c_id0=mm5date2numeric(c_start_date) c_ide=mm5date2numeric(c_end_date) if ((id0c_ide)) then print *, 'prep_chimere : coarse-mesh file does not contain required time slots.' call exit1('Exiting.') end if delta_time_idx=interdat(c_id0,id0) end subroutine check_time !---------------------------------------------------------------------------------------- subroutine get_facet(k, endtxt, dim1, dim2) ! Returns netcdf variable extension and dimensions for facet k=1..5 integer, intent(in) :: k character(len=*), intent(out) :: endtxt integer, intent(out) :: dim1, dim2 select case(k) case(1) ! Top endtxt='_top' dim1=f_we_dimid dim2=f_sn_dimid case(2) ! West endtxt='_west' dim1=f_sn_dimid dim2=f_bt_dimid case(3) ! East endtxt='_east' dim1=f_sn_dimid dim2=f_bt_dimid case(4) ! South endtxt='_south' dim1=f_we_dimid dim2=f_bt_dimid case(5) ! North endtxt='_north' dim1=f_we_dimid dim2=f_bt_dimid end select end subroutine get_facet !---------------------------------------------------------------------------------------- subroutine prepare_out_file ! open fine-mesh (output) file integer :: j, k character(len=splen) :: endtxt ! Create boundary conditions file ncstat=nf90_create(fno, cmode, f_id) NCERR(__LINE__) ! create dimensions ncstat=nf90_def_dim(f_id,'Time', NF90_UNLIMITED, f_time_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'DateStrLen', dlen, f_date_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'SpStrLen', splen, f_spstr_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'west_east', nc_geog%nx, f_we_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'south_north',nc_geog%ny, f_sn_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'bottom_top', nverti, f_bt_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'h_boundary', nhbound, f_hb_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(f_id,'Species', nspecies, f_species_dimid) NCERR(__LINE__) ! create variables ncstat=nf90_def_var(f_id,'Times',NF90_CHAR,(/f_date_dimid, f_time_dimid/),f_times_varid) NCERR(__LINE__) ncstat=nf90_def_var(f_id,'species',NF90_CHAR,(/f_spstr_dimid, f_species_dimid/),f_species_varid) NCERR(__LINE__) ncstat=nf90_def_var(f_id,'lon',NF90_DOUBLE,(/f_we_dimid, f_sn_dimid/),f_lon_varid) NCERR(__LINE__) ncstat=nf90_def_var(f_id,'lat',NF90_DOUBLE,(/f_we_dimid, f_sn_dimid/),f_lat_varid) NCERR(__LINE__) ! CONCENTRATION FIELDS j=0 do k=1,5 call get_facet(k, endtxt, dim1, dim2) print*,'Define gaseous and aerosols species:' do ispec=1,nspecies j=j+1 tmpname = trim(specname(ispec))//endtxt print*,ispec,'/',nspecies,' ',trim(tmpname) ncstat=nf90_def_var(f_id,tmpname,NF90_DOUBLE,(/dim1,dim2,f_time_dimid/),ibvarid(j)) NCERR(__LINE__) ncstat=nf90_put_att(f_id,ibvarid(j),'units','ppb') NCERR(__LINE__) ncstat=nf90_put_att(f_id,ibvarid(j),'long_name',trim(tmpname)//' Concentration') NCERR(__LINE__) enddo enddo ! of k ncstat=nf90_put_att(f_id,NF90_GLOBAL,'Title',title) NCERR(__LINE__) ncstat=nf90_put_att(f_id,NF90_GLOBAL,'Sub-title',subtitle) NCERR(__LINE__) ncstat=nf90_put_att(f_id,NF90_GLOBAL,'Generating_process',generator) NCERR(__LINE__) ncstat=nf90_put_att(f_id,NF90_GLOBAL,'Conventions',conventions) NCERR(__LINE__) ncstat=nf90_put_att(f_id,NF90_GLOBAL,'Domain',domain) NCERR(__LINE__) call update_history(c_history,cwd,fno,systime,usrname,hname,fni) ncstat=nf90_put_att(f_id,NF90_GLOBAL,'history', & c_history(1:len_trim(c_history))) NCERR(__LINE__) ! End netCDF definition mode ncstat=nf90_enddef(f_id) NCERR(__LINE__) end subroutine prepare_out_file !---------------------------------------------------------------------------------------- subroutine prepare_ini_file ! create "ini" file ncstat=nf90_create(fnini, cmode, ini_ncid) ; NCERR(__LINE__) ! create dimensions in ini file ncstat=nf90_def_dim(ini_ncid,'Time', NF90_UNLIMITED, ini_time_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'DateStrLen', dlen, ini_date_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'SpStrLen', splen, ini_splen_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'west_east', nc_geog%nx, ini_we_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'south_north', nc_geog%ny, ini_sn_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'bottom_top', nverti, ini_bt_dimid) NCERR(__LINE__) ncstat=nf90_def_dim(ini_ncid,'Species', nspecies, ini_species_dimid) NCERR(__LINE__) ! create variables in ini file ! TIMES ncstat=nf90_def_var(ini_ncid,'Times',NF90_CHAR,(/ini_date_dimid,ini_time_dimid/),ini_times_varid) NCERR(__LINE__) ! LON ncstat=nf90_def_var(ini_ncid,'lon',NF90_DOUBLE,(/ini_we_dimid, ini_sn_dimid/),ini_lon_varid) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_lon_varid,'units','degrees_east') NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_lon_varid,'long_name','Longitude') NCERR(__LINE__) ! LAT ncstat=nf90_def_var(ini_ncid,'lat',NF90_DOUBLE,(/ini_we_dimid, ini_sn_dimid/),ini_lat_varid) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_lat_varid,'units','degrees_north') NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_lat_varid,'long_name','Latitude') NCERR(__LINE__) ! SPECIES ncstat=nf90_def_var(ini_ncid,'species',NF90_CHAR, & & (/ini_splen_dimid, ini_species_dimid/),ini_species_varid) NCERR(__LINE__) ! CONCENTRATION FIELDS do ispec=1,nspecies ncstat=nf90_def_var(ini_ncid,specname(ispec),NF90_DOUBLE, & & (/ini_we_dimid, ini_sn_dimid, ini_bt_dimid, ini_time_dimid/),ini_spec_varid(ispec)) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_spec_varid(ispec),'units','ppb') NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,ini_spec_varid(ispec),'long_name', & & specname(ispec)(1:len_trim(specname(ispec))) // ' Concentration') NCERR(__LINE__) enddo ! create global attributes ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Title',title) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Sub-title',subtitle) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Chimere_type','ini') NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Generating_process',generator) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Conventions',conventions) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'Domain',domain) NCERR(__LINE__) ncstat=nf90_put_att(ini_ncid,NF90_GLOBAL,'history', c_history(1:len_trim(c_history))) NCERR(__LINE__) ncstat=nf90_enddef(ini_ncid) NCERR(__LINE__) end subroutine prepare_ini_file !------------------------------------------------------------------------------------------- end program prep_chimere !------------------------------------------------------------------------------------------- subroutine i2str(i,str) implicit none integer,intent(in) :: i character(len=2),intent(out) :: str write(str,'(I2.2)')i end subroutine i2str !-------------------------------------------------------------------------------------------