MODULE volcano_mod IMPLICIT NONE ! saved,protected variables: (local to the core) integer,save,protected :: nvolc ! number of active volcanoes !$OMP THREADPRIVATE(nvolc) real,save,allocatable,protected :: latlonvals(:,:) ! dimension (nvolc,2) ! lat, lon position of each volcano (deg) !$OMP THREADPRIVATE(latlonvals) integer,save,allocatable,protected :: heightvals(:) ! dimension (nvolc) ! altitude grid point index ! of eruption for each volcano !$OMP THREADPRIVATE(heightvals) integer,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,2) ! for each volcano: interval ! between each eruption in days, ! day offset from 0 !$OMP THREADPRIVATE(eruptfreqs) real,save,allocatable,protected :: volfluxes(:,:) ! dimension (nvolc,nq) ! outgassing flux of each tracer for ! each volcano, in kg/s !$OMP THREADPRIVATE(volfluxes) integer,save,allocatable,protected :: ivolc(:) ! dimension (nvolc) ! horizontal grid index of volcano !$OMP THREADPRIVATE(ivolc) CONTAINS SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,zzlev,ssource,nq, & massarea,zday,ptimestep) use time_phylmdz_mod, only: daysec IMPLICIT NONE !======================================================================= ! Volcanic Activity (updated to run in parallel) ! ! Description : ! Author : Lionel Wilson ! Adapted for the LMD/GCM by Laura Kerber, J.-B. Madeleine, Saira Hamid ! Adapted to be more generalisable by Ashwin Braude (02.11.2023) ! ! Reference : ! @ARTICLE{2007JVGR..163...83W, ! author = {{Wilson}, L. and {Head}, J.~W.}, ! title = "{Explosive volcanic eruptions on Mars: Tephra and ! accretionary lapilli formation, dispersal and recognition in the ! geologic record}", ! journal = {Journal of Volcanology and Geothermal Research}, ! year = 2007, ! month = jun, ! volume = 163} !======================================================================= ! Variable statement ! __________________________________________________________________ ! Parameters ! ---------- ! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW : ! ============================================= ! Source flux (kg/s) ! REAL, parameter :: mmsource = 1E8 !Mastin et al. 2009 ! REAL, parameter :: wsource = 1.E6 !1wt% magma water content(Greeley 1987) ! REAL, parameter :: h2so4source = 1.E8 ! ============================================= ! Local variables ! --------------- INTEGER :: l INTEGER :: iq ! Tracer identifier INTEGER :: ig INTEGER :: volci LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) ! Inputs ! ------ INTEGER, intent(in) :: nq ! Number of tracers INTEGER, intent(in) :: ngrid ! Number of grid points INTEGER,intent(in) :: nlay ! Number of vertical levels REAL,intent(in) :: pplev(ngrid,nlay+1) ! Pressure between the layers (Pa) REAL,intent(in) :: zzlev(ngrid,nlay+1) ! height between the layers (km) REAL,intent(in) :: wv(ngrid,nlay+1) ! wind REAL,intent(in) :: wu(ngrid,nlay+1) ! wind REAL,intent(in) :: pt(ngrid,nlay+1) ! temp REAL,intent(in) :: massarea(ngrid,nlay) REAL,intent(in) :: zday ! "Universal time": in sols (and fraction thereof). ! since the North. Spring equinox (Ls=0) REAL,intent(in) :: ptimestep ! Physics timestep (s). ! Outputs ! ------- REAL, intent(out) :: ssource(ngrid,nlay,nq) ! Source tendency (kg.kg-1.s-1) ! Initialization ! __________________________________________________________________ if (firstcall) then ! initialize everything the very first time this routine is called call read_volcano(nq,ngrid) firstcall=.false. endif ssource(1:ngrid,1:nlay,1:nq)=0 !all arrays=zero since it's a local variable within routine that need to be initialized do volci=1,nvolc !If time to erupt volcano if(modulo(nint(zday)-eruptfreqs(volci,2),eruptfreqs(volci,1)).eq.0)then !and then emit source flux at that point WRITE(*,*) 'Volcano ',volci,' Time: ', nint(zday), eruptfreqs(volci,:) WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',heightvals(volci) l = heightvals(volci) ig=ivolc(volci) do iq=1,nq ssource(ig,l,iq) = ssource(ig,l,iq)+& volfluxes(volci,iq)/massarea(ig,l) enddo endif enddo ! of do volci=1,nvolc END SUBROUTINE volcano !======================================================================= SUBROUTINE read_volcano(nq,ngrid) ! this routine reads in volcano.def and initializes module variables USE mod_phys_lmdz_para, only: is_master, bcast use mod_phys_lmdz_transfert_para, only: reduce_sum use comcstfi_mod, only: pi use tracer_h, only: is_known_tracer, tracer_index use mod_grid_phy_lmdz, only: nvertex ! number of grid cell "corners" use geometry_mod, only: boundslon, boundslat ! grid cell boundaries (rad) use geometry_mod, only: longitude_deg, latitude_deg ! grid cell coordinates (deg) IMPLICIT NONE integer, intent(in) :: nq ! number of tracers INTEGER, intent(in) :: ngrid ! number of atmospheric columns real, allocatable :: latlonvals_total(:,:) !lat, lon position of each volcano integer, allocatable :: heightvals_total(:) ! alt index of eruption real, allocatable :: eruptfreqs_total(:,:) ! date and freq of eruption real, allocatable :: volfluxes_total(:,:) ! injected tracer rate integer :: ierr ! error return code integer :: iq,l,volci,ivoltrac,nvoltrac,blank integer :: eruptfreq,eruptoffset integer :: nvolctrac ! number of outgased tracers integer :: ig,i integer :: nvolc_total ! total (planet-wide) number of volcanoes integer :: nvolc_check real :: lat_volc,lon_volc,volflux character(len=500) :: voltrac(nq) character(len=500) :: tracline ! to read volcano.def line real :: boundslon_deg(nvertex) ! cell longitude boundaries (deg) real :: boundslat_deg(nvertex) ! cell latitude boundaries (deg) real :: minlat, minlon, maxlat, maxlon logical,allocatable :: belongs_here(:) integer,allocatable :: volc_index(:) character(len=*),parameter :: rname="read_volcano" ! 1. Only the master reads and processes the volcano.def file if (is_master) then ! 1.1 Open file and get total number of volcanoes open(407, form = 'formatted', status = 'old', & file = 'volcano.def', iostat=ierr) if (ierr /=0) then print*,trim(rname)//': Problem in opening volcano.def' stop end if write(*,*)trim(rname)//' Reading file volcano.def' DO READ(407,'(A)',iostat=ierr) tracline IF (ierr==0) THEN IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header !Read in the number of point sources/volcanoes in volcano.def READ(tracline,*) nvolc_total EXIT ENDIF ELSE ! If pb, or if reached EOF without having found nqtot ! call abort_physic('initracer','Unable to read numbers '// & ! 'of tracers in traceur.def',1) print*,trim(rname)//": unable to read numbers of tracer in volcano.def" stop ENDIF ENDDO endif ! of if (is_master) ! tell the world about nvolc_total CALL bcast(nvolc_total) ! 1.2. allocate arrays (all cores) to store all the information ALLOCATE(latlonvals_total(nvolc_total,2),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for latlonvals_total()" call abort_physic(rname,"failled allocation",1) endif ALLOCATE(heightvals_total(nvolc_total),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for heightvals_total()" call abort_physic(rname,"failled allocation",1) endif ALLOCATE(eruptfreqs_total(nvolc_total,2),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptfreqs_total()" call abort_physic(rname,"failled allocation",1) endif ALLOCATE(volfluxes_total(nvolc_total,nq),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for volfluxes_total()" call abort_physic(rname,"failled allocation",1) endif ! initialize volfluxes_total volfluxes_total(:,:)=0.0 ! 1.3. continue (master only) to read and process volcano.def if (is_master) then volci = 0 !Iterate for each volcano do while (volci.lt.nvolc_total) DO READ(407,'(A)',iostat=ierr) tracline IF (ierr==0) THEN IF (index(tracline,'#').ne.1) THEN READ(tracline,*) volci EXIT ENDIF ELSE ! If pb, or if reached EOF without having found nqtot ! call abort_physic('initracer','Unable to read numbers '// & ! 'of tracers in traceur.def',1) print*,trim(rname)//": unable to read line in volcano.def" stop ENDIF ENDDO write(*,*)trim(rname)//': volcano number:',volci READ(407,*,iostat=ierr) lat_volc,lon_volc,l write(*,*)trim(rname)//' latitude:',lat_volc write(*,*)trim(rname)//' longitude:',lon_volc latlonvals_total(volci,1) = lat_volc latlonvals_total(volci,2) = lon_volc heightvals_total(volci) = l !Read in number of active tracers being outgassed READ(407,*,iostat=ierr) nvoltrac write(*,*)trim(rname)//' number of outgased tracers:',nvoltrac !Read in eruption frequency (in multiples of iphysiq) and offset of eruption READ(407,'(A)',iostat=ierr) tracline READ(tracline,*) eruptfreq, eruptoffset write(*,*)trim(rname)//' eruption frequency (sols):',eruptfreq write(*,*)trim(rname)//' eruption offset (sols):',eruptoffset eruptfreqs_total(volci,1) = eruptfreq eruptfreqs_total(volci,2) = eruptoffset do ivoltrac=1,nvoltrac !Read in tracer name and check it's in traceur.def read(407,'(A)') tracline blank = index(tracline,' ') ! Find position of 1st blank in tracline voltrac(ivoltrac) = tracline(1:blank) ! check it is a valid tracer: if (is_known_tracer(voltrac(ivoltrac))) then !Read in flux of tracer out of volcano (kg/s) READ(407,'(A)',iostat=ierr) tracline READ(tracline,*) volflux iq=tracer_index(voltrac(ivoltrac)) volfluxes_total(volci,iq) = volflux else write(*,*) trim(rname)//": problem with ",trim(voltrac(ivoltrac)) call abort_physic(rname, 'tracer '//& trim(voltrac(ivoltrac))//' not found', 1) endif enddo ! of do ivoltrac=1,nvoltrac enddo ! of do while (volci.lt.nvolc_total) close(407) endif ! of if (is_master) ! 2. Share data with all cores call bcast(latlonvals_total) call bcast(heightvals_total) call bcast(eruptfreqs_total) call bcast(volfluxes_total) ! 3. Each core extracts and stores data relevant to its domain only ! 3.1. find out how many volcanoes for this core nvolc=0 allocate(belongs_here(nvolc_total),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for belong_here()" call abort_physic(rname,"failled allocation",1) endif allocate(volc_index(nvolc_total),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for volc_index()" call abort_physic(rname,"failled allocation",1) endif do volci=1,nvolc_total belongs_here(volci)=.false. volc_index(volci)=0 lat_volc=latlonvals_total(volci,1) lon_volc=latlonvals_total(volci,2) ! find out if this volcano is in the ig cell do ig=1,ngrid boundslon_deg(:)=boundslon(ig,:)/pi*180. boundslat_deg(:)=boundslat(ig,:)/pi*180. minlon=minval(boundslon_deg) maxlon=maxval(boundslon_deg) minlat=minval(boundslat_deg) maxlat=maxval(boundslat_deg) if ((minlat<=lat_volc).and.(lat_volc<=maxlat).and. & (((minlon<=lon_volc).and.(lon_volc<=maxlon)).or. & ((minlon+360.<=lon_volc).and.(lon_volc<=maxlon+360.)).or. & ((minlon-360.<=lon_volc).and.(lon_volc<=maxlon-360.)))) then nvolc=nvolc+1 belongs_here(volci)=.true. volc_index(volci)=ig ! the job is done move on to next volcano exit endif enddo ! of do ig=1,ngrid enddo ! of volci=1,nvolc_total ! 3.2. allocate arrays, now that nvolc is known if (nvolc>0) then allocate(latlonvals(nvolc,2),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for latlonvals()" call abort_physic(rname,"failled allocation",1) endif allocate(heightvals(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for heightvals()" call abort_physic(rname,"failled allocation",1) endif allocate(eruptfreqs(nvolc,2),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptfreqs()" call abort_physic(rname,"failled allocation",1) endif allocate(volfluxes(nvolc,nq),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for volfluxes()" call abort_physic(rname,"failled allocation",1) endif allocate(ivolc(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for ivolc()" call abort_physic(rname,"failled allocation",1) endif endif ! of if (nvolc>0) ! 3.3. copy over global data to local arrays do volci=1,nvolc do i=1,nvolc_total if (belongs_here(i)) then ! copy data to local arrays latlonvals(volci,1:2)=latlonvals_total(i,1:2) heightvals(volci)=heightvals_total(i) eruptfreqs(volci,1:2)=eruptfreqs_total(i,1:2) volfluxes(volci,1:nq)=volfluxes_total(i,1:nq) ivolc(volci)=volc_index(i) ! the job is done, move on to next one belongs_here(i)=.false. exit endif ! of if (belongs_here(i)) enddo ! of do i=1,nvolc_total enddo ! of do volci=1,nvolc ! 4. Sanity check ! verify that the sum of volcanoes over all cores is indeed nvolc_total if (is_master) nvolc_check=0 call bcast(nvolc_check) call reduce_sum(nvolc,nvolc_check) call bcast(nvolc_check) if (nvolc_check/=nvolc_total) then write(*,*) trim(rname)//" error. Sum of volcanoes over all cores ", & nvolc_check," differs from nvolc_total=",nvolc_total call abort_physic(rname," error. Check and fix volcano.F90",1) endif ! give a summary, for each core: write(*,*) trim(rname)//": nvolc=",nvolc," /",nvolc_total do volci=1,nvolc write(*,*) trim(rname)//": latlonvals(1:2)=",latlonvals(volci,1:2) write(*,*) trim(rname)//": heightvals=",heightvals(volci) write(*,*) trim(rname)//": eruptfreqs(1:2)=",eruptfreqs(volci,1:2) write(*,*) trim(rname)//": volfluxes(1:nq)=",volfluxes(volci,1:nq) write(*,*) trim(rname)//": ivolc=",ivolc(volci) write(*,*) trim(rname)//": longitude(ivolc)=",longitude_deg(ivolc(volci)) write(*,*) trim(rname)//": latitude(ivolc)=",latitude_deg(ivolc(volci)) enddo END SUBROUTINE read_volcano END MODULE volcano_mod