Index: /trunk/LMDZ.GENERIC/changelog.txt =================================================================== --- /trunk/LMDZ.GENERIC/changelog.txt (revision 3435) +++ /trunk/LMDZ.GENERIC/changelog.txt (revision 3436) @@ -1981,2 +1981,6 @@ Bug fix in the Thermal Plume Model for generic tracers. Error in computation of potential temperature. + +== 24/09/2024 == AB+EM +Updated version of "volcano.F90" : add setting duration of eruption +also added a vocano.def example in deftank Index: /trunk/LMDZ.GENERIC/deftank/volcano.def =================================================================== --- /trunk/LMDZ.GENERIC/deftank/volcano.def (revision 3436) +++ /trunk/LMDZ.GENERIC/deftank/volcano.def (revision 3436) @@ -0,0 +1,73 @@ +#volcano.def - define parameters of volcanoes on the surface of Mars +#Structure file as follows (comments are preceded by a hash sign): +# +#Number of volcanoes +#Volcano 1 +#Number of volcano in sequence (in this case 1) +#lat lon altitude +#number of tracers outgassed +#eruption frequency, offset and length (in days) +#name of tracer 1 +#flux of tracer 1 (kg/s) +#name of tracer 2 +#flux of tracer 2 (kg/s) etc. +#Volcano 2 etc +2 +#Volcano 1: Hadriacus Mons +1 +-30.44 92.18 2.0 +8 +61 3.25 3.0 +h2o_vap +2.39e8 +h2 +3.95e6 +h2s_vap +2.81e7 +so2 +7.35e7 +s2_vap +5.66e7 +co +1.73e7 +hcl_vap +4.41e7 +co2 +5.83e7 +#Volcano 2: Hadriacus Mons +2 +-8.5 174.4 20.0 +7 +31 10.25 1.0 +h2o_vap +2.39e8 +h2 +3.95e6 +h2s_vap +2.81e7 +so2 +7.35e7 +s2_vap +5.66e7 +co +1.73e7 +co2 +5.83e7 + +# ex : Albor Tholus lon=150.4 and lat= 19.0 (Robbins2011) +# ex : Amphritites lon=52.66 and lat= -58.00 +# ex : Apollinaris Patera lon=174.4 and lat=-8.5 +# ex : Arsia Mons lon=-120.46 and lat= -9.14 +# ex : Ascraeus Mons lon=-104.37 and lat= 11.1 +# ex : Cerberus Fossae lon= 160.0 and lat= 10.0 +# ex : Electris volcano lon=-173.21 and lat = -37.35 +# ex : Elysium Mons lon=147 and lat= 24.8 +# ex : Hadriacus Mons lon=92.18 and lat= -30.44 +# ex : Hecates Tholus lon=150.08 and lat= 31.68 +# ex : Malea Patera lon=50.96 and lat= -63.09 +# ex : Olympus Mons lon=-133.9 and lat =18.7 +# ex : Pavonis Mons lon=-112.85 and lat= 0.662 +# ex : Peneus Patera lon=60.76 and lat= -58.05 +# ex : Pityusa Patera lon= 36.87 and lat= -66.77 +# ex : Syrtis Major lon=66.4 and lat= 9.85 +# ex : Tyrrhenia Mons lon=106.55 and lat= -21.32 Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3435) +++ /trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90 (revision 3436) @@ -1509,5 +1509,5 @@ ! ----------------------------------------------------- if (callvolcano) then - call volcano(ngrid,nlayer,pplev,pu,pv,pt,zzlev,zdqvolc,nq,massarea,& + call volcano(ngrid,nlayer,pplev,pu,pv,pt,zdqvolc,nq,massarea,& zday,ptimestep) pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + & Index: /trunk/LMDZ.GENERIC/libf/phystd/volcano.F90 =================================================================== --- /trunk/LMDZ.GENERIC/libf/phystd/volcano.F90 (revision 3435) +++ /trunk/LMDZ.GENERIC/libf/phystd/volcano.F90 (revision 3436) @@ -16,8 +16,9 @@ !$OMP THREADPRIVATE(heightvals) - integer,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,2) + real,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,3) ! for each volcano: interval ! between each eruption in days, - ! day offset from 0 + ! day offset from 0, + ! and duration of eruption in days. !$OMP THREADPRIVATE(eruptfreqs) @@ -31,10 +32,19 @@ !$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,zzlev,ssource,nq, & + SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,ssource,nq, & massarea,zday,ptimestep) - use time_phylmdz_mod, only: daysec + use time_phylmdz_mod, only: daysec ! day length (s) + use vertical_layers_mod, only: pseudoalt !atm. layer pseudo-altitude (km) IMPLICIT NONE @@ -47,4 +57,5 @@ ! 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 : @@ -68,5 +79,5 @@ ! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW : -! ============================================= +! ============================================ ! Source flux (kg/s) ! REAL, parameter :: mmsource = 1E8 !Mastin et al. 2009 @@ -92,5 +103,4 @@ 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 @@ -110,5 +120,21 @@ if (firstcall) then ! initialize everything the very first time this routine is called - call read_volcano(nq,ngrid) + 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 @@ -118,17 +144,38 @@ 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) + !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)+& - volfluxes(volci,iq)/massarea(ig,l) + (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 @@ -139,5 +186,5 @@ !======================================================================= - SUBROUTINE read_volcano(nq,ngrid) + 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 @@ -148,20 +195,25 @@ 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 - integer :: eruptfreq,eruptoffset + 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 + real :: lat_volc,lon_volc,volflux,alt_volc character(len=500) :: voltrac(nq) character(len=500) :: tracline ! to read volcano.def line @@ -169,4 +221,5 @@ 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(:) @@ -207,20 +260,30 @@ if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for latlonvals_total()" - call abort_physic(rname,"failled allocation",1) + 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,"failled allocation",1) - endif - ALLOCATE(eruptfreqs_total(nvolc_total,2),stat=ierr) + 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,"failled allocation",1) + 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,"failled allocation",1) + 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 @@ -247,7 +310,26 @@ ENDDO write(*,*)trim(rname)//': volcano number:',volci - READ(407,*,iostat=ierr) lat_volc,lon_volc,l + 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 @@ -258,9 +340,13 @@ !Read in eruption frequency (in multiples of iphysiq) and offset of eruption READ(407,'(A)',iostat=ierr) tracline - READ(tracline,*) eruptfreq, eruptoffset + 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 @@ -295,4 +381,6 @@ 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 @@ -345,25 +433,35 @@ if (ierr/=0) then write(*,*) trim(rname)//" failed allocation for latlonvals()" - call abort_physic(rname,"failled allocation",1) + 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,"failled allocation",1) - endif - allocate(eruptfreqs(nvolc,2),stat=ierr) + 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,"failled allocation",1) + 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,"failled allocation",1) + 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,"failled allocation",1) + 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) @@ -376,7 +474,10 @@ latlonvals(volci,1:2)=latlonvals_total(i,1:2) heightvals(volci)=heightvals_total(i) - eruptfreqs(volci,1:2)=eruptfreqs_total(i,1:2) + 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. @@ -403,5 +504,5 @@ 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)//": eruptfreqs(1:3)=",eruptfreqs(volci,1:3) write(*,*) trim(rname)//": volfluxes(1:nq)=",volfluxes(volci,1:nq) write(*,*) trim(rname)//": ivolc=",ivolc(volci)