source: trunk/LMDZ.GENERIC/libf/phystd/volcano.F90 @ 3300

Last change on this file since 3300 was 3299, checked in by emillour, 8 months ago

Generic PCM:

  • add some auxiliary functions in module tracer_h.F90: is_known_tracer(tname) returns ".true." if "tname" is a known tracer tracer_index(tname) returns the tracer index of tracer "tname"
  • add the possibility of "volcanic outgasing" as point sources. controled by a "callvolcano" flag (default .false.) details about source location and injection details (frequency, length) need be provided in a "volcano.def" ASCII file read at runtime. (see routine "read_volcano" in volcano.F90)

AB+EM

File size: 16.0 KB
RevLine 
[3299]1MODULE volcano_mod
2
3IMPLICIT NONE
4  ! saved,protected variables: (local to the core)
5 
6  integer,save,protected :: nvolc ! number of active volcanoes
7!$OMP THREADPRIVATE(nvolc)
8
9  real,save,allocatable,protected :: latlonvals(:,:) ! dimension (nvolc,2)
10                                     ! lat, lon position of each volcano (deg)
11!$OMP THREADPRIVATE(latlonvals)
12 
13  integer,save,allocatable,protected :: heightvals(:)   ! dimension (nvolc)
14                                        ! altitude grid point index
15                                        ! of eruption for each volcano
16!$OMP THREADPRIVATE(heightvals)
17
18  integer,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,2)
19                                        ! for each volcano: interval
20                                        ! between each eruption in days,
21                                        ! day offset from 0
22!$OMP THREADPRIVATE(eruptfreqs)
23 
24  real,save,allocatable,protected :: volfluxes(:,:) ! dimension (nvolc,nq)
25                                     ! outgassing flux of each tracer for
26                                     ! each volcano, in kg/s
27!$OMP THREADPRIVATE(volfluxes)
28
29  integer,save,allocatable,protected :: ivolc(:) ! dimension (nvolc)
30                                        ! horizontal grid index of volcano
31!$OMP THREADPRIVATE(ivolc)
32
33CONTAINS
34
35  SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,zzlev,ssource,nq, &
36                     massarea,zday,ptimestep)
37
38      use time_phylmdz_mod, only: daysec
39
40      IMPLICIT NONE
41
42!=======================================================================
43!   Volcanic Activity (updated to run in parallel)
44!
45!   Description :
46!     Author : Lionel Wilson
47!     Adapted for the LMD/GCM by Laura Kerber, J.-B. Madeleine, Saira Hamid
48!     Adapted to be more generalisable by Ashwin Braude (02.11.2023)
49!
50!   Reference :
51!       @ARTICLE{2007JVGR..163...83W,
52!       author = {{Wilson}, L. and {Head}, J.~W.},
53!       title = "{Explosive volcanic eruptions on Mars: Tephra and
54!       accretionary lapilli formation, dispersal and recognition in the
55!       geologic record}",
56!       journal = {Journal of Volcanology and Geothermal Research},
57!       year = 2007,
58!       month = jun,
59!       volume = 163}
60!=======================================================================
61
62!     Variable statement
63!     __________________________________________________________________
64
65
66! Parameters
67! ----------
68
69! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW :
70! =============================================
71!     Source flux (kg/s)
72!       REAL, parameter :: mmsource = 1E8  !Mastin et al. 2009
73!       REAL, parameter :: wsource = 1.E6 !1wt% magma water content(Greeley 1987)
74!       REAL, parameter :: h2so4source = 1.E8
75! =============================================
76! Local variables
77! ---------------
78
79      INTEGER :: l
80      INTEGER :: iq                     ! Tracer identifier
81      INTEGER :: ig
82      INTEGER :: volci
83
84      LOGICAL,SAVE :: firstcall=.true.
85!$OMP THREADPRIVATE(firstcall)
86
87! Inputs
88! ------
89
90      INTEGER, intent(in) :: nq               ! Number of tracers
91      INTEGER, intent(in) :: ngrid            ! Number of grid points
92      INTEGER,intent(in) ::  nlay            ! Number of vertical levels
93      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)
95      REAL,intent(in) :: wv(ngrid,nlay+1)    ! wind
96      REAL,intent(in) :: wu(ngrid,nlay+1)    ! wind
97      REAL,intent(in) :: pt(ngrid,nlay+1)    ! temp
98      REAL,intent(in) :: massarea(ngrid,nlay)
99      REAL,intent(in) :: zday ! "Universal time": in sols (and fraction thereof).
100                              ! since the North. Spring equinox (Ls=0)
101      REAL,intent(in) :: ptimestep             ! Physics timestep (s).
102
103! Outputs
104! -------
105      REAL, intent(out) :: ssource(ngrid,nlay,nq)    ! Source tendency (kg.kg-1.s-1)
106
107!     Initialization
108!     __________________________________________________________________
109   
110      if (firstcall) then
111        ! initialize everything the very first time this routine is called
112        call read_volcano(nq,ngrid)
113        firstcall=.false.
114      endif
115     
116      ssource(1:ngrid,1:nlay,1:nq)=0 !all arrays=zero since it's a local variable within routine that need to be initialized
117     
118      do volci=1,nvolc
119       
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)
126          l = heightvals(volci)
127          ig=ivolc(volci)
128          do iq=1,nq
129            ssource(ig,l,iq) = ssource(ig,l,iq)+&
130                                          volfluxes(volci,iq)/massarea(ig,l)
131          enddo
132
133        endif
134
135      enddo ! of do volci=1,nvolc
136
137  END SUBROUTINE volcano
138 
139  !=======================================================================
140 
141  SUBROUTINE read_volcano(nq,ngrid)
142  ! this routine reads in volcano.def and initializes module variables
143    USE mod_phys_lmdz_para, only: is_master, bcast
144    use mod_phys_lmdz_transfert_para, only: reduce_sum
145    use comcstfi_mod, only: pi
146    use tracer_h, only: is_known_tracer, tracer_index
147    use mod_grid_phy_lmdz,  only: nvertex ! number of grid cell "corners"
148    use geometry_mod, only: boundslon, boundslat ! grid cell boundaries (rad)
149    use geometry_mod, only: longitude_deg, latitude_deg ! grid cell coordinates (deg)
150    IMPLICIT NONE
151
152    integer, intent(in) :: nq ! number of tracers
153    INTEGER, intent(in) :: ngrid ! number of atmospheric columns
154    real, allocatable :: latlonvals_total(:,:) !lat, lon position of each volcano
155    integer, allocatable :: heightvals_total(:) ! alt index of eruption
156    real, allocatable :: eruptfreqs_total(:,:) ! date and freq of eruption
157    real, allocatable :: volfluxes_total(:,:) ! injected tracer rate
158    integer :: ierr ! error return code
159    integer :: iq,l,volci,ivoltrac,nvoltrac,blank
160    integer :: eruptfreq,eruptoffset
161    integer :: nvolctrac ! number of outgased tracers
162    integer :: ig,i
163    integer :: nvolc_total ! total (planet-wide) number of volcanoes
164    integer :: nvolc_check
165    real    :: lat_volc,lon_volc,volflux
166    character(len=500) :: voltrac(nq)
167    character(len=500) :: tracline   ! to read volcano.def line
168    real :: boundslon_deg(nvertex) ! cell longitude boundaries (deg)
169    real :: boundslat_deg(nvertex) ! cell latitude boundaries (deg)
170    real :: minlat, minlon, maxlat, maxlon
171    logical,allocatable :: belongs_here(:)
172    integer,allocatable :: volc_index(:)
173    character(len=*),parameter :: rname="read_volcano"
174
175    ! 1. Only the master reads and processes the volcano.def file
176    if (is_master) then
177      ! 1.1 Open file and get total number of volcanoes
178        open(407, form = 'formatted', status = 'old', &
179                         file = 'volcano.def', iostat=ierr)
180        if (ierr /=0) then
181          print*,trim(rname)//': Problem in opening volcano.def'
182          stop
183        end if
184       
185        write(*,*)trim(rname)//' Reading file volcano.def'
186        DO
187          READ(407,'(A)',iostat=ierr) tracline
188          IF (ierr==0) THEN
189            IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
190              !Read in the number of point sources/volcanoes in volcano.def
191              READ(tracline,*) nvolc_total
192              EXIT
193            ENDIF
194          ELSE ! If pb, or if reached EOF without having found nqtot
195          !        call abort_physic('initracer','Unable to read numbers '// &
196          !  'of tracers in traceur.def',1)
197              print*,trim(rname)//": unable to read numbers of tracer in volcano.def"
198              stop
199          ENDIF
200        ENDDO
201    endif ! of if (is_master)
202    ! tell the world about nvolc_total
203    CALL bcast(nvolc_total)
204   
205    ! 1.2. allocate arrays (all cores) to store all the information
206    ALLOCATE(latlonvals_total(nvolc_total,2),stat=ierr)
207    if (ierr/=0) then
208      write(*,*) trim(rname)//" failed allocation for latlonvals_total()"
209      call abort_physic(rname,"failled allocation",1)
210    endif
211    ALLOCATE(heightvals_total(nvolc_total),stat=ierr)
212    if (ierr/=0) then
213      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)
217    if (ierr/=0) then
218      write(*,*) trim(rname)//" failed allocation for eruptfreqs_total()"
219      call abort_physic(rname,"failled allocation",1)
220    endif
221    ALLOCATE(volfluxes_total(nvolc_total,nq),stat=ierr)
222    if (ierr/=0) then
223      write(*,*) trim(rname)//" failed allocation for volfluxes_total()"
224      call abort_physic(rname,"failled allocation",1)
225    endif
226    ! initialize volfluxes_total
227    volfluxes_total(:,:)=0.0
228   
229    ! 1.3. continue (master only) to read and process volcano.def
230    if (is_master) then
231      volci = 0
232      !Iterate for each volcano
233      do while (volci.lt.nvolc_total)
234        DO
235          READ(407,'(A)',iostat=ierr) tracline
236          IF (ierr==0) THEN
237            IF (index(tracline,'#').ne.1) THEN
238              READ(tracline,*) volci
239              EXIT
240            ENDIF
241          ELSE ! If pb, or if reached EOF without having found nqtot
242        !        call abort_physic('initracer','Unable to read numbers '// &
243        !  'of tracers in traceur.def',1)
244            print*,trim(rname)//": unable to read line in volcano.def"
245            stop
246          ENDIF
247        ENDDO
248        write(*,*)trim(rname)//': volcano number:',volci
249        READ(407,*,iostat=ierr) lat_volc,lon_volc,l
250        write(*,*)trim(rname)//' latitude:',lat_volc
251        write(*,*)trim(rname)//' longitude:',lon_volc
252        latlonvals_total(volci,1) = lat_volc
253        latlonvals_total(volci,2) = lon_volc
254        heightvals_total(volci) = l
255        !Read in number of active tracers being outgassed
256        READ(407,*,iostat=ierr) nvoltrac
257        write(*,*)trim(rname)//' number of outgased tracers:',nvoltrac
258        !Read in eruption frequency (in multiples of iphysiq) and offset of eruption
259        READ(407,'(A)',iostat=ierr) tracline
260        READ(tracline,*) eruptfreq, eruptoffset
261        write(*,*)trim(rname)//' eruption frequency (sols):',eruptfreq
262        write(*,*)trim(rname)//' eruption offset (sols):',eruptoffset
263        eruptfreqs_total(volci,1) = eruptfreq
264        eruptfreqs_total(volci,2) = eruptoffset
265       
266        do ivoltrac=1,nvoltrac
267
268          !Read in tracer name and check it's in traceur.def
269          read(407,'(A)') tracline
270          blank = index(tracline,' ') ! Find position of 1st blank in tracline
271          voltrac(ivoltrac) = tracline(1:blank)
272          ! check it is a valid tracer:
273          if (is_known_tracer(voltrac(ivoltrac))) then
274            !Read in flux of tracer out of volcano (kg/s)
275            READ(407,'(A)',iostat=ierr) tracline
276            READ(tracline,*) volflux
277            iq=tracer_index(voltrac(ivoltrac))
278            volfluxes_total(volci,iq) = volflux
279          else
280            write(*,*) trim(rname)//": problem with ",trim(voltrac(ivoltrac))
281            call abort_physic(rname, 'tracer '//&
282                              trim(voltrac(ivoltrac))//' not found', 1)
283          endif
284         
285        enddo ! of do ivoltrac=1,nvoltrac
286       
287      enddo ! of do while (volci.lt.nvolc_total)
288   
289      close(407)
290    endif ! of if (is_master)
291   
292    ! 2. Share data with all cores
293    call bcast(latlonvals_total)
294    call bcast(heightvals_total)
295    call bcast(eruptfreqs_total)
296    call bcast(volfluxes_total)
297 
298    ! 3. Each core extracts and stores data relevant to its domain only
299    ! 3.1. find out how many volcanoes for this core
300    nvolc=0
301    allocate(belongs_here(nvolc_total),stat=ierr)
302    if (ierr/=0) then
303      write(*,*) trim(rname)//" failed allocation for belong_here()"
304      call abort_physic(rname,"failled allocation",1)
305    endif
306    allocate(volc_index(nvolc_total),stat=ierr)
307    if (ierr/=0) then
308      write(*,*) trim(rname)//" failed allocation for volc_index()"
309      call abort_physic(rname,"failled allocation",1)
310    endif
311   
312    do volci=1,nvolc_total
313      belongs_here(volci)=.false.
314      volc_index(volci)=0
315     
316      lat_volc=latlonvals_total(volci,1)
317      lon_volc=latlonvals_total(volci,2)
318
319      ! find out if this volcano is in the ig cell
320      do ig=1,ngrid
321        boundslon_deg(:)=boundslon(ig,:)/pi*180.
322        boundslat_deg(:)=boundslat(ig,:)/pi*180.
323       
324        minlon=minval(boundslon_deg)
325        maxlon=maxval(boundslon_deg)
326        minlat=minval(boundslat_deg)
327        maxlat=maxval(boundslat_deg)
328       
329        if ((minlat<=lat_volc).and.(lat_volc<=maxlat).and. &
330            (((minlon<=lon_volc).and.(lon_volc<=maxlon)).or. &
331             ((minlon+360.<=lon_volc).and.(lon_volc<=maxlon+360.)).or. &
332             ((minlon-360.<=lon_volc).and.(lon_volc<=maxlon-360.)))) then
333          nvolc=nvolc+1
334          belongs_here(volci)=.true.
335          volc_index(volci)=ig
336          ! the job is done move on to next volcano
337          exit
338        endif
339      enddo ! of do ig=1,ngrid
340    enddo ! of volci=1,nvolc_total
341   
342    ! 3.2. allocate arrays, now that nvolc is known
343    if (nvolc>0) then
344      allocate(latlonvals(nvolc,2),stat=ierr)
345      if (ierr/=0) then
346        write(*,*) trim(rname)//" failed allocation for latlonvals()"
347        call abort_physic(rname,"failled allocation",1)
348      endif
349      allocate(heightvals(nvolc),stat=ierr)
350      if (ierr/=0) then
351        write(*,*) trim(rname)//" failed allocation for heightvals()"
352        call abort_physic(rname,"failled allocation",1)
353      endif
354      allocate(eruptfreqs(nvolc,2),stat=ierr)
355      if (ierr/=0) then
356        write(*,*) trim(rname)//" failed allocation for eruptfreqs()"
357        call abort_physic(rname,"failled allocation",1)
358      endif
359      allocate(volfluxes(nvolc,nq),stat=ierr)
360      if (ierr/=0) then
361        write(*,*) trim(rname)//" failed allocation for volfluxes()"
362        call abort_physic(rname,"failled allocation",1)
363      endif
364      allocate(ivolc(nvolc),stat=ierr)
365      if (ierr/=0) then
366        write(*,*) trim(rname)//" failed allocation for ivolc()"
367        call abort_physic(rname,"failled allocation",1)
368      endif
369    endif ! of if (nvolc>0)
370   
371    ! 3.3. copy over global data to local arrays
372    do volci=1,nvolc
373      do i=1,nvolc_total
374        if (belongs_here(i)) then
375          ! copy data to local arrays
376          latlonvals(volci,1:2)=latlonvals_total(i,1:2)
377          heightvals(volci)=heightvals_total(i)
378          eruptfreqs(volci,1:2)=eruptfreqs_total(i,1:2)
379          volfluxes(volci,1:nq)=volfluxes_total(i,1:nq)
380          ivolc(volci)=volc_index(i)
381          ! the job is done, move on to next one
382          belongs_here(i)=.false.
383          exit
384        endif ! of if (belongs_here(i))
385      enddo ! of do i=1,nvolc_total
386    enddo ! of do volci=1,nvolc
387   
388    ! 4. Sanity check
389    ! verify that the sum of volcanoes over all cores is indeed nvolc_total
390    if (is_master) nvolc_check=0
391    call bcast(nvolc_check)
392    call reduce_sum(nvolc,nvolc_check)
393    call bcast(nvolc_check)
394    if (nvolc_check/=nvolc_total) then
395      write(*,*) trim(rname)//" error. Sum of volcanoes over all cores ", &
396       nvolc_check," differs from nvolc_total=",nvolc_total
397      call abort_physic(rname," error. Check and fix volcano.F90",1)
398    endif
399   
400    ! give a summary, for each core:
401    write(*,*) trim(rname)//": nvolc=",nvolc," /",nvolc_total
402    do volci=1,nvolc
403      write(*,*) trim(rname)//": latlonvals(1:2)=",latlonvals(volci,1:2)
404      write(*,*) trim(rname)//": heightvals=",heightvals(volci)
405      write(*,*) trim(rname)//": eruptfreqs(1:2)=",eruptfreqs(volci,1:2)
406      write(*,*) trim(rname)//": volfluxes(1:nq)=",volfluxes(volci,1:nq)
407      write(*,*) trim(rname)//": ivolc=",ivolc(volci)
408      write(*,*) trim(rname)//": longitude(ivolc)=",longitude_deg(ivolc(volci))
409      write(*,*) trim(rname)//": latitude(ivolc)=",latitude_deg(ivolc(volci))
410    enddo
411   
412  END SUBROUTINE read_volcano
413
414END MODULE volcano_mod
Note: See TracBrowser for help on using the repository browser.