source: LMDZ6/trunk/libf/phylmd/StratAer/strataer_mod.F90 @ 4514

Last change on this file since 4514 was 4513, checked in by tlurton, 20 months ago

Strataer: a bug correction in strataer_mod.F90, and an implementation of injection duration for SAI case.
  • Property svn:keywords set to Id
File size: 10.8 KB
Line 
1! $Id: strataer_mod.F90 4513 2023-04-20 07:57:22Z oboucher $
2MODULE strataer_mod
3! This module contains information about strato microphysic model parameters
4 
5  IMPLICIT NONE
6
7  ! flag to constraint nucleation rate in a lat/pres box
8  LOGICAL,SAVE :: flag_nuc_rate_box      ! Nucleation rate limit or not to a lat/pres
9  !$OMP THREADPRIVATE(flag_nuc_rate_box)
10  REAL,SAVE    :: nuclat_min             ! min lat to activate nuc rate
11  REAL,SAVE    :: nuclat_max             ! max lat to activate nuc rate
12  REAL,SAVE    :: nucpres_min            ! min pres to activate nuc rate
13  REAL,SAVE    :: nucpres_max            ! max pres to activate nuc rate
14  !$OMP THREADPRIVATE(nuclat_min, nuclat_max, nucpres_min, nucpres_max)
15 
16  ! flag for sulfur emission scenario: (0) background aer ; (1) volcanic eruption ; (2) strato aer injections (SAI)
17  INTEGER,SAVE :: flag_sulf_emit
18  !$OMP THREADPRIVATE(flag_sulf_emit)
19
20  ! flag for sulfur emission altitude distribution: (0) gaussian; (1) uniform
21  INTEGER,SAVE :: flag_sulf_emit_distrib
22  !$OMP THREADPRIVATE(flag_sulf_emit_distrib)
23 
24  !--flag_sulf_emit=1 -- Volcanic eruption(s)
25  INTEGER,SAVE             :: nErupt                    ! number of eruptions specs
26  REAL,SAVE                :: injdur                    ! volcanic injection duration
27  !$OMP THREADPRIVATE(nErupt, injdur)
28  INTEGER,ALLOCATABLE,SAVE :: year_emit_vol(:)          ! year of emission date
29  INTEGER,ALLOCATABLE,SAVE :: mth_emit_vol(:)           ! month of emission date
30  INTEGER,ALLOCATABLE,SAVE :: day_emit_vol(:)           ! day of emission date
31  !$OMP THREADPRIVATE(year_emit_vol, mth_emit_vol, day_emit_vol)
32  REAL,ALLOCATABLE,SAVE    :: m_aer_emiss_vol(:)        ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
33  REAL,ALLOCATABLE,SAVE    :: altemiss_vol(:)           ! emission altitude in m
34  REAL,ALLOCATABLE,SAVE    :: sigma_alt_vol(:)          ! standard deviation of emission altitude in m
35  !$OMP THREADPRIVATE(m_aer_emiss_vol, altemiss_vol, sigma_alt_vol)
36  INTEGER,ALLOCATABLE,SAVE :: ponde_lonlat_vol(:)       ! lon/lat ponderation factor
37  REAL,ALLOCATABLE,SAVE    :: xlat_min_vol(:)           ! min latitude of volcano in degree
38  REAL,ALLOCATABLE,SAVE    :: xlat_max_vol(:)           ! max latitude of volcano in degree
39  REAL,ALLOCATABLE,SAVE    :: xlon_min_vol(:)           ! min longitude of volcano in degree
40  REAL,ALLOCATABLE,SAVE    :: xlon_max_vol(:)           ! max longitude of volcano in degree
41  !$OMP THREADPRIVATE(ponde_lonlat_vol, xlat_min_vol, xlat_max_vol, xlon_min_vol, xlon_max_vol)
42 
43  !--flag_sulf_emit=2 --SAI
44  REAL,SAVE    :: m_aer_emiss_sai        ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
45  REAL,SAVE    :: altemiss_sai           ! emission altitude in m
46  REAL,SAVE    :: sigma_alt_sai          ! standard deviation of emission altitude in m
47  !$OMP THREADPRIVATE(m_aer_emiss_sai, altemiss_sai, sigma_alt_sai)
48  REAL,SAVE    :: xlat_sai               ! latitude of SAI in degree
49  REAL,SAVE    :: xlon_sai               ! longitude of SAI in degree
50  REAL,SAVE    :: dlat, dlon             ! delta latitude and d longitude of grid in degree
51  !$OMP THREADPRIVATE(xlat_sai, xlon_sai, dlat, dlon)
52 
53  !--flag_sulf_emit=3 -- SAI
54  REAL,SAVE    :: xlat_max_sai           ! maximum latitude of SAI in degrees
55  REAL,SAVE    :: xlat_min_sai           ! minimum latitude of SAI in degrees
56  !$OMP THREADPRIVATE(xlat_min_sai,xlat_max_sai)
57 
58CONTAINS
59   
60  SUBROUTINE strataer_init()
61    USE ioipsl_getin_p_mod, ONLY : getin_p
62    USE print_control_mod, ONLY : lunout
63    USE mod_phys_lmdz_para, ONLY : is_master
64
65    ! Local var
66    INTEGER       :: ieru
67 
68    INTEGER :: i
69   
70    WRITE(lunout,*) 'IN STRATAER INIT WELCOME!'
71   
72    !Config Key  = flag_sulf_emit
73    !Config Desc = aerosol emission mode
74    ! - 0 = background aerosol
75    ! - 1 = volcanic eruption
76    ! - 2 = geo-ingeneering design
77    ! - 3 = geo-engineering between two latitudes
78    !Config Def  = 0
79    !Config Help = Used in physiq.F
80    !
81    flag_sulf_emit = 0 ! Background (default)
82    flag_sulf_emit_distrib = 0 ! Gaussian (default)
83    nErupt = 0 ! eruption number
84    injdur = 0 ! init injection duration
85    CALL getin_p('flag_sulf_emit',flag_sulf_emit)
86
87    IF (flag_sulf_emit==1) THEN ! Volcano
88       CALL getin_p('nErupt',nErupt)
89       CALL getin_p('injdur',injdur)
90       CALL getin_p('flag_sulf_emit_distrib',flag_sulf_emit_distrib)
91
92       year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0
93       m_aer_emiss_vol=0. ; altemiss_vol=0. ; sigma_alt_vol=0.
94       xlon_min_vol=0. ; xlon_max_vol=0.
95       xlat_min_vol=0. ; xlat_max_vol=0.
96
97       IF (nErupt.GT.0) THEN
98          ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt))
99          ALLOCATE(m_aer_emiss_vol(nErupt),altemiss_vol(nErupt),sigma_alt_vol(nErupt))
100          ALLOCATE(xlat_min_vol(nErupt),xlon_min_vol(nErupt))
101          ALLOCATE(xlat_max_vol(nErupt),xlon_max_vol(nErupt))
102       ELSE
103          WRITE(lunout,*) 'ERROR : Using flag_sulf_emit=1 (ie Volcanic eruption) but nErupt (',nErupt,') <=0 !'
104          CALL abort_physic('strataer_mod','No eruption define in physiq_def (nErupt=0). Add one or use background condition.',1)
105       ENDIF
106
107       CALL getin_p('year_emit_vol',year_emit_vol)
108       CALL getin_p('mth_emit_vol',mth_emit_vol)
109       CALL getin_p('day_emit_vol',day_emit_vol)
110       CALL getin_p('m_aer_emiss_vol',m_aer_emiss_vol)
111       CALL getin_p('altemiss_vol',altemiss_vol)
112       CALL getin_p('sigma_alt_vol',sigma_alt_vol)
113       CALL getin_p('xlon_min_vol',xlon_min_vol)
114       CALL getin_p('xlon_max_vol',xlon_max_vol)
115       CALL getin_p('xlat_min_vol',xlat_min_vol)
116       CALL getin_p('xlat_max_vol',xlat_max_vol)
117    ELSEIF (flag_sulf_emit == 2) THEN ! SAI
118       CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
119       CALL getin_p('altemiss_sai',altemiss_sai)
120       CALL getin_p('sigma_alt_sai',sigma_alt_sai)
121       CALL getin_p('xlat_sai',xlat_sai)
122       CALL getin_p('xlon_sai',xlon_sai)
123    ELSEIF (flag_sulf_emit == 3) THEN ! SAI between latitudes
124       CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
125       CALL getin_p('altemiss_sai',altemiss_sai)
126       CALL getin_p('sigma_alt_sai',sigma_alt_sai)
127       CALL getin_p('xlon_sai',xlon_sai)
128       CALL getin_p('xlat_max_sai',xlat_max_sai)
129       CALL getin_p('xlat_min_sai',xlat_min_sai)
130    ENDIF
131
132    !Config Key  = flag_nuc_rate_box
133    !Config Desc = define or not a box for nucleation rate
134    ! - F = global nucleation
135    ! - T = 2D-box for nucleation need nuclat_min, nuclat_max, nucpres_min and
136    ! nucpres_max
137    !       to define its bounds.
138    !Config Def  = F
139    !Config Help = Used in physiq.F
140    !
141    flag_nuc_rate_box = .FALSE.
142    nuclat_min=0  ; nuclat_max=0
143    nucpres_min=0 ; nucpres_max=0
144
145    CALL getin_p('flag_nuc_rate_box',flag_nuc_rate_box)
146    CALL getin_p('nuclat_min',nuclat_min)
147    CALL getin_p('nuclat_max',nuclat_max)
148    CALL getin_p('nucpres_min',nucpres_min)
149    CALL getin_p('nucpres_max',nucpres_max)
150
151    !IF (is_master) THEN
152    WRITE(lunout,*) 'flag_sulf_emit = ',flag_sulf_emit
153    IF (flag_sulf_emit == 1) THEN
154       WRITE(lunout,*) 'IN STRATAER nErupt: ',nErupt
155       WRITE(lunout,*) 'IN STRATAER injdur: ',injdur
156       WRITE(lunout,*) 'IN STRATAER : year_emit_vol',year_emit_vol
157       WRITE(lunout,*) 'IN STRATAER : mth_emit_vol',mth_emit_vol
158       WRITE(lunout,*) 'IN STRATAER : day_emit_vol',day_emit_vol
159       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_vol',m_aer_emiss_vol
160       WRITE(lunout,*) 'IN STRATAER : altemiss_vol',altemiss_vol
161       WRITE(lunout,*) 'IN STRATAER : sigma_alt_vol',sigma_alt_vol
162       WRITE(lunout,*) 'IN STRATAER : xlat_min_vol',xlat_min_vol
163       WRITE(lunout,*) 'IN STRATAER : xlat_max_vol',xlat_max_vol
164       WRITE(lunout,*) 'IN STRATAER : xlon_min_vol',xlon_min_vol
165       WRITE(lunout,*) 'IN STRATAER : xlon_max_vol',xlon_max_vol
166    ELSEIF (flag_sulf_emit == 2) THEN
167       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
168       WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
169       WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
170       WRITE(lunout,*) 'IN STRATAER : xlat_sai',xlat_sai
171       WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
172    ELSEIF (flag_sulf_emit == 3) THEN
173       WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
174       WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
175       WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
176       WRITE(lunout,*) 'IN STRATAER : xlat_min_sai',xlat_min_sai
177       WRITE(lunout,*) 'IN STRATAER : xlat_max_sai',xlat_max_sai
178       WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
179    ENDIF
180    WRITE(lunout,*) 'flag_sulf_emit_distrib = ',flag_sulf_emit_distrib
181    WRITE(lunout,*) 'IN STRATAER : flag_nuc_rate_box = ',flag_nuc_rate_box
182    IF (flag_nuc_rate_box) THEN
183       WRITE(lunout,*) 'IN STRATAER : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max
184       WRITE(lunout,*) 'IN STRATAER : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max
185    ENDIF
186    !ENDIF
187
188    IF (flag_sulf_emit > 0) THEN
189       CALL strataer_ponde_init
190       WRITE(lunout,*) 'IN STRATAER INT2 END'
191    END IF
192    WRITE(lunout,*) 'IN STRATAER END'
193
194  END SUBROUTINE strataer_init
195 
196  ! Compute the ponderation to applicate in each grid point for all eruptions and init
197  ! dlat & dlon variables
198  SUBROUTINE strataer_ponde_init()
199   
200    USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
201    USE dimphy, ONLY: klon
202    USE mod_grid_phy_lmdz, ONLY: nbp_lat, nbp_lon
203    USE print_control_mod, ONLY : lunout
204
205    INCLUDE "YOMCST.h"  !--RPI
206
207    ! local var
208    REAL                :: pi,lat_reg_deg,lon_reg_deg! latitude and d longitude of grid in degree
209    INTEGER             :: ieru, i, j
210   
211    ALLOCATE(ponde_lonlat_vol(nErupt))
212   
213    !Compute lon/lat ponderation for injection
214    dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
215    dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
216    WRITE(lunout,*) 'IN STRATAER_INIT dlat=',dlat,'dlon=',dlon
217    WRITE(lunout,*) 'IN STRATAER_INIT nErupt=',nErupt
218    WRITE(lunout,*) 'IN STRATAER_INIT xlat_min=',xlat_min_vol,'xlat_max=',xlat_max_vol
219    WRITE(lunout,*) 'IN STRATAER_INIT xlon_min=',xlon_min_vol,'xlon_max=',xlon_max_vol
220    DO ieru=1, nErupt
221       ponde_lonlat_vol(ieru) = 0
222       DO i=1,nbp_lon
223          lon_reg_deg = lon_reg(i)*180./RPI
224          DO j=1,nbp_lat
225             lat_reg_deg = lat_reg(j)*180./RPI
226             IF  ( lat_reg_deg.GE.xlat_min_vol(ieru)-dlat .AND. lat_reg_deg.LT.xlat_max_vol(ieru)+dlat .AND. &
227                  lon_reg_deg.GE.xlon_min_vol(ieru)-dlon .AND. lon_reg_deg.LT.xlon_max_vol(ieru)+dlon ) THEN
228                ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1
229             ENDIF
230          ENDDO
231       ENDDO
232       IF(ponde_lonlat_vol(ieru) == 0) THEN
233          WRITE(lunout,*) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=',ieru
234       ENDIF
235    ENDDO !ieru
236    WRITE(lunout,*) 'IN STRATAER_INIT ponde_lonlat: ',ponde_lonlat_vol
237   
238  END SUBROUTINE strataer_ponde_init
239 
240END MODULE strataer_mod
Note: See TracBrowser for help on using the repository browser.