!###############################################################################! ! CCATT-BRAMS/MCGA-CPTEC/WRF emission model CPTEC/INPE ! ! Version 1.0.0: 12/nov/2010 ! ! Coded by Saulo Freitas and Karla Longo ! ! Contact: gmai@cptec.inpe.br - http://meioambiente.cptec.inpe.br ! !###############################################################################! module volc_degassing_emissions ! ------------------------------------------------ integer, parameter :: nvolcanoes=10000 integer, parameter :: maxnspecies = 200 , nspecies = 1 integer, parameter :: & SO2 = 1 character(LEN=25),dimension(nspecies),parameter :: spc_name= & ! '1234567890123456789012345' (/ & 'SO2 '& /) !--- type degassing_vars real, pointer, dimension(:,:,:) :: src end type degassing_vars type (degassing_vars), allocatable :: degassing_g(:) !--- type volc_vars real, pointer, dimension(:) :: volcano end type volc_vars type (volc_vars), allocatable :: volc_prop(:) integer, parameter :: nprop = 9 integer, parameter :: & LATI =1 & ,LONG =2 & ,ELEV =3 & ,BEGT =4 & ,CRAT =5 & ,HEAT =6 & ,plum_heigth =7 & ,DURA =8 & ,ESO2 =9 !--- type volc_vars_P real, pointer, dimension(:,:) :: prop end type volc_vars_P type (volc_vars_P), allocatable :: degassingP_g(:) contains !--------------------------------------------------------------- subroutine alloc_volcanoes(degassing_g,degassingP_g,n1,n2,n3) implicit none type (degassing_vars),dimension(nspecies) :: degassing_g type (volc_vars_P),dimension(nprop) :: degassingP_g integer,intent(in) :: n1,n2,n3 integer ispc,iprop do ispc=1,nspecies allocate (degassing_g(ispc)%src(n1,n2,n3)) degassing_g(ispc)%src(:,:,:)=0. enddo do iprop=1,nprop allocate (degassingP_g(iprop)%prop(n1,n2)) degassingP_g(iprop)%prop(:,:)=0. enddo end subroutine alloc_volcanoes !--------------------------------------------------------------- subroutine nullify_volcanoes(degassing_g,degassingP_g) implicit none type (degassing_vars),dimension(nspecies) :: degassing_g type (volc_vars_P) ,dimension(nprop) :: degassingP_g integer ispc,iprop do ispc=1,nspecies if (associated(degassing_g(ispc)%src)) nullify (degassing_g(ispc)%src) enddo do iprop=1,nprop if (associated(degassingP_g(iprop)%prop)) nullify (degassingP_g(iprop)%prop) enddo end subroutine nullify_volcanoes !--------------------------------------------------------------- subroutine alloc_volc_prop(volc_prop) implicit none type (volc_vars),dimension(nprop) :: volc_prop integer i do i=1,nprop allocate (volc_prop(i)%volcano(nvolcanoes)) enddo end subroutine alloc_volc_prop !--------------------------------------------------------------- subroutine nullify_volc_prop(volc_prop) implicit none type (volc_vars),dimension(nprop) :: volc_prop integer i do i=1,nprop if (associated(volc_prop(i)%volcano)) nullify (volc_prop(i)%volcano) enddo end subroutine nullify_volc_prop !--------------------------------------------------------------- subroutine mem_volc_prop() implicit none integer i if(.not. allocated(volc_prop)) allocate(volc_prop(nprop)) call nullify_volc_prop(volc_prop(:)) call alloc_volc_prop (volc_prop(:)) end subroutine mem_volc_prop !--------------------------------------------------------------- subroutine get_dimension_per_file_volcano(volc_dir,iyear,imon,iday& ,nevents) use netcdf implicit none integer, intent(in) ::iyear,imon,iday character(len=*), intent(in) ::volc_dir integer, intent(out) ::nevents character (len=240) :: filename integer ncid integer current_date,neventsID character(len=4) :: cyear logical :: there nevents=0 current_date=iyear*10000+imon*100+iday write(cyear(1:4),10) iyear 10 format (i4.4) filename=trim(volc_dir)//'/volc_so2_'//cyear//'.nc' inquire(file=trim(filename),exist=there) if (.not. there) then print*,'==============================================================' print*,'Volcanic data not found', TRIM(filename) print*,'null data will set in the emission files' print*,'==============================================================' return endif print*,' ============= volcanoes =====================================' print*,' opening degassing volc data:',current_date,trim(filename) ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. call check(NF90_OPEN(filename, NF90_NOWRITE, ncid)) call check(nf90_inq_dimid(ncid, "nevents", neventsID)) call check(nf90_inquire_dimension(ncid, neventsID, len = nevents)) print*,'nevents found:', nevents end subroutine get_dimension_per_file_volcano !--------------------------------------------------------------- subroutine get_emission_per_volcano(volc_dir,actual_number_volc,iyear,imon,iday& ,qsc,qlong,qlati,qelev,qBEGT,qcrat,qheat,qplum_heigth & ,qdura,QESO2,nevents) use netcdf implicit none !include 'netcdf.inc' integer, intent(in) ::iyear,imon,iday,nevents character(len=*), intent(in) ::volc_dir integer, intent(out) ::actual_number_volc real, dimension(nvolcanoes) , intent(out):: qlong,qlati,qelev,qBEGT,qcrat,qheat,qplum_heigth & ,qdura ,QESO2 real, dimension(nvolcanoes,nspecies), intent(out) :: qsc character (len=240) :: filename integer i1,i2,iesp,ifoc,ncid,var_id integer,allocatable, dimension (:) :: date real, allocatable, dimension (:) :: cloud_column_height,so2,lon,lat,elevation integer nd1,nd2,current_date character(len=4) :: cyear logical :: there allocate (cloud_column_height(nevents),so2(nevents),& lon(nevents),lat(nevents),elevation(nevents),date(nevents)) cloud_column_height=0.;so2=0.;lon=0.;lat=0.;elevation=0.;date=0 current_date=iyear*10000+imon*100+iday write(cyear(1:4),10) iyear 10 format (i4.4) filename=trim(volc_dir)//'/volc_so2_'//cyear//'.nc' inquire(file=trim(filename),exist=there) if (.not. there) then print*,'==============================================================' print*,'Volcanic data not found', TRIM(filename) print*,'null data will set in the emission files' print*,'==============================================================' return endif print*,' ============= volcanoes =====================================' print*,' opening degassing volc data:',current_date,trim(filename) ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. call check(NF90_OPEN(filename, NF90_NOWRITE, ncid)) ! Get the varid of the data variable, based on its name. call check( nf90_inq_varid(ncid, 'date', var_id) ) call check( nf90_get_var(ncid, var_id, date ) ) !- get the number of volcanos emitting on the current day do i1=1,nevents if(date(i1) == current_date) then nd1=i1;exit endif enddo do i1=nevents,1,-1 if(date(i1) == current_date) then nd2=i1;exit endif enddo actual_number_volc=nd2-nd1+1 if(actual_number_volc>nvolcanoes) stop 'actual_number_volc>nvolcanoes' print*,'Amount of volcanoes found:', actual_number_volc,date(nd1),date(nd2) !- get data call check( nf90_inq_varid(ncid, 'so2', var_id) ) call check( nf90_get_var(ncid, var_id, so2 ) ) call check( nf90_inq_varid(ncid, 'lat', var_id) ) call check( nf90_get_var(ncid, var_id, lat ) ) call check( nf90_inq_varid(ncid, 'lon', var_id) ) call check( nf90_get_var(ncid, var_id, lon ) ) call check( nf90_inq_varid(ncid, 'elevation', var_id) ) !meters call check( nf90_get_var(ncid, var_id, elevation ) )!meters call check( nf90_inq_varid(ncid, 'cloud_column_height', var_id) )!meters call check( nf90_get_var(ncid, var_id, cloud_column_height ) ) !meters !- calculates the emission by volcano do ifoc = 1,actual_number_volc QLONG(ifoc)=lon (nd1+ifoc-1) QLATI(ifoc)=lat (nd1+ifoc-1) QELEV(ifoc)=elevation(nd1+ifoc-1)!meters asl QESO2(ifoc)=so2 (nd1+ifoc-1)*1.e+6 !from kiloton to kg Qplum_heigth(ifoc)=cloud_column_height(nd1+ifoc-1)!meters asl !not defined QBEGT(ifoc)=-999. QCRAT(ifoc)=-999. QHEAT(ifoc)=-999. QDURA(ifoc)=86400. ! emission per day !not defined if(ifoc==1) print*,'Volcano --- LONG --- LATI ---- ELEV(m) --- plum_heigth(m) --- ESO2' write(*,100) ifoc ,QLONG(ifoc), QLATI(ifoc), QELEV(ifoc), Qplum_heigth(ifoc), QESO2(ifoc) 100 format(1x,i6,5F13.2) do iesp=1,nspecies qsc(ifoc,iesp) = QESO2(ifoc) !already in kg enddo !print*,qsc(ifoc,ash) enddo print*,'------------- end volcanoes list -----------------' deallocate( date,cloud_column_height,so2,lon,lat,elevation) end subroutine get_emission_per_volcano !--------------------------------------------------------------- subroutine interp_volc_to_model(ng,n1,n2,n3,xt,yt,xm,ym,deltax,deltay,plat,plon,rlat& ,rlon,rland,nspecies_actual& ,actual_number_volc,nvolcanoes & ,qsc,qlon,qlat,qelev,qBEGT,qcrat,qheat,qplum_heigth & ,qdura,QESO2) use grid_dims_out, only: grid_type use mem_grid, only : grid_g implicit none integer, intent(in) ::actual_number_volc,nspecies_actual,nvolcanoes integer, intent(in) :: ng,n1,n2,n3 integer :: igbox,jgbox real deltax,deltay,plat,plon real, dimension(nvolcanoes) , intent(in):: & qlon,qlat,qelev,qBEGT& ,qcrat,qheat,qplum_heigth ,qdura ,QESO2 real, dimension(nvolcanoes,nspecies_actual), intent(in) :: qsc real xt(n1),yt(n2),xm(n1),ym(n2) real, dimension(n1,n2):: rlat,rlon,rland real, allocatable, dimension(:,:):: grid_area integer,dimension(n1,n2) :: num_volc_per_gbox integer,dimension(actual_number_volc) :: igbox_1d,jgbox_1d integer i,j,iesp,ifoc real wlon,elon,xlon,splon !tks num_volc_per_gbox=0 wlon=minval(rlon) !tks elon=maxval(rlon) !tks do ifoc = 1,actual_number_volc if(grid_type == 'rams' .or. grid_type == 'polar') then call get_ij_rams(n1,n2,xt,yt,xm,ym,plat,plon,deltax,deltay& ,qlon(ifoc),qlat(ifoc),igbox,jgbox) elseif(grid_type == 'lambert'.or. grid_type == 'mercator' ) then !tks xlon=qlon(ifoc) !tks if(wlon.lt.-180)then splon=wlon+360. !tks if(xlon.ge.splon)xlon=xlon-360. !tks endif if(elon.gt.180)then splon=elon-360. !tks if(xlon.le.splon)xlon=xlon+360. !tks endif call get_ij_rams(n1,n2,xt,yt,xm,ym,plat,plon,deltax,deltay& ,xlon,qlat(ifoc),igbox,jgbox) elseif(grid_type == 'll') then call get_ij_ll(n1,n2,rlat,rlon,qlon(ifoc),qlat(ifoc),igbox,jgbox) elseif(grid_type == 'fim') then call get_ij_fim(n1,n2,rlat,rlon,qlon(ifoc),qlat(ifoc),igbox,jgbox) elseif(grid_type == 'gg') then call get_ij_gg(n1,n2,rlat,rlon,qlon(ifoc),qlat(ifoc),igbox,jgbox) else stop 'unknown grid type at volcanoes degass emissions ' endif ! print*,'fire=',ifoc,igbox,jgbox,rlat(igbox,jgbox),rlon(igbox,jgbox) !print *,'ycit=',rlat(igbox,jgbox) !print *,'xcit=',rlon(igbox,jgbox)!,area(nprg),test !lat,lon !print*,'dummy=main(xcit,ycit)' igbox_1d(ifoc)=igbox jgbox_1d(ifoc)=jgbox if(igbox == -1 .or. jgbox == -1 ) cycle !print*,'volc=',ifoc,igbox,jgbox,rlat(igbox,jgbox),rlon(igbox,jgbox) !- number of volcanos per grid box num_volc_per_gbox(igbox,jgbox)=num_volc_per_gbox(igbox,jgbox)+1 !- here plum_heigth is the depth of the plume above the local terrain ! degassingP_g(plum_heigth )%prop(igbox,jgbox) = degassingP_g(plum_heigth )%prop(igbox,jgbox)+ & ! ! volc_prop(plum_heigth)%volcano(ifoc) - volc_prop(elev)%volcano(ifoc) do iesp=1,nspecies_actual !- total mass emitted per gridbox degassing_g(iesp)%src(igbox,jgbox,1) = degassing_g(iesp)%src(igbox,jgbox,1) & + qsc(ifoc,iesp) enddo ! iesp loop enddo ! ifoc loop ! !- get the mean cloud column weighted by SO2 mass for cases with more than one volcano per grid box do ifoc = 1,actual_number_volc igbox= igbox_1d(ifoc) jgbox= jgbox_1d(ifoc) if(igbox == -1 .or. jgbox == -1 ) cycle !! !!- here plum_heigth is the depth of the plume above the local terrain !! degassingP_g(plum_heigth )%prop(igbox,jgbox) = degassingP_g(plum_heigth )%prop(igbox,jgbox) +& !! (volc_prop(plum_heigth)%volcano(ifoc) - volc_prop(elev)%volcano(ifoc))*& !! qsc(ifoc,SO2)/degassing_g(SO2)%src(igbox,jgbox,1) !- we will keep the the depth of the plume above the sea level degassingP_g(plum_heigth )%prop(igbox,jgbox) = degassingP_g(plum_heigth )%prop(igbox,jgbox) +& volc_prop(plum_heigth)%volcano(ifoc)*& qsc(ifoc,SO2)/degassing_g(SO2)%src(igbox,jgbox,1) !- elevation degassingP_g(elev )%prop(igbox,jgbox) = degassingP_g(elev )%prop(igbox,jgbox) +& volc_prop(elev)%volcano(ifoc)*& qsc(ifoc,SO2)/degassing_g(SO2)%src(igbox,jgbox,1) !print*,ifoc,igbox,jgbox !print*,degassingP_g(plum_heigth )%prop(igbox,jgbox),(volc_prop(plum_heigth)%volcano(ifoc) - & !volc_prop(elev)%volcano(ifoc)),qsc(ifoc,SO2)/degassing_g(SO2)%src(igbox,jgbox,1) enddo ! ifoc loop ! !- get the mean cloud column for cases with more than one volcano per grid box !do i=1,n1 ! do j=1,n2 ! if(num_volc_per_gbox(i,j)>1) then ! degassingP_g(plum_heigth )%prop(i,j) = degassingP_g(plum_heigth )%prop(i,j) & ! /float(num_volc_per_gbox(i,j)) ! print*,'ng2=',i,j,num_volc_per_gbox(i,j),degassingP_g(plum_heigth )%prop(i,j) ! endif ! enddo !enddo print*,'Max plume height:',maxval(degassingP_g(plum_heigth )%prop(:,:)) !don't change of place these lines below ... !-------------- convert to surface flux (kg/m^2) !- calculate the grib box area allocate(grid_area(n1,n2)) if(grid_type == 'rams' .or. grid_type == 'polar') then call get_area_rams(grid_area,n1,n2,xt,yt,xm,ym)! elseif(grid_type == 'lambert' .or. grid_type == 'mercator' ) then grid_area(:,:)=1./(grid_g(ng)%dxt(:,:)*grid_g(ng)%dyt(:,:)) elseif(grid_type == 'll') then call get_area_ll(grid_area,n1,n2) elseif(grid_type == 'fim') then call get_area_fim(grid_area,n1,n2) elseif(grid_type == 'gg') then call get_area_gg(grid_area,n1,n2,rlat,rlon) else stop 'unknown grid type at volc degass emissions ' endif ! convert from mass [kg/day] to flux [kg/m^2/day] do iesp=1,nspecies_actual degassing_g(iesp)%src(:,:,1) = degassing_g(iesp)%src(:,:,1)/grid_area(:,:) print*,'volc_degass=',iesp,maxval(degassing_g(iesp)%src(:,:,1))!,maxval(grid_area) enddo ! deallocate(grid_area) end subroutine interp_volc_to_model !--------------------------------------------------------------- end module volc_degassing_emissions !--------------------------------------------------------------- !--------------------------------------------------------------- !--------------------------------------------------------------- subroutine process_degassing_volc(iyear,imon,iday,ihour,ng,ngrids,n1,n2,n3,rlat,rlon,rland& ,deltax,deltay,xt,yt,xm,ym,plat,plon) use grid_dims_out, only : degass_volc_data_dir use volc_degassing_emissions, only : volc_prop,nvolcanoes,nspecies,& LATI,LONG,ELEV,BEGT,CRAT,HEAT,plum_heigth,DURA,ESO2,mem_volc_prop,& get_emission_per_volcano,interp_volc_to_model, get_dimension_per_file_volcano implicit none integer, intent (in) :: iyear,imon,iday,ihour,ng,n1,n2,n3,ngrids real, intent (in), dimension(n1,n2) :: rlat,rlon,rland real, intent (in) :: deltax,deltay real, intent(in) :: xt(n1), yt(n2),xm(n1), ym(n2),plat,plon real, dimension(nvolcanoes,nspecies) :: qsc integer, save :: actual_number_volc integer :: nevents if(ng == 1) then ! just for the first grid is needed !- allocate memory for fires properties. call mem_volc_prop () call get_dimension_per_file_volcano(degass_volc_data_dir& ,iyear,imon,iday& , nevents ) !-- process the source emission estimate for each valid fire count and accumulate per !-- model grid box call get_emission_per_volcano(degass_volc_data_dir& , actual_number_volc,iyear,imon,iday& , qsc & , volc_prop(LONG)%volcano & , volc_prop(LATI)%volcano & , volc_prop(ELEV)%volcano & , volc_prop(BEGT)%volcano & , volc_prop(CRAT)%volcano & , volc_prop(HEAT)%volcano & , volc_prop(plum_heigth)%volcano & , volc_prop(DURA)%volcano & , volc_prop(ESO2)%volcano & , nevents ) endif !ng==1 !- performs the "gridding" call interp_volc_to_model(ng,n1,n2,n3,xt,yt,xm,ym,deltax,deltay,plat,plon& , rlat,rlon,rland & , nspecies,actual_number_volc,nvolcanoes & , qsc & , volc_prop(LONG)%volcano & , volc_prop(LATI)%volcano & , volc_prop(ELEV)%volcano & , volc_prop(BEGT)%volcano & , volc_prop(CRAT)%volcano & , volc_prop(HEAT)%volcano & , volc_prop(plum_heigth)%volcano & , volc_prop(DURA)%volcano & , volc_prop(ESO2)%volcano ) end subroutine process_degassing_volc !--------------------------------------------------------------- subroutine mem_degassing_volcanoes(n1,n2,n3) use volc_degassing_emissions implicit none integer i,j integer, intent(in) :: n1,n2,n3 if(.not. allocated(degassing_g )) allocate(degassing_g (nspecies)) if(.not. allocated(degassingP_g)) allocate(degassingP_g(nprop)) call nullify_volcanoes(degassing_g(:),degassingP_g(:)) call alloc_volcanoes (degassing_g(:),degassingP_g(:),n1,n2,n3) end subroutine mem_degassing_volcanoes !--------------------------------------------------------------- subroutine get_degass_volc_indentity(spc_name,ident) !use chem1_list use volc_degassing_emissions, only : volcanoes_nspecies=>nspecies& ,volcanoes_spc_name=>spc_name implicit none integer isp character (len=*), intent(in) :: spc_name integer , intent(out) :: ident do isp = 1,volcanoes_nspecies ident=-1 if(spc_name == volcanoes_spc_name(isp)) then print*,'==>volcanoes found for ',spc_name ident=isp return endif enddo !print*,'chem1-list specie ',trim(spc_name), ' does not match of any one of bbbem' !print*,'ident=',ident !stop 444 end subroutine get_degass_volc_indentity !--------------------------------------------------------------- subroutine check(status) use netcdf integer, intent ( in) :: status !if(status /= 0) then if(status /= nf90_noerr) then print *, trim(nf90_strerror(status)) stop 'netcdf error' !stop 2 end if end subroutine check !---------------------------------------------------------------