module prep_aemis use chimere_common use pncvar use chimere_params, only: iprec implicit none ! fix to bug emis point michele 13/4/2018 #define NCERR(lnum) if (ncstat/=NF90_NOERR) call pnc_err(ncstat,lnum,__FILE__) contains subroutine ini_aemis_params ! Read aemis namelist ! Namelist integer :: fireemis character(len=256) :: fngtm,fno,fnof character(len=*),parameter :: fniARGS = 'prepemis.nml' integer :: ifn_args, ierr namelist /args/ forecast,fngtm,fniVCO,asurfemis,apointemis,fireemis,emission_dir,emission_dir_fire,fno,fnof ! Read input arguments from namelist call opfi(ifn_args,fniARGS,'f','o') read(ifn_args,nml=args,iostat=ierr) if(ierr/=0)call exit1('prepemis: problem reading arguments namelist') close (ifn_args) end subroutine ini_aemis_params subroutine ini_aemis(fname) ! Initialize new AEMISSIONS file character(len=*),intent(in) :: fname integer :: ia_emis, ncstat, ierr character(len=128) :: usrname,hname,systime character(len=256) :: cwd integer,parameter :: datestrlen = 19 ! length of a date string character(len=1024) :: history character(len=*),parameter :: title='CHIMERE SUITE' character(len=*),parameter :: subtitle='Emissions file - Surface + Point sources' character(len=*),parameter :: generator='Generated by prepemis' character(len=*),parameter :: conventions='' character(len=11) :: time11, DateStrLen11, west_east11, south_north11, bottom_top11, SpStrLen11, Species11 type(ncvar_type) :: var_Times, var_lon, var_lat, var_species !------------------------------------------------------------------------------------------ ! Get user/system information call get_system(usrname,hname,systime,cwd) ! open and start to write the AEMISSION.nc file ! The pncvar_create function needs an exact character chain length: here forced to 11 time11='Time' DateStrLen11='DateStrLen' west_east11='west_east' south_north11='south_north' bottom_top11='bottom_top' SpStrLen11='SpStrLen' Species11='Species' call pncvar_create(emisa_nc, fname, & & (/time11, DateStrLen11, west_east11, south_north11, bottom_top11, SpStrLen11, Species11/), & & (/int(NF90_UNLIMITED), dlen, nzonal_entire_domain, nmerid_entire_domain, nverti, splen, nemisa /), wrk_comm) ! Coord vars call pncvar_create_var(emisa_nc, 'Times', emisa_nc%did((/2,1/)), NF90_CHAR, var_Times, 'Date string', 'WRF-style Times') call pncvar_create_var(emisa_nc, 'species', emisa_nc%did((/6,7/)), NF90_CHAR, var_species, '-', 'Species names') call pncvar_create_var(emisa_nc, 'lon', emisa_nc%did((/3,4/)), NF90_DOUBLE, var_lon, 'degrees_east', 'Longitude') call pncvar_create_var(emisa_nc, 'lat', emisa_nc%did((/3,4/)), NF90_DOUBLE, var_lat, 'degrees_north', 'Latitude') ! Emissions do ia_emis=1,nemisa call pncvar_create_var(emisa_nc, emisa_species(ia_emis)%name, emisa_nc%did((/3,4,5,1/)), NF90_DOUBLE, var_emisa(ia_emis), & 'molecule/cm2/s', trim(emisa_species(ia_emis)%name)//' Emission') enddo ! Global attrs ncstat=nf90mpi_put_att(emisa_nc%id, NF90_GLOBAL, 'Title', title) NCERR(__LINE__) ncstat=nf90mpi_put_att(emisa_nc%id,NF90_GLOBAL,'Sub-title',subtitle) NCERR(__LINE__) ncstat=nf90mpi_put_att(emisa_nc%id,NF90_GLOBAL,'Generating_process',generator) NCERR(__LINE__) ncstat=nf90mpi_put_att(emisa_nc%id,NF90_GLOBAL,'Conventions',conventions) NCERR(__LINE__) ncstat=nf90mpi_put_att(emisa_nc%id,NF90_GLOBAL,'Domain',domain) NCERR(__LINE__) #if (defined(IFORT)+defined(G95)+defined(PGI)+defined(GFORTRAN)) if (rank == 0) then history = & 'File '//trim(fnemisa) // & ' was generated on ' // & trim(systime) // & ' by ' // & trim(usrname) // & ' on ' // & trim(hname) end if call mpi_bcast(history, 1024, mpi_character, 0, wrk_comm, ierr) #else history='unknown' #endif ncstat = nf90mpi_put_att(emisa_nc%id, NF90_GLOBAL, 'history', trim(history)) NCERR(__LINE__) ncstat = nf90mpi_enddef(emisa_nc%id) NCERR(__LINE__) ! Write lon/lat call pncvar_put_var(domains, var_lon, xlong(1:nzonal,1:nmerid), rank, (/nzonal,nmerid/)) call pncvar_put_var(domains, var_lat, xlati(1:nzonal,1:nmerid), rank, (/nzonal,nmerid/)) ! Write species call pncvar_write_char_array(emisa_nc, var_species, rank, emisa_species(1:nemisa)%name) !-------------------------------------------------- end subroutine ini_aemis !-------------------------------------------------- subroutine update_aemis(ihourrun, ksens) ! If needed open AEMISSIONS file & init vars integer, intent(in) :: ihourrun, ksens integer :: ne, ncstat, delta_time character(len=512) :: file_name ! External function character(len=512) :: get_out_file ! hours since idatestart in the current record delta_time = ihourrun + ksens - 1 ! Exit unless at the beginning of AEMISSIONS file if (mod(delta_time, frames_per_file) /= 0) return ! Close old .nc files if (delta_time > 0) then ncstat = nf90mpi_close(emisa_nc%id) NCERR(__LINE__) end if ! Name of the new file to open file_name = get_out_file(fnemisa, idate(delta_time), frames_per_file, idateend, '') ! Open new file call pncvar_open(emisa_nc, trim(file_name), NF90_NOWRITE, wrk_comm) ! Check dim sizes if (rank == 0) then if (findlen(emisa_nc,'west_east') /= nzonal_entire_domain) & & call exit1( 'update_aemis : west_east dimension inconsistency') if (findlen(emisa_nc,'south_north') /= nmerid_entire_domain) & & call exit1( 'update_aemis : south_north dimension inconsistency') if (findlen(emisa_nc,'SpStrLen') /= splen) & & call exit1( 'update_aemis : species name format error in A-emissions file') if (findlen(emisa_nc,'DateStrLen') /= dlen) & & call exit1( 'update_aemis : date format error in A-emissions file') if (findlen(emisa_nc,'Species') < nemisa) & & call exit1( 'update_aemis : Species dimension inconsistency') end if ! Get the var IDs of each species in AEMISSIONS file do ne=1,nemisa call pncvar_get_var_info(emisa_nc, emisa_species(ne)%name, var_emisa(ne)) enddo end subroutine update_aemis !-------------------------------------------------- subroutine read_aemissions(ihourrun, ksens) integer, intent(in) :: ihourrun, ksens integer, save :: irec = -9999 integer :: ne, ip, i, ik, nl, ime, izo, delta_time ! hours since idatestart in the current record delta_time = ihourrun + ksens - 1 ! Reset counter when new file if (mod(delta_time, frames_per_file) == 0) irec = 0 irec = irec + 1 ! Read current time & check consistency with ihourrun/idatestart call pncvar_get_check_time(emisa_nc, rank, idatestart, delta_time, irec) emisa(:,:,:,:,ksens) = dzero ! Gas emissions for each active anthropic species do ne=1,nemisa call pncvar_get_var(domains, var_emisa(ne), emisa(ne,1:nzonal,1:nmerid,1:nlevemis,ksens), rank, & & (/nzonal,nmerid,nlevemis,1/), (/irec/)) enddo ! Particulate emissions if (bins == 1) then emipa(:,:,:,:,:,ksens) = dzero do ne=1,nemisa ip = iprofemipa(ne) if(ip.ne.0) then ik = ispecemip(ip) do i=1,ms do ime=1,nmerid do izo=1,nzonal do nl=1,nlevemis emipa(ik,i,izo,ime,nl,ksens) = emipa(ik,i,izo,ime,nl,ksens) & & + emisa(ne,izo,ime,nl,ksens)*emiprof(ip,i) enddo enddo enddo enddo endif enddo endif end subroutine read_aemissions !-------------------------------------------------- subroutine process_aemis(ks,itimeslot) integer,intent(in) :: ks,itimeslot integer :: ia_emis integer :: iemis_point(nemisa),iemis_surf(nemisa) real(kind=iprec),allocatable,dimension(:,:,:) :: buf3d character(len=256) :: emisa_surf_file,emisa_point_file character(len=2) :: month integer :: ncstat,s_ncid,ime,izo,iv,i,ik,ip type(ncvar_type),dimension(:),allocatable :: var_emisa_tab integer, save :: irec character(len=512) :: get_out_file emisa(:,:,:,:,ks) = dzero ! Allocation of pncvar variables array if (asurfemis==1.and..not. allocated(var_surf_emisa_spec)) allocate(var_surf_emisa_spec(nemisa,0:nhourrun)) if (asurfemis==1.and..not. allocated(var_surf_emisa_emepheight)) allocate(var_surf_emisa_emepheight(nemisa,0:nhourrun)) if (asurfemis==1.and..not. allocated(emisa_surf_out)) allocate(emisa_surf_out(nemisa,0:nhourrun)) if (apointemis==1.and..not. allocated(var_emisa_point_sourceilat_out)) allocate(var_emisa_point_sourceilat_out(nemisa,0:nhourrun)) if (apointemis==1.and..not. allocated(var_emisa_point_sourceilon_out)) allocate(var_emisa_point_sourceilon_out(nemisa,0:nhourrun)) if (apointemis==1.and..not. allocated(var_emisa_point_sourceheight_out)) allocate(var_emisa_point_sourceheight_out(nemisa,0:nhourrun)) if (apointemis==1.and..not. allocated(var_emisa_point_spec_out)) allocate(var_emisa_point_spec_out(nemisa,0:nhourrun)) if (apointemis==1.and..not. allocated(emisa_point_out)) allocate(emisa_point_out(nemisa,0:nhourrun)) if (imonth(itimeslot) < 10) then write(month,'(I1)'), imonth(itimeslot) month = "0"//month else write(month,'(I2)'), imonth(itimeslot) end if ! Output to AEMISSIONS..nc if (mod(itimeslot, frames_per_file) == 0) then ! Close old file if opened if (itimeslot > 0) then ncstat=nf90mpi_close(emisa_nc%id) NCERR(__LINE__) endif ! Create new file & variables call ini_aemis(trim(get_out_file(fnemisa, idate(itimeslot), frames_per_file, idateend, simlab))) ! reset counter irec = 0 endif irec = irec + 1 ! Save time (starting from idatestart) call pncvar_write_time(emisa_nc, rank, idatestart, itimeslot, irec) allocate(buf3d(nzonal,nmerid,nverti)) do ia_emis=1,nemisa buf3d=dzero ! Read surface emission file if(asurfemis.eq.1) then ! File to be read emisa_surf_file=emission_dir(1:len_trim(emission_dir))//"/EMIS."//domain(1:len_trim(domain))//"."//month//"."//emisa_species(ia_emis)%name(1:len_trim(emisa_species(ia_emis)%name))//".s.nc" iemis_surf(ia_emis)=1 if (trim(emisa_surf_file) /= trim(emisa_surf_file_prev(ia_emis)).or.itimeslot==0) then ! Check if file exists ncstat = nf90mpi_open(wrk_comm, trim(emisa_surf_file), NF90_NOWRITE, MPI_INFO_NULL, s_ncid) if (ncstat == nf90_noerr) then ncstat=nf90mpi_close(s_ncid) NCERR(__LINE__) ! Note that for now only 2 variables are read from the surface emission files and that the order of these variable is hard coded. allocate(var_emisa_tab(2)) call read_surface_emission(emisa_surf_file,emisa_species(ia_emis)%name,nzonal_entire_domain,nmerid_entire_domain,wrk_comm,var_emisa_tab,emisa_surf_out(ia_emis,itimeslot)) var_surf_emisa_spec(ia_emis,itimeslot)=var_emisa_tab(1) var_surf_emisa_emepheight(ia_emis,itimeslot)=var_emisa_tab(2) deallocate(var_emisa_tab) else iemis_surf(ia_emis)=0 if (rank==0 .and. chimverb > 3) print *,"No surface emissions for : ",emisa_species(ia_emis)%name(1:len_trim(emisa_species(ia_emis)%name))," Month : ", month end if else var_surf_emisa_spec(ia_emis,itimeslot) = var_surf_emisa_spec(ia_emis,itimeslot-1) var_surf_emisa_emepheight(ia_emis,itimeslot) = var_surf_emisa_emepheight(ia_emis,itimeslot-1) emisa_surf_out(ia_emis,itimeslot) = emisa_surf_out(ia_emis,itimeslot-1) iemis_surf(ia_emis) = iemis_surf_prev(ia_emis) end if end if ! Read point emission file if(apointemis.eq.1) then ! File to be read emisa_point_file=emission_dir(1:len_trim(emission_dir))//"/EMIS."//domain(1:len_trim(domain))//"."//month//"."//emisa_species(ia_emis)%name(1:len_trim(emisa_species(ia_emis)%name))//".p.nc" iemis_point(ia_emis)=1 if (trim(emisa_point_file) /= trim(emisa_point_file_prev(ia_emis)).or.itimeslot==0) then ! Check if file exists ncstat = nf90mpi_open(wrk_comm, trim(emisa_point_file), NF90_NOWRITE, MPI_INFO_NULL, s_ncid) if (ncstat == nf90_noerr) then ncstat=nf90mpi_close(s_ncid) NCERR(__LINE__) ! Note that for now only 4 variables are read from the point emission files and that the order of these variable is hard coded. allocate(var_emisa_tab(4)) call read_point_emission(emisa_point_file,emisa_species(ia_emis)%name,nzonal_entire_domain,nmerid_entire_domain,wrk_comm,var_emisa_tab,emisa_point_out(ia_emis,itimeslot)) var_emisa_point_sourceilat_out(ia_emis,itimeslot) = var_emisa_tab(1) var_emisa_point_sourceilon_out(ia_emis,itimeslot) = var_emisa_tab(2) var_emisa_point_sourceheight_out(ia_emis,itimeslot) = var_emisa_tab(3) var_emisa_point_spec_out(ia_emis,itimeslot) = var_emisa_tab(4) deallocate(var_emisa_tab) else iemis_point(ia_emis)=0 end if else var_emisa_point_sourceilat_out(ia_emis,itimeslot) = var_emisa_point_sourceilat_out(ia_emis,itimeslot-1) var_emisa_point_sourceilon_out(ia_emis,itimeslot) = var_emisa_point_sourceilon_out(ia_emis,itimeslot-1) var_emisa_point_sourceheight_out(ia_emis,itimeslot) = var_emisa_point_sourceheight_out(ia_emis,itimeslot-1) var_emisa_point_spec_out(ia_emis,itimeslot) = var_emisa_point_spec_out(ia_emis,itimeslot-1) emisa_point_out(ia_emis,itimeslot)=emisa_point_out(ia_emis,itimeslot-1) iemis_point(ia_emis) = iemis_point_prev(ia_emis) end if end if ! Process surface and point emissions if(asurfemis.eq.1.and.iemis_surf(ia_emis).eq.1) call process_surface_emission(ia_emis,itimeslot,ks,buf3d) if(apointemis.eq.1.and.iemis_point(ia_emis).eq.1) call process_point_emission(ia_emis,itimeslot,ks,buf3d) ! Write anthropogenic emission file call pncvar_put_var(domains,var_emisa(ia_emis), buf3d(1:nzonal,1:nmerid,1:nverti), rank, (/nzonal,nmerid,nverti,1/), (/irec/)) ! Gas emissions for each active anthropic species emisa(ia_emis,1:nzonal,1:nmerid,1:nverti,ks)=buf3d(1:nzonal,1:nmerid,1:nverti) emisa_surf_file_prev(ia_emis) = emisa_surf_file emisa_point_file_prev(ia_emis) = emisa_point_file enddo iemis_surf_prev = iemis_surf iemis_point_prev = iemis_point ! Particulate emissions if(bins.eq.1)then do ia_emis=1,nemisa ip = iprofemipa(ia_emis) if(ip.ne.0) then ik = ispecemip(ip) do i=1,ms do ime=1,nmerid do izo=1,nzonal do iv=1,nverti emipa(ik,i,izo,ime,iv,ks) = emipa(ik,i,izo,ime,iv,ks) & & + emisa(ia_emis,izo,ime,iv,ks)*emiprof(ip,i) enddo enddo enddo enddo endif enddo endif ! Synchronize disk to avoid data loss ncstat=nf90mpi_sync(emisa_nc%id) NCERR(__LINE__) end subroutine process_aemis subroutine read_surface_emission(emis_file,species_name,nzonal_entire_domain,nmerid_entire_domain,communicator,var_emisa_out,emisa_surf_out) ! This routine creates and fills an ncfile structure for a given ! input file name. It does not use the chimere_common module to ! simplify the adaptation to others input file formats. ! Input character(len=256),intent(in) :: emis_file integer,intent(in) :: nzonal_entire_domain,nmerid_entire_domain integer,intent(in) :: communicator ! used for parallel netcdf reading character(len=*),intent(in) :: species_name ! Output (order of variable in the table is hard coded) type(ncvar_type),intent(out),dimension(2) :: var_emisa_out type(ncfile),intent(out) :: emisa_surf_out call pncvar_open(emisa_surf_out,emis_file,NF90_NOWRITE,communicator) if(findlen(emisa_surf_out,'Time')/=24) call exit1('prepemis : NHOURPDAY dimension inconsistency. Exiting') if(findlen(emisa_surf_out,'west_east')/=nzonal_entire_domain) call exit1('prepemis : WE dimension inconsistency. Exiting') if(findlen(emisa_surf_out,'south_north')/=nmerid_entire_domain) call exit1('prepemis : SN dimension inconsistency. Exiting') if(findlen(emisa_surf_out,'type_day')/=7) call exit1('prepemis : TYPE_DAY dimension inconsistency. Exiting') call pncvar_get_var_info(emisa_surf_out, trim(species_name), var_emisa_out(1)) call pncvar_get_var_info(emisa_surf_out, 'EMEP_levels', var_emisa_out(2)) end subroutine read_surface_emission subroutine read_point_emission(emis_file,species_name,nzonal_entire_domain,nmerid_entire_domain,communicator,var_emisa_out,emisa_point_out) implicit none ! This routine creates and fills an ncfile structure for a given ! input file name. It does not use the chimere_common module to ! simplify the adaptation to others input file formats. ! Input character(len=256),intent(in) :: emis_file integer,intent(in) :: nzonal_entire_domain,nmerid_entire_domain integer,intent(in) :: communicator ! used for parallel netcdf reading character(len=*),intent(in) :: species_name ! Output (order of variable in the table is hard coded) type(ncvar_type),intent(out),dimension(4) :: var_emisa_out type(ncfile),intent(out) :: emisa_point_out call pncvar_open(emisa_point_out,emis_file,NF90_NOWRITE,communicator) if(findlen(emisa_point_out,'Time')/=24) call exit1('prepemis : NHOURPDAY dimension inconsistency. Exiting') if(findlen(emisa_point_out,'west_east')/=nzonal_entire_domain) call exit1('prepemis : WE dimension inconsistency. Exiting') if(findlen(emisa_point_out,'south_north')/=nmerid_entire_domain) call exit1('prepemis : SN dimension inconsistency. Exiting') if(findlen(emisa_point_out,'type_day')/=7) call exit1('prepemis : TYPE_DAY dimension inconsistency. Exiting') call pncvar_get_var_info(emisa_point_out,'Sources_ilat',var_emisa_out(1)) call pncvar_get_var_info(emisa_point_out,'Sources_ilon',var_emisa_out(2)) call pncvar_get_var_info(emisa_point_out,'Sources_Height',var_emisa_out(3)) call pncvar_get_var_info(emisa_point_out, trim(species_name),var_emisa_out(4)) end subroutine read_point_emission !========================================================================= subroutine process_surface_emission(ia_emis,itimeslot,ks,buf3d) ! Inputs integer,intent(in) :: ia_emis ! index aemis integer,intent(in) :: itimeslot ! index time slot integer,intent(in) :: ks ! Output real(kind=iprec),dimension(nzonal,nmerid,nverti),intent(out) :: buf3d real(kind=iprec),allocatable,dimension(:,:) :: dis1todis2 real(kind=iprec),allocatable,dimension(:) :: hlayext,emepheight integer :: iv,izo,ime real(kind=iprec),dimension(:,:,:),allocatable :: buf3d_tmp integer :: nlevel_emep,idayofweek call dayofweek(idate(itimeslot),idayofweek) nlevel_emep=findlen(emisa_surf_out(ia_emis,itimeslot),'bottom_top') allocate(buf3d_tmp(nzonal,nmerid,nlevel_emep)) call pncvar_get_var(domains,var_surf_emisa_spec(ia_emis,itimeslot), buf3d_tmp(1:nzonal,1:nmerid,1:nlevel_emep), rank,(/nzonal,nmerid,nlevel_emep,1,1/),(/1,idayofweek,ihour(itimeslot)+1/)) ! Get emep height allocate(emepheight(0:nlevel_emep)) emepheight=0.0 call pncvar_get_var(domains,var_surf_emisa_emepheight(ia_emis,itimeslot),emepheight(1:nlevel_emep),rank,(/nlevel_emep/),(/1/)) allocate(hlayext(0:nverti)) allocate(dis1todis2(nlevel_emep,nverti)) do ime=1,nmerid do izo=1,nzonal do iv=1,nverti hlayext(iv)=hlay(izo,ime,iv,ks) enddo hlayext(0)=0. call modisec_alt(emepheight(0:nlevel_emep),nlevel_emep+1,hlayext,nverti+1,dis1todis2) buf3d(izo,ime,:)=matmul(transpose(dis1todis2(:,:)),buf3d_tmp(izo,ime,:)) enddo enddo end subroutine process_surface_emission subroutine modisec_alt(dist1,nsec1,dist2,nsec2,dis1todis2) ! modisec w/o logarithmic ponderation, which is adapted only for aerosols integer nsec1,nsec2 integer loca(nsec2) real(kind=iprec) :: dist1(nsec1),dis1todis2(nsec1-1,nsec2-1) real(kind=iprec) :: dist2(nsec2) integer :: i,j dis1todis2=0.0 if(dist2(1).lt.dist1(1)) loca(1)=0 do i=1,nsec2-1 if(dist2(i+1).lt.dist1(1)) loca(i+1)=0 do j=1,nsec1-1 if(dist2(i).lt.dist1(j+1).and.dist2(i).ge.dist1(j))loca(i)=j enddo if(dist2(i).ge.dist1(nsec1)) loca(i)=nsec1 enddo if(dist2(nsec2).gt.dist1(nsec1)) loca(nsec2)=nsec1 if(dist2(nsec2).lt.dist1(nsec1)) loca(nsec2)=nsec1-1 do j=1,nsec1-1 do i=1,nsec2-1 if(loca(i).eq.j.and.loca(i+1).eq.j) then dis1todis2(j,i)=((dist2(i+1))-(dist2(i)))/((dist1(j+1))-(dist1(j))) elseif(loca(i).eq.j.and.loca(i+1).gt.j) then dis1todis2(j,i)=((dist1(j+1))-(dist2(i)))/((dist1(j+1))-(dist1(j))) elseif(loca(i).lt.j.and.loca(i+1).eq.j) then dis1todis2(j,i)=((dist2(i+1))-(dist1(j)))/((dist1(j+1))-(dist1(j))) elseif(loca(i).lt.j.and.loca(i+1).gt.j) then dis1todis2(j,i)=1.0 endif enddo enddo end subroutine modisec_alt subroutine process_point_emission(ia_emis,itimeslot,ks,buf3d) ! Inputs integer,intent(in) :: ia_emis ! index aemis integer,intent(in) :: itimeslot ! index time slot integer,intent(in) :: ks ! Output real(kind=iprec),dimension(nzonal,nmerid,nverti),intent(out) :: buf3d real(kind=iprec),dimension(:),allocatable :: buf1d,height_sources integer :: idayofweek,p_nsources,ivert,ns_p integer :: ilon_chim,ilat_chim ! ilon and ilat on CHIMERE subdomain logical :: found real(kind=iprec) :: height1,height2 integer,dimension(:),allocatable :: ivert_sources,ilon_sources,ilat_sources call dayofweek(idate(itimeslot),idayofweek) p_nsources=findlen(emisa_point_out(ia_emis,itimeslot),'Sources') ! Allocatations allocate(buf1d(p_nsources),height_sources(p_nsources),ivert_sources(p_nsources),ilon_sources(p_nsources),ilat_sources(p_nsources)) !original error point call pncvar_get_var(domains,var_surf_emisa_spec(ia_emis,itimeslot), buf1d(1:p_nsources), rank,(/p_nsources,1,1/),(/1,idayofweek,ihour(itimeslot)+1/)) ! fix to bug point error michele 13/4/2018 call pncvar_get_var(domains,var_emisa_point_spec_out(ia_emis,itimeslot), buf1d(1:p_nsources), rank,(/p_nsources,1,1/),(/1,idayofweek,ihour(itimeslot)+1/)) ! ! Read sources height call pncvar_get_var(domains,var_emisa_point_sourceheight_out(ia_emis,itimeslot), height_sources(1:p_nsources), rank,(/p_nsources/),(/1/)) ! Read sources ilon call pncvar_get_var_int(domains,var_emisa_point_sourceilon_out(ia_emis,itimeslot), ilon_sources(1:p_nsources), rank,(/p_nsources/),(/1/)) ! Read sources ilat call pncvar_get_var_int(domains,var_emisa_point_sourceilat_out(ia_emis,itimeslot), ilat_sources(1:p_nsources), rank,(/p_nsources/),(/1/)) do ns_p=1,p_nsources ilon_chim = ilon_sources(ns_p) - domains%dom(rank+1)%izstart + 1 ilat_chim = ilat_sources(ns_p) - domains%dom(rank+1)%imstart + 1 !fix to ensure point emission inside domain michele 13/4/2018 if (ilon_chim > 0.and.ilon_chim <= domains%dom(rank+1)%nzcount.and.ilat_chim > 0.and. & & ilat_chim <= domains%dom(rank+1)%nmcount) then ! found=.false. do ivert=1,nverti if (.not.found) then if(ivert.eq.1)then height1=0. else height1=hlay(ilon_chim,ilat_chim,ivert-1,ks) end if height2=hlay(ilon_chim,ilat_chim,ivert,ks) if(height_sources(ns_p).ge.height1.and.height_sources(ns_p).lt.height2) then ivert_sources(ns_p)=ivert found=.true. endif end if enddo if(.not.found) call exit1('*** ERROR in '//__FILE__//' **** Check your emission source point') enddo do ns_p=1,p_nsources ilon_chim = ilon_sources(ns_p) - domains%dom(rank+1)%izstart + 1 ilat_chim = ilat_sources(ns_p) - domains%dom(rank+1)%imstart + 1 !fix to ensure that source is located in the current subdomain michele 13/4/2018 if (ilon_chim > 0.and.ilon_chim <= domains%dom(rank+1)%nzcount.and.ilat_chim > 0.and. & & ilat_chim <= domains%dom(rank+1)%nmcount) then buf3d(ilon_chim,ilat_chim,ivert_sources(ns_p))=buf3d(ilon_chim,ilat_chim,ivert_sources(ns_p))+buf1d(ns_p) end if enddo end subroutine process_point_emission end module prep_aemis