Context Navigation

ch4n2o_correction_mod.F90 @ 5906

Last change on this file since 5906 was 5906, checked in by Laurent Fairhead, 38 hours ago
Forgotten export of variable (strange as it compiled with ifort on irene rome) LF
File size: 15.2 KB

Line
1	!=======================================================================
2	! ch4n2o_correction_mod.F90 - Real-world ==> effective CH4 & N2O
3	!-----------------------------------------------------------------------
4	!
5	! Author: Patricia Cadule
6	!
7	! Concept:
8	! Melnikova, I., Ciais, P., Tanaka, K., Shiogama, H., Tachiiri, K., Yokohata, T., and Boucher, O.:
9	! Carbon cycle and climate feedbacks under CO2 and non-CO2 overshoot pathways,
10	! Earth Syst. Dynam., 16, 257–273,
11	! https://doi.org/10.5194/esd-16-257-2025, 2025
12	!
13	! Calibration inputs:
14	! Olivier Boucher, Jean-Louis Dufresne
15	!
16	! Purpose
17	! -------
18	! Provide effective methane (CH4) and nitrous oxide (N2O)
19	! concentrations to the radiation scheme such that the model's
20	! radiative forcing is consistent with:
21	!
22	! * Melnikova et al. (2025, ESD), Eq. (A4), Table A2 ("m2025"), and
23	! * the Etminan (2016, GRL) based optimisation by O. Boucher &
24	! J.-L. Dufresne ("mbd")
25	!
26	! Both families are implemented here in the same direction:
27	!
28	! C_real (input4MIP / RFMIP) [ppb]
29	! \| (transfer function)
30	! v
31	! C_eff (model-effective) [ppb]
32	!
33	! Modes
34	! -----
35	! ghg_mode = "m2025"
36	! Melnikova et al. (2025), inverse of Eq. (A4):
37	!
38	! CH4: C_eff = CH4_PI + [ (C_real − CH4_PI)/A_CH4_M25 ]^(1/C_CH4_M25)
39	! N2O: C_eff = N2O_PI + [ (C_real − N2O_PI)/B_N2O_M25 ]^(1/D_N2O_M25)
40	!
41	! For C_real <= PI we use the identity mapping:
42	! C_eff = C_real
43	!
44	! ghg_mode = "mbd"
45	! Melnikova-Boucher–Dufresne mapping:
46	!
47	! CH4_mod2rw:
48	! C_real = CH4_PI
49	! + A_CH4_MBD * (C_eff − CH4_MIN)^P_CH4_MBD
50	! − A_CH4_MBD * (CH4_PI − CH4_MIN)^P_CH4_MBD
51	!
52	! Inverse (used here):
53	! inner = (C_real − CH4_PI
54	! + A_CH4_MBD * (CH4_PI − CH4_MIN)^P_CH4_MBD) / A_CH4_MBD
55	! if inner > 0:
56	! C_eff = CH4_MIN + inner^(1/P_CH4_MBD)
57	! else:
58	! C_eff = CH4_MIN ! saturate at lower effective bound
59	!
60	! N2O is analogous, with N2O_MIN, A_N2O_MBD, P_N2O_MBD.
61	!
62	! The valid real-world domain is C_real >= C_min_valid_rw, where
63	! C_min_valid_rw = C_PI - A*(C_PI - C_min)^P. Below this the
64	! analytic inverse is not defined; in that regime we clip to
65	! C_eff = C_min and optionally print a single debug warning.
66	!
67	! ghg_mode = "identity"
68	! Pure pass-through for diagnostics:
69	! C_eff = C_real
70	!
71	! Notes
72	! -------------------
73	! * For very deep sub-PI values one may enforce a floor in effective
74	! concentrations via:
75	! ghg_min_eff_ppb_ch4, ghg_min_eff_ppb_n2o
76	!
77	! Public API
78	! ----------
79	! CALL ch4n2o_init()
80	! CALL ch4n2o_compute(year, ch4_in_ppb, n2o_in_ppb, ch4_eff_ppb, n2o_eff_ppb)
81	! CALL ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
82	! CALL n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)
83	!
84	!=======================================================================
85	MODULE ch4n2o_correction_mod
86	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
87	USE mod_phys_lmdz_para, ONLY: bcast
88	USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
89	USE ioipsl_getin_p_mod, ONLY: getin_p
90	IMPLICIT NONE
91	PRIVATE
92
93	INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 300)
94
95	!---------------- Public routines ------------------------------------
96	PUBLIC :: ch4n2o_init
97	PUBLIC :: ch4n2o_compute
98	PUBLIC :: ch4_eff_compute_only
99	PUBLIC :: n2o_eff_compute_only
100
101	! Optional routing of effective ppb to outputs elsewhere
102	LOGICAL, PUBLIC :: ok_CH4_eff_ppb = .FALSE.
103	!$OMP THREADPRIVATE(ok_CH4_eff_ppb)
104	LOGICAL, PUBLIC :: ok_N2O_eff_ppb = .FALSE.
105	!$OMP THREADPRIVATE(ok_N2O_eff_ppb)
106
107	INTEGER, PUBLIC :: ghg_year = 1850 ! for diagnostics only
108	CHARACTER(LEN=16), PUBLIC :: ghg_mode = 'm2025' ! 'm2025' \| 'mbd' \| 'identity'
109
110	!---------------- Configuration (defaults) ---------------------------
111
112	LOGICAL :: ghg_debug = .FALSE.
113
114	! Pre-industrial pivots
115	REAL(dp) :: ghg_ch4_zero_ppb = 808.25_dp
116	REAL(dp) :: ghg_n2o_zero_ppb = 273.02_dp
117
118	! Melnikova et al. (2025) parameters (Table A2)
119	REAL(dp) :: ghg_A_ch4_m25 = 19.44701978_dp ! A_CH4_M25
120	REAL(dp) :: ghg_C_ch4_m25 = 0.49593024_dp ! C_CH4_M25
121	REAL(dp) :: ghg_B_n2o_m25 = 0.84856644_dp ! B_N2O_M25
122	REAL(dp) :: ghg_D_n2o_m25 = 1.13865386_dp ! D_N2O_M25
123
124	! Boucher–Dufresne parameters
125	REAL(dp) :: ghg_ch4_min_ppb = 808.25_dp - 100.0_dp ! CH4_MIN
126	REAL(dp) :: ghg_n2o_min_ppb = 273.02_dp - 80.0_dp ! N2O_MIN
127
128	REAL(dp) :: ghg_A_ch4_mbd = 77.74824074_dp ! A_CH4_MBD
129	REAL(dp) :: ghg_P_ch4_mbd = 0.36559975_dp ! P_CH4_MBD
130	REAL(dp) :: ghg_A_n2o_mbd = 0.17792515_dp ! A_N2O_MBD
131	REAL(dp) :: ghg_P_n2o_mbd = 1.38444104_dp ! P_N2O_MBD
132
133	! Real-world minima for which the MBD inverse is defined
134	REAL(dp) :: ghg_ch4_min_valid_rw = 0.0_dp
135	REAL(dp) :: ghg_n2o_min_valid_rw = 0.0_dp
136
137	! Optional floors for extreme sub-PI robustness
138	REAL(dp) :: ghg_min_eff_ppb_ch4 = 0.0_dp
139	REAL(dp) :: ghg_min_eff_ppb_n2o = 0.0_dp
140
141	! One-time flags to avoid flooding logs in debug mode
142	LOGICAL :: warned_ch4_mbd_below = .FALSE.
143	LOGICAL :: warned_n2o_mbd_below = .FALSE.
144
145	CONTAINS
146	!=======================================================================
147	SUBROUTINE ch4n2o_init()
148	CHARACTER(LEN=32) :: mode_in
149	REAL(dp) :: rbuf(18)
150	INTEGER :: ibuf(4), mode_code
151
152	mode_in = ghg_mode
153	CALL getin_p('ghg_mode', mode_in)
154	CALL getin_p('ghg_debug', ghg_debug)
155	CALL getin_p('ghg_year', ghg_year)
156
157	! Allow overriding of parameters if needed
158	CALL getin_p('ghg_ch4_zero_ppb', ghg_ch4_zero_ppb)
159	CALL getin_p('ghg_n2o_zero_ppb', ghg_n2o_zero_ppb)
160
161	CALL getin_p('ghg_A_ch4_m25', ghg_A_ch4_m25)
162	CALL getin_p('ghg_C_ch4_m25', ghg_C_ch4_m25)
163	CALL getin_p('ghg_B_n2o_m25', ghg_B_n2o_m25)
164	CALL getin_p('ghg_D_n2o_m25', ghg_D_n2o_m25)
165
166	CALL getin_p('ghg_ch4_min_ppb', ghg_ch4_min_ppb)
167	CALL getin_p('ghg_n2o_min_ppb', ghg_n2o_min_ppb)
168	CALL getin_p('ghg_A_ch4_mbd', ghg_A_ch4_mbd)
169	CALL getin_p('ghg_P_ch4_mbd', ghg_P_ch4_mbd)
170	CALL getin_p('ghg_A_n2o_mbd', ghg_A_n2o_mbd)
171	CALL getin_p('ghg_P_n2o_mbd', ghg_P_n2o_mbd)
172
173	CALL getin_p('ghg_min_eff_ppb_ch4', ghg_min_eff_ppb_ch4)
174	CALL getin_p('ghg_min_eff_ppb_n2o', ghg_min_eff_ppb_n2o)
175
176	CALL to_lower_inplace(mode_in)
177	SELECT CASE (TRIM(mode_in))
178	CASE ('m2025','mbd','identity')
179	ghg_mode = TRIM(mode_in)
180	CASE DEFAULT
181	IF (is_mpi_root .AND. is_omp_root) THEN
182	WRITE(,) '[ch4n2o] unknown ghg_mode="',TRIM(mode_in),'" => using m2025'
183	END IF
184	ghg_mode = 'm2025'
185	END SELECT
186
187	! Broadcast scalars across MPI tasks
188	rbuf = (/ ghg_ch4_zero_ppb, ghg_n2o_zero_ppb, &
189	ghg_A_ch4_m25, ghg_C_ch4_m25, ghg_B_n2o_m25, ghg_D_n2o_m25, &
190	ghg_ch4_min_ppb, ghg_n2o_min_ppb, &
191	ghg_A_ch4_mbd, ghg_P_ch4_mbd, ghg_A_n2o_mbd, ghg_P_n2o_mbd, &
192	ghg_min_eff_ppb_ch4, ghg_min_eff_ppb_n2o, &
193	0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp /)
194	CALL bcast(rbuf)
195
196	ghg_ch4_zero_ppb = rbuf(1)
197	ghg_n2o_zero_ppb = rbuf(2)
198	ghg_A_ch4_m25 = rbuf(3)
199	ghg_C_ch4_m25 = rbuf(4)
200	ghg_B_n2o_m25 = rbuf(5)
201	ghg_D_n2o_m25 = rbuf(6)
202	ghg_ch4_min_ppb = rbuf(7)
203	ghg_n2o_min_ppb = rbuf(8)
204	ghg_A_ch4_mbd = rbuf(9)
205	ghg_P_ch4_mbd = rbuf(10)
206	ghg_A_n2o_mbd = rbuf(11)
207	ghg_P_n2o_mbd = rbuf(12)
208	ghg_min_eff_ppb_ch4 = rbuf(13)
209	ghg_min_eff_ppb_n2o = rbuf(14)
210
211	ibuf(1) = MERGE(1,0,ghg_debug)
212	ibuf(2) = ghg_year
213	ibuf(3) = 0
214	ibuf(4) = 0
215	CALL bcast(ibuf)
216	ghg_debug = (ibuf(1) == 1)
217	ghg_year = ibuf(2)
218
219	! Encode mode as integer for robustness
220	SELECT CASE (TRIM(ghg_mode))
221	CASE ('m2025'); mode_code = 1
222	CASE ('mbd'); mode_code = 2
223	CASE ('identity'); mode_code = 3
224	END SELECT
225	CALL bcast(mode_code)
226	SELECT CASE (mode_code)
227	CASE (1); ghg_mode = 'm2025'
228	CASE (2); ghg_mode = 'mbd'
229	CASE (3); ghg_mode = 'identity'
230	CASE DEFAULT
231	ghg_mode = 'm2025'
232	END SELECT
233
234	! Compute minimum real-world values for which MBD inverse is defined
235	ghg_ch4_min_valid_rw = ghg_ch4_zero_ppb - &
236	ghg_A_ch4_mbd * (ghg_ch4_zero_ppb - ghg_ch4_min_ppb)**ghg_P_ch4_mbd
237	ghg_n2o_min_valid_rw = ghg_n2o_zero_ppb - &
238	ghg_A_n2o_mbd * (ghg_n2o_zero_ppb - ghg_n2o_min_ppb)**ghg_P_n2o_mbd
239
240	IF (is_mpi_root .AND. is_omp_root) THEN
241	WRITE(*,'(A,1X,A)') '[ch4n2o] mode =',TRIM(ghg_mode)
242	WRITE(*,'(A,2(F9.2,1X))') ' pivots (ppb) CH4,N2O =', &
243	ghg_ch4_zero_ppb, ghg_n2o_zero_ppb
244	SELECT CASE (TRIM(ghg_mode))
245	CASE ('m2025')
246	WRITE(*,'(A,2(F10.5,1X))') ' CH4 A,C =',ghg_A_ch4_m25,ghg_C_ch4_m25
247	WRITE(*,'(A,2(F10.5,1X))') ' N2O B,D =',ghg_B_n2o_m25,ghg_D_n2o_m25
248	CASE ('mbd')
249	WRITE(*,'(A,2(F10.6,1X))') ' CH4 A,P =',ghg_A_ch4_mbd,ghg_P_ch4_mbd
250	WRITE(*,'(A,2(F10.6,1X))') ' N2O A,P =',ghg_A_n2o_mbd,ghg_P_n2o_mbd
251	WRITE(*,'(A,2(F9.2,1X))') ' mins (ppb) CH4,N2O =', &
252	ghg_ch4_min_ppb, ghg_n2o_min_ppb
253	WRITE(*,'(A,2(F9.2,1X))') ' valid C_real >= (CH4,N2O) =', &
254	ghg_ch4_min_valid_rw, ghg_n2o_min_valid_rw
255	END SELECT
256	WRITE(*,'(A,I6)') ' ghg_year =', ghg_year
257	END IF
258	END SUBROUTINE ch4n2o_init
259	!-----------------------------------------------------------------------
260
261	SUBROUTINE ch4n2o_compute(year, ch4_in_ppb, n2o_in_ppb, ch4_eff_ppb, n2o_eff_ppb)
262	INTEGER, INTENT(IN) :: year
263	REAL(dp), INTENT(IN) :: ch4_in_ppb, n2o_in_ppb
264	REAL(dp), INTENT(OUT) :: ch4_eff_ppb, n2o_eff_ppb
265
266	CALL ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
267	CALL n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)
268
269	IF (ghg_debug .AND. is_mpi_root .AND. is_omp_root) THEN
270	WRITE(*,'(A,I6,2(A,F12.3))') 'CH4: y=',year,' real=',ch4_in_ppb,' eff=',ch4_eff_ppb
271	WRITE(*,'(A,I6,2(A,F12.3))') 'N2O: y=',year,' real=',n2o_in_ppb,' eff=',n2o_eff_ppb
272	END IF
273	END SUBROUTINE ch4n2o_compute
274	!-----------------------------------------------------------------------
275
276	SUBROUTINE ch4_eff_compute_only(year, ch4_in_ppb, ch4_eff_ppb)
277	INTEGER, INTENT(IN) :: year
278	REAL(dp), INTENT(IN) :: ch4_in_ppb
279	REAL(dp), INTENT(OUT) :: ch4_eff_ppb
280
281	SELECT CASE (TRIM(ghg_mode))
282	CASE ('identity')
283	ch4_eff_ppb = ch4_in_ppb
284
285	CASE ('m2025')
286	ch4_eff_ppb = ceff_ch4_m2025(ch4_in_ppb, ghg_ch4_zero_ppb, &
287	ghg_A_ch4_m25, ghg_C_ch4_m25)
288
289	CASE ('mbd')
290	ch4_eff_ppb = ceff_ch4_mbd(ch4_in_ppb, ghg_ch4_zero_ppb, ghg_ch4_min_ppb, &
291	ghg_A_ch4_mbd, ghg_P_ch4_mbd)
292
293	IF (ghg_debug .AND. .NOT. warned_ch4_mbd_below) THEN
294	IF (ch4_in_ppb < ghg_ch4_min_valid_rw .AND. is_mpi_root .AND. is_omp_root) THEN
295	WRITE(*,'(A,F10.3,A,F10.3,A)') &
296	'[ch4n2o] CH4 MBD: C_real below valid domain (', ch4_in_ppb, &
297	' < ', ghg_ch4_min_valid_rw, ' ppb), C_eff clipped to CH4_MIN.'
298	warned_ch4_mbd_below = .TRUE.
299	END IF
300	END IF
301	END SELECT
302
303	IF (ghg_min_eff_ppb_ch4 > 0.0_dp) THEN
304	ch4_eff_ppb = MAX(ch4_eff_ppb, ghg_min_eff_ppb_ch4)
305	END IF
306	END SUBROUTINE ch4_eff_compute_only
307	!-----------------------------------------------------------------------
308
309	SUBROUTINE n2o_eff_compute_only(year, n2o_in_ppb, n2o_eff_ppb)
310	INTEGER, INTENT(IN) :: year
311	REAL(dp), INTENT(IN) :: n2o_in_ppb
312	REAL(dp), INTENT(OUT) :: n2o_eff_ppb
313
314	SELECT CASE (TRIM(ghg_mode))
315	CASE ('identity')
316	n2o_eff_ppb = n2o_in_ppb
317
318	CASE ('m2025')
319	n2o_eff_ppb = ceff_n2o_m2025(n2o_in_ppb, ghg_n2o_zero_ppb, &
320	ghg_B_n2o_m25, ghg_D_n2o_m25)
321
322	CASE ('mbd')
323	n2o_eff_ppb = ceff_n2o_mbd(n2o_in_ppb, ghg_n2o_zero_ppb, ghg_n2o_min_ppb, &
324	ghg_A_n2o_mbd, ghg_P_n2o_mbd)
325
326	IF (ghg_debug .AND. .NOT. warned_n2o_mbd_below) THEN
327	IF (n2o_in_ppb < ghg_n2o_min_valid_rw .AND. is_mpi_root .AND. is_omp_root) THEN
328	WRITE(*,'(A,F10.3,A,F10.3,A)') &
329	'[ch4n2o] N2O MBD: C_real below valid domain (', n2o_in_ppb, &
330	' < ', ghg_n2o_min_valid_rw, ' ppb), C_eff clipped to N2O_MIN.'
331	warned_n2o_mbd_below = .TRUE.
332	END IF
333	END IF
334	END SELECT
335
336	IF (ghg_min_eff_ppb_n2o > 0.0_dp) THEN
337	n2o_eff_ppb = MAX(n2o_eff_ppb, ghg_min_eff_ppb_n2o)
338	END IF
339	END SUBROUTINE n2o_eff_compute_only
340
341	!-----------------------------------------------------------------------
342	! Maths kernels (PURE ELEMENTAL)
343	!-----------------------------------------------------------------------
344
345	PURE ELEMENTAL FUNCTION ceff_ch4_m2025(c_real, c_pi, A, Cexp) RESULT(c_eff)
346	REAL(dp), INTENT(IN) :: c_real, c_pi, A, Cexp
347	REAL(dp) :: c_eff, delta
348	#ifdef USE_OMP_SIMD
349	!$OMP DECLARE SIMD(ceff_ch4_m2025)
350	#endif
351	IF (c_real <= c_pi) THEN
352	c_eff = c_real
353	ELSE
354	delta = (c_real - c_pi) / A
355	c_eff = c_pi + delta**(1.0_dp/Cexp)
356	END IF
357	END FUNCTION ceff_ch4_m2025
358	!-----------------------------------------------------------------------
359
360	PURE ELEMENTAL FUNCTION ceff_n2o_m2025(c_real, c_pi, B, Dexp) RESULT(c_eff)
361	REAL(dp), INTENT(IN) :: c_real, c_pi, B, Dexp
362	REAL(dp) :: c_eff, delta
363	#ifdef USE_OMP_SIMD
364	!$OMP DECLARE SIMD(ceff_n2o_m2025)
365	#endif
366	IF (c_real <= c_pi) THEN
367	c_eff = c_real
368	ELSE
369	delta = (c_real - c_pi) / B
370	c_eff = c_pi + delta**(1.0_dp/Dexp)
371	END IF
372	END FUNCTION ceff_n2o_m2025
373	!-----------------------------------------------------------------------
374
375	PURE ELEMENTAL FUNCTION ceff_ch4_mbd(c_real, c_pi, c_min, A, Pexp) RESULT(c_eff)
376	! Boucher–Dufresne CH4 inverse mapping
377	REAL(dp), INTENT(IN) :: c_real, c_pi, c_min, A, Pexp
378	REAL(dp) :: c_eff, inner
379	#ifdef USE_OMP_SIMD
380	!$OMP DECLARE SIMD(ceff_ch4_mbd)
381	#endif
382	inner = (c_real - c_pi + A(c_pi - c_min)*Pexp) / A
383	IF (inner > 0.0_dp) THEN
384	c_eff = c_min + inner**(1.0_dp/Pexp)
385	ELSE
386	c_eff = c_min
387	END IF
388	END FUNCTION ceff_ch4_mbd
389	!-----------------------------------------------------------------------
390
391	PURE ELEMENTAL FUNCTION ceff_n2o_mbd(c_real, c_pi, c_min, A, Pexp) RESULT(c_eff)
392	! Boucher–Dufresne N2O inverse mapping
393	REAL(dp), INTENT(IN) :: c_real, c_pi, c_min, A, Pexp
394	REAL(dp) :: c_eff, inner
395	#ifdef USE_OMP_SIMD
396	!$OMP DECLARE SIMD(ceff_n2o_mbd)
397	#endif
398	inner = (c_real - c_pi + A(c_pi - c_min)*Pexp) / A
399	IF (inner > 0.0_dp) THEN
400	c_eff = c_min + inner**(1.0_dp/Pexp)
401	ELSE
402	c_eff = c_min
403	END IF
404	END FUNCTION ceff_n2o_mbd
405
406	!-----------------------------------------------------------------------
407	! Utility: in-place lower-casing for mode strings
408	!-----------------------------------------------------------------------
409	SUBROUTINE to_lower_inplace(s)
410	CHARACTER(*), INTENT(INOUT) :: s
411	INTEGER :: i, ia
412	DO i = 1, LEN_TRIM(s)
413	ia = IACHAR(s(i:i))
414	IF (ia >= IACHAR('A') .AND. ia <= IACHAR('Z')) s(i:i) = ACHAR(ia+32)
415	END DO
416	END SUBROUTINE to_lower_inplace
417
418	END MODULE ch4n2o_correction_mod
419

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format