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) real,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,3) ! for each volcano: interval ! between each eruption in days, ! day offset from 0, ! and duration of eruption in days. !$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) real,save,allocatable,protected :: dtmp(:) ! dimension (nvolc) ! time counter for eruption !$OMP THREADPRIVATE(dtmp) logical,save,allocatable,protected :: eruptbool(:) ! dimension (nvolc) ! flag to signal whether to erupt volcano or not. !$OMP THREADPRIVATE(eruptbool) CONTAINS SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,ssource,nq, & massarea,zday,ptimestep) use time_phylmdz_mod, only: daysec ! day length (s) use vertical_layers_mod, only: pseudoalt !atm. layer pseudo-altitude (km) 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) ! Further optimised for parallelisation by Ehouarn Millour (12.04.2024) ! ! 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) :: 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,nlay) !Make sure eruptfreqs values are multiples of timestep ! within roundoffs of modulo hence comparison to 1.e-6*(ptimestep/daysec) and not 0 do volci=1,nvolc ! if((modulo(eruptfreqs(volci,1),(ptimestep/daysec)).ne.0).or.(modulo(eruptfreqs(volci,2),(ptimestep/daysec)).ne.0))then if ((modulo(eruptfreqs(volci,1),(ptimestep/daysec)).gt.1.e-6*(ptimestep/daysec)).or. & (modulo(eruptfreqs(volci,2),(ptimestep/daysec)).gt.1.e-6*(ptimestep/daysec))) then write(*,*) volci,eruptfreqs(volci,1)/(ptimestep/daysec),eruptfreqs(volci,2)/(ptimestep/daysec) call abort_physic('volcano','frequency and offset of eruption must both be multiples of timestep',1) endif if(eruptfreqs(volci,3)-eruptfreqs(volci,1).ge.0.0)then write(*,*) volci, eruptfreqs(volci,1), eruptfreqs(volci,3) call abort_physic('volcano','duration of eruption must be lower than frequency',1) endif write(*,*)'volcano: First eruption of volcano ',volci, ' scheduled in ', & modulo(zday-eruptfreqs(volci,2),eruptfreqs(volci,1)), 'days' enddo 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 start erupting volcano if(modulo(nint((zday-eruptfreqs(volci,2))*daysec/ptimestep),nint(eruptfreqs(volci,1)*daysec/ptimestep)).eq.0)then dtmp(volci) = eruptfreqs(volci,3) eruptbool(volci) = .true. !emit source flux over the duration of the eruption if(dtmp(volci).lt.ptimestep/daysec)then WRITE(*,*) 'Volcano ',volci,' Time: ', zday, eruptfreqs(volci,:) l = heightvals(volci) ig=ivolc(volci) WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',pseudoalt(l), ' for ', dtmp(volci), ' days' do iq=1,nq ssource(ig,l,iq) = ssource(ig,l,iq)+& (dtmp(volci)*daysec/ptimestep)*volfluxes(volci,iq)/massarea(ig,l) enddo eruptbool(volci) = .false. endif endif !If eruption lasts for more than a single timestep, need to continue erupting. if(eruptbool(volci))then WRITE(*,*) 'Volcano ',volci,' Time: ', zday, eruptfreqs(volci,:) l = heightvals(volci) ig=ivolc(volci) WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',pseudoalt(l), ' for ', dtmp(volci), ' more days' do iq=1,nq ssource(ig,l,iq) = ssource(ig,l,iq)+& (MIN(dtmp(volci),ptimestep/daysec)*daysec/ptimestep)*volfluxes(volci,iq)/massarea(ig,l) enddo dtmp(volci) = dtmp(volci) - ptimestep/daysec if(dtmp(volci).lt.0.0) eruptbool(volci) = .false. endif enddo ! of do volci=1,nvolc END SUBROUTINE volcano !======================================================================= SUBROUTINE read_volcano(nq,ngrid,nlay) ! 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) use vertical_layers_mod, only: pseudoalt !atm. layer pseudo-altitude (km) IMPLICIT NONE integer, intent(in) :: nq ! number of tracers INTEGER, intent(in) :: ngrid ! number of atmospheric columns INTEGER, intent(in) :: nlay ! real :: zlev(nlay+1) ! inter-layer altitudes (km) 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 real, allocatable :: dtmp_total(:) logical, allocatable :: eruptbool_total(:) integer :: ierr ! error return code integer :: iq,l,volci,ivoltrac,nvoltrac,blank real :: eruptfreq,eruptoffset,eruptdur 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,alt_volc 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 real tmp(nlay+1) 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,"failed 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,"failed allocation",1) endif ALLOCATE(eruptfreqs_total(nvolc_total,3),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptfreqs_total()" call abort_physic(rname,"failed 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,"failed allocation",1) endif ALLOCATE(dtmp_total(nvolc_total),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for dtmp_total()" call abort_physic(rname,"failed allocation",1) endif ALLOCATE(eruptbool_total(nvolc_total),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptbool_total()" call abort_physic(rname,"failed 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,alt_volc !Convert altitude in pressure coords to altitude in grid coords. ! First compute inter-layer pseudo altitudes (m) zlev(1)=0. do i=2,nlay ! inter-layer i is half way between mid-layers i and i-1 ! and convert km to m zlev(i)=1.e3*(pseudoalt(i)+pseudoalt(i-1))/2. enddo !extrapolate topmost layer top boundary zlev(nlay+1)= 2*zlev(nlay)-zlev(nlay-1) ! Then identify eruption layer index do i=1,nlay+1 tmp(i) = abs(zlev(i)-alt_volc*1000.0) enddo l = minloc(tmp,1) if(zlev(l).gt.alt_volc)then l = l-1 endif write(*,*)trim(rname)//' latitude:',lat_volc write(*,*)trim(rname)//' longitude:',lon_volc write(*,*)trim(rname)//' altitude snapped to between',zlev(l), ' and ', zlev(l+1), 'm' 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, eruptdur write(*,*)trim(rname)//' eruption frequency (sols):',eruptfreq write(*,*)trim(rname)//' eruption offset (sols):',eruptoffset write(*,*)trim(rname)//' eruption duration (sols):',eruptdur eruptfreqs_total(volci,1) = eruptfreq eruptfreqs_total(volci,2) = eruptoffset eruptfreqs_total(volci,3) = eruptdur dtmp_total(volci) = eruptdur eruptbool_total(volci) = .false. 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) call bcast(dtmp_total) call bcast(eruptbool_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,"failed allocation",1) endif allocate(heightvals(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for heightvals()" call abort_physic(rname,"failed allocation",1) endif allocate(eruptfreqs(nvolc,3),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptfreqs()" call abort_physic(rname,"failed allocation",1) endif allocate(volfluxes(nvolc,nq),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for volfluxes()" call abort_physic(rname,"failed allocation",1) endif allocate(ivolc(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for ivolc()" call abort_physic(rname,"failed allocation",1) endif allocate(dtmp(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for dtmp()" call abort_physic(rname,"failed allocation",1) endif allocate(eruptbool(nvolc),stat=ierr) if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for eruptbool()" call abort_physic(rname,"failed 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:3)=eruptfreqs_total(i,1:3) volfluxes(volci,1:nq)=volfluxes_total(i,1:nq) ivolc(volci)=volc_index(i) dtmp(volci) = dtmp_total(i) eruptbool(volci) = eruptbool_total(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:3)=",eruptfreqs(volci,1:3) 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