Changeset 3436 for trunk


Ignore:
Timestamp:
Sep 24, 2024, 2:36:41 PM (2 months ago)
Author:
emillour
Message:

Generic PCM:
Updated version of "volcano.F90" : add setting duration of eruption
also added a vocano.def example in deftank
AB

Location:
trunk/LMDZ.GENERIC
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/changelog.txt

    r3435 r3436  
    19811981Bug fix in the Thermal Plume Model for generic tracers.
    19821982Error in computation of potential temperature.
     1983
     1984== 24/09/2024 == AB+EM
     1985Updated version of "volcano.F90" : add setting duration of eruption
     1986also added a vocano.def example in deftank

  • trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

    r3429 r3436  
    15091509  ! -----------------------------------------------------
    15101510        if (callvolcano) then
    1511           call volcano(ngrid,nlayer,pplev,pu,pv,pt,zzlev,zdqvolc,nq,massarea,&
     1511          call volcano(ngrid,nlayer,pplev,pu,pv,pt,zdqvolc,nq,massarea,&
    15121512                       zday,ptimestep)
    15131513          pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &

  • trunk/LMDZ.GENERIC/libf/phystd/volcano.F90

    r3299 r3436  
    1616!$OMP THREADPRIVATE(heightvals)
    1717
    18   integer,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,2)
     18  real,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,3)
    1919                                        ! for each volcano: interval
    2020                                        ! between each eruption in days,
    21                                         ! day offset from 0
     21                                        ! day offset from 0,
     22                                        ! and duration of eruption in days.
    2223!$OMP THREADPRIVATE(eruptfreqs)
    2324 
    … …  
    3132!$OMP THREADPRIVATE(ivolc)
    3233
     34  real,save,allocatable,protected :: dtmp(:) ! dimension (nvolc)
     35                                        ! time counter for eruption
     36!$OMP THREADPRIVATE(dtmp)
     37
     38  logical,save,allocatable,protected :: eruptbool(:) ! dimension (nvolc)
     39                                        ! flag to signal whether to erupt volcano or not.
     40!$OMP THREADPRIVATE(eruptbool)
     41
    3342CONTAINS
    3443
    35   SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,zzlev,ssource,nq, &
     44  SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,ssource,nq, &
    3645                     massarea,zday,ptimestep)
    3746
    38       use time_phylmdz_mod, only: daysec
     47      use time_phylmdz_mod, only: daysec ! day length (s)
     48      use vertical_layers_mod, only: pseudoalt !atm. layer pseudo-altitude (km)
    3949
    4050      IMPLICIT NONE
    … …  
    4757!     Adapted for the LMD/GCM by Laura Kerber, J.-B. Madeleine, Saira Hamid
    4858!     Adapted to be more generalisable by Ashwin Braude (02.11.2023)
     59!     Further optimised for parallelisation by Ehouarn Millour (12.04.2024)
    4960!
    5061!   Reference :
    … …  
    6879
    6980! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW :
    70 ! =============================================
     81! ============================================
    7182!     Source flux (kg/s)
    7283!       REAL, parameter :: mmsource = 1E8  !Mastin et al. 2009
    … …  
    92103      INTEGER,intent(in) ::  nlay            ! Number of vertical levels
    93104      REAL,intent(in) :: pplev(ngrid,nlay+1) ! Pressure between the layers (Pa)
    94       REAL,intent(in) :: zzlev(ngrid,nlay+1) ! height between the layers (km)
    95105      REAL,intent(in) :: wv(ngrid,nlay+1)    ! wind
    96106      REAL,intent(in) :: wu(ngrid,nlay+1)    ! wind
    … …  
    110120      if (firstcall) then
    111121        ! initialize everything the very first time this routine is called
    112         call read_volcano(nq,ngrid)
     122        call read_volcano(nq,ngrid,nlay)
     123        !Make sure eruptfreqs values are multiples of timestep
     124        ! within roundoffs of modulo hence comparison to 1.e-6*(ptimestep/daysec) and not 0
     125        do volci=1,nvolc
     126!          if((modulo(eruptfreqs(volci,1),(ptimestep/daysec)).ne.0).or.(modulo(eruptfreqs(volci,2),(ptimestep/daysec)).ne.0))then
     127           if ((modulo(eruptfreqs(volci,1),(ptimestep/daysec)).gt.1.e-6*(ptimestep/daysec)).or. &
     128               (modulo(eruptfreqs(volci,2),(ptimestep/daysec)).gt.1.e-6*(ptimestep/daysec))) then
     129            write(*,*) volci,eruptfreqs(volci,1)/(ptimestep/daysec),eruptfreqs(volci,2)/(ptimestep/daysec)
     130            call abort_physic('volcano','frequency and offset of eruption must both be multiples of timestep',1)
     131          endif
     132          if(eruptfreqs(volci,3)-eruptfreqs(volci,1).ge.0.0)then
     133            write(*,*) volci, eruptfreqs(volci,1), eruptfreqs(volci,3)
     134            call abort_physic('volcano','duration of eruption must be lower than frequency',1)
     135          endif
     136          write(*,*)'volcano: First eruption of volcano ',volci, ' scheduled in ', &
     137          modulo(zday-eruptfreqs(volci,2),eruptfreqs(volci,1)), 'days'
     138        enddo
    113139        firstcall=.false.
    114140      endif
    … …  
    118144      do volci=1,nvolc
    119145       
    120         !If time to erupt volcano
    121         if(modulo(nint(zday)-eruptfreqs(volci,2),eruptfreqs(volci,1)).eq.0)then
    122 
    123           !and then emit source flux at that point
    124           WRITE(*,*) 'Volcano ',volci,' Time: ', nint(zday), eruptfreqs(volci,:)
    125           WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',heightvals(volci)
     146        !If time to start erupting volcano
     147        if(modulo(nint((zday-eruptfreqs(volci,2))*daysec/ptimestep),nint(eruptfreqs(volci,1)*daysec/ptimestep)).eq.0)then
     148
     149          dtmp(volci) = eruptfreqs(volci,3)
     150          eruptbool(volci) = .true.
     151
     152          !emit source flux over the duration of the eruption
     153          if(dtmp(volci).lt.ptimestep/daysec)then
     154            WRITE(*,*) 'Volcano ',volci,' Time: ', zday, eruptfreqs(volci,:)
     155            l = heightvals(volci)
     156            ig=ivolc(volci)
     157            WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',pseudoalt(l), ' for ', dtmp(volci), ' days'
     158            do iq=1,nq
     159              ssource(ig,l,iq) = ssource(ig,l,iq)+&
     160                                          (dtmp(volci)*daysec/ptimestep)*volfluxes(volci,iq)/massarea(ig,l)
     161            enddo
     162            eruptbool(volci) = .false.
     163          endif
     164        endif
     165
     166        !If eruption lasts for more than a single timestep, need to continue erupting.
     167        if(eruptbool(volci))then
     168          WRITE(*,*) 'Volcano ',volci,' Time: ', zday, eruptfreqs(volci,:)
    126169          l = heightvals(volci)
    127170          ig=ivolc(volci)
     171          WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',pseudoalt(l), ' for ', dtmp(volci), ' more days'
     172
    128173          do iq=1,nq
    129174            ssource(ig,l,iq) = ssource(ig,l,iq)+&
    130                                           volfluxes(volci,iq)/massarea(ig,l)
     175                                        (MIN(dtmp(volci),ptimestep/daysec)*daysec/ptimestep)*volfluxes(volci,iq)/massarea(ig,l)
    131176          enddo
    132177
     178          dtmp(volci) = dtmp(volci) - ptimestep/daysec
     179          if(dtmp(volci).lt.0.0) eruptbool(volci) = .false.
    133180        endif
    134181
    … …  
    139186  !=======================================================================
    140187 
    141   SUBROUTINE read_volcano(nq,ngrid)
     188  SUBROUTINE read_volcano(nq,ngrid,nlay)
    142189  ! this routine reads in volcano.def and initializes module variables
    143190    USE mod_phys_lmdz_para, only: is_master, bcast
    … …  
    148195    use geometry_mod, only: boundslon, boundslat ! grid cell boundaries (rad)
    149196    use geometry_mod, only: longitude_deg, latitude_deg ! grid cell coordinates (deg)
     197    use vertical_layers_mod, only: pseudoalt !atm. layer pseudo-altitude (km)
    150198    IMPLICIT NONE
    151199
    152200    integer, intent(in) :: nq ! number of tracers
    153201    INTEGER, intent(in) :: ngrid ! number of atmospheric columns
     202    INTEGER, intent(in) :: nlay !
     203    real :: zlev(nlay+1) ! inter-layer altitudes (km)
    154204    real, allocatable :: latlonvals_total(:,:) !lat, lon position of each volcano
    155205    integer, allocatable :: heightvals_total(:) ! alt index of eruption
    156206    real, allocatable :: eruptfreqs_total(:,:) ! date and freq of eruption
    157207    real, allocatable :: volfluxes_total(:,:) ! injected tracer rate
     208    real, allocatable :: dtmp_total(:)
     209    logical, allocatable :: eruptbool_total(:)
    158210    integer :: ierr ! error return code
    159211    integer :: iq,l,volci,ivoltrac,nvoltrac,blank
    160     integer :: eruptfreq,eruptoffset
     212    real :: eruptfreq,eruptoffset,eruptdur
    161213    integer :: nvolctrac ! number of outgased tracers
    162214    integer :: ig,i
    163215    integer :: nvolc_total ! total (planet-wide) number of volcanoes
    164216    integer :: nvolc_check
    165     real    :: lat_volc,lon_volc,volflux
     217    real    :: lat_volc,lon_volc,volflux,alt_volc
    166218    character(len=500) :: voltrac(nq)
    167219    character(len=500) :: tracline   ! to read volcano.def line
    … …  
    169221    real :: boundslat_deg(nvertex) ! cell latitude boundaries (deg)
    170222    real :: minlat, minlon, maxlat, maxlon
     223    real tmp(nlay+1)
    171224    logical,allocatable :: belongs_here(:)
    172225    integer,allocatable :: volc_index(:)
    … …  
    207260    if (ierr/=0) then
    208261      write(*,*) trim(rname)//" failed allocation for latlonvals_total()"
    209       call abort_physic(rname,"failled allocation",1)
     262      call abort_physic(rname,"failed allocation",1)
    210263    endif
    211264    ALLOCATE(heightvals_total(nvolc_total),stat=ierr)
    212265    if (ierr/=0) then
    213266      write(*,*) trim(rname)//" failed allocation for heightvals_total()"
    214       call abort_physic(rname,"failled allocation",1)
    215     endif
    216     ALLOCATE(eruptfreqs_total(nvolc_total,2),stat=ierr)
     267      call abort_physic(rname,"failed allocation",1)
     268    endif
     269    ALLOCATE(eruptfreqs_total(nvolc_total,3),stat=ierr)
    217270    if (ierr/=0) then
    218271      write(*,*) trim(rname)//" failed allocation for eruptfreqs_total()"
    219       call abort_physic(rname,"failled allocation",1)
     272      call abort_physic(rname,"failed allocation",1)
    220273    endif
    221274    ALLOCATE(volfluxes_total(nvolc_total,nq),stat=ierr)
    222275    if (ierr/=0) then
    223276      write(*,*) trim(rname)//" failed allocation for volfluxes_total()"
    224       call abort_physic(rname,"failled allocation",1)
     277      call abort_physic(rname,"failed allocation",1)
     278    endif
     279    ALLOCATE(dtmp_total(nvolc_total),stat=ierr)
     280    if (ierr/=0) then
     281      write(*,*) trim(rname)//" failed allocation for dtmp_total()"
     282      call abort_physic(rname,"failed allocation",1)
     283    endif
     284    ALLOCATE(eruptbool_total(nvolc_total),stat=ierr)
     285    if (ierr/=0) then
     286      write(*,*) trim(rname)//" failed allocation for eruptbool_total()"
     287      call abort_physic(rname,"failed allocation",1)
    225288    endif
    226289    ! initialize volfluxes_total
    … …  
    247310        ENDDO
    248311        write(*,*)trim(rname)//': volcano number:',volci
    249         READ(407,*,iostat=ierr) lat_volc,lon_volc,l
     312        READ(407,*,iostat=ierr) lat_volc,lon_volc,alt_volc
     313        !Convert altitude in pressure coords to altitude in grid coords.
     314        ! First compute inter-layer pseudo altitudes (m)
     315        zlev(1)=0.
     316        do i=2,nlay
     317          ! inter-layer i is half way between mid-layers i and i-1
     318          ! and convert km to m
     319          zlev(i)=1.e3*(pseudoalt(i)+pseudoalt(i-1))/2.
     320        enddo
     321        !extrapolate topmost layer top boundary
     322        zlev(nlay+1)= 2*zlev(nlay)-zlev(nlay-1)
     323        ! Then identify eruption layer index
     324        do i=1,nlay+1
     325          tmp(i) = abs(zlev(i)-alt_volc*1000.0)
     326        enddo
     327        l = minloc(tmp,1)
     328        if(zlev(l).gt.alt_volc)then
     329          l = l-1
     330        endif
    250331        write(*,*)trim(rname)//' latitude:',lat_volc
    251332        write(*,*)trim(rname)//' longitude:',lon_volc
     333        write(*,*)trim(rname)//' altitude snapped to between',zlev(l), ' and ', zlev(l+1), 'm'
    252334        latlonvals_total(volci,1) = lat_volc
    253335        latlonvals_total(volci,2) = lon_volc
    … …  
    258340        !Read in eruption frequency (in multiples of iphysiq) and offset of eruption
    259341        READ(407,'(A)',iostat=ierr) tracline
    260         READ(tracline,*) eruptfreq, eruptoffset
     342        READ(tracline,*) eruptfreq, eruptoffset, eruptdur
    261343        write(*,*)trim(rname)//' eruption frequency (sols):',eruptfreq
    262344        write(*,*)trim(rname)//' eruption offset (sols):',eruptoffset
     345        write(*,*)trim(rname)//' eruption duration (sols):',eruptdur
    263346        eruptfreqs_total(volci,1) = eruptfreq
    264347        eruptfreqs_total(volci,2) = eruptoffset
     348        eruptfreqs_total(volci,3) = eruptdur
     349        dtmp_total(volci) = eruptdur
     350        eruptbool_total(volci) = .false.
    265351       
    266352        do ivoltrac=1,nvoltrac
    … …  
    295381    call bcast(eruptfreqs_total)
    296382    call bcast(volfluxes_total)
     383    call bcast(dtmp_total)
     384    call bcast(eruptbool_total)
    297385 
    298386    ! 3. Each core extracts and stores data relevant to its domain only
    … …  
    345433      if (ierr/=0) then
    346434        write(*,*) trim(rname)//" failed allocation for latlonvals()"
    347         call abort_physic(rname,"failled allocation",1)
     435        call abort_physic(rname,"failed allocation",1)
    348436      endif
    349437      allocate(heightvals(nvolc),stat=ierr)
    350438      if (ierr/=0) then
    351439        write(*,*) trim(rname)//" failed allocation for heightvals()"
    352         call abort_physic(rname,"failled allocation",1)
    353       endif
    354       allocate(eruptfreqs(nvolc,2),stat=ierr)
     440        call abort_physic(rname,"failed allocation",1)
     441      endif
     442      allocate(eruptfreqs(nvolc,3),stat=ierr)
    355443      if (ierr/=0) then
    356444        write(*,*) trim(rname)//" failed allocation for eruptfreqs()"
    357         call abort_physic(rname,"failled allocation",1)
     445        call abort_physic(rname,"failed allocation",1)
    358446      endif
    359447      allocate(volfluxes(nvolc,nq),stat=ierr)
    360448      if (ierr/=0) then
    361449        write(*,*) trim(rname)//" failed allocation for volfluxes()"
    362         call abort_physic(rname,"failled allocation",1)
     450        call abort_physic(rname,"failed allocation",1)
    363451      endif
    364452      allocate(ivolc(nvolc),stat=ierr)
    365453      if (ierr/=0) then
    366454        write(*,*) trim(rname)//" failed allocation for ivolc()"
    367         call abort_physic(rname,"failled allocation",1)
     455        call abort_physic(rname,"failed allocation",1)
     456      endif
     457      allocate(dtmp(nvolc),stat=ierr)
     458      if (ierr/=0) then
     459        write(*,*) trim(rname)//" failed allocation for dtmp()"
     460        call abort_physic(rname,"failed allocation",1)
     461      endif
     462      allocate(eruptbool(nvolc),stat=ierr)
     463      if (ierr/=0) then
     464        write(*,*) trim(rname)//" failed allocation for eruptbool()"
     465        call abort_physic(rname,"failed allocation",1)
    368466      endif
    369467    endif ! of if (nvolc>0)
    … …  
    376474          latlonvals(volci,1:2)=latlonvals_total(i,1:2)
    377475          heightvals(volci)=heightvals_total(i)
    378           eruptfreqs(volci,1:2)=eruptfreqs_total(i,1:2)
     476          eruptfreqs(volci,1:3)=eruptfreqs_total(i,1:3)
    379477          volfluxes(volci,1:nq)=volfluxes_total(i,1:nq)
    380478          ivolc(volci)=volc_index(i)
     479          dtmp(volci) = dtmp_total(i)
     480          eruptbool(volci) = eruptbool_total(i)
     481
    381482          ! the job is done, move on to next one
    382483          belongs_here(i)=.false.
    … …  
    403504      write(*,*) trim(rname)//": latlonvals(1:2)=",latlonvals(volci,1:2)
    404505      write(*,*) trim(rname)//": heightvals=",heightvals(volci)
    405       write(*,*) trim(rname)//": eruptfreqs(1:2)=",eruptfreqs(volci,1:2)
     506      write(*,*) trim(rname)//": eruptfreqs(1:3)=",eruptfreqs(volci,1:3)
    406507      write(*,*) trim(rname)//": volfluxes(1:nq)=",volfluxes(volci,1:nq)
    407508      write(*,*) trim(rname)//": ivolc=",ivolc(volci)
Note: See TracChangeset for help on using the changeset viewer.