Changeset 3436 for trunk/LMDZ.GENERIC/libf

Timestamp:

Sep 24, 2024, 2:36:41 PM (4 weeks 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/libf/phystd

Files:

• physiq_mod.F90 (modified) (1 diff)
• volcano.F90 (modified) (17 diffs)

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

@@ -1509 +1509 @@
   ! -----------------------------------------------------
         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) + &

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

@@ -16 +16 @@
 !$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 +32 @@
 !$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 +57 @@
 !     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 +79 @@
 
 ! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW :
-! =============================================
+! ============================================
 !     Source flux (kg/s)
 !       REAL, parameter :: mmsource = 1E8  !Mastin et al. 2009
@@ -92 +103 @@
       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 +120 @@
       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 +144 @@
       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 +186 @@
   !=======================================================================
   
-  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 +195 @@
     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 +221 @@
     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 +260 @@
     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 +310 @@
         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 +340 @@
         !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 +381 @@
     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 +433 @@
       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 +474 @@
           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 +504 @@
       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)