Context Navigation

interp_sulf_input.F90 @ 5086

Last change on this file since 5086 was 5082, checked in by abarral, 16 months ago
(lint) Fix obsolete boolean operators
Property svn:keywords set to `Id`
File size: 13.2 KB

Line
1	!
2	! $Id: interp_sulf_input.F90 5082 2024-07-19 15:41:58Z abarral $
3	!
4	SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
5
6	USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
7	nf95_inq_varid, nf95_inquire_dimension, nf95_open
8	USE lmdz_netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
9
10	USE mod_grid_phy_lmdz
11	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
12	USE mod_phys_lmdz_omp_data, ONLY : is_omp_root
13	USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2
14	USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim
15	USE mod_phys_lmdz_para
16	USE dimphy
17	USE phys_cal_mod
18	USE infotrac_phy
19	USE aerophys
20	USE YOMCST
21	USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer
22
23	IMPLICIT NONE
24
25	include "dimensions.h"
26
27	! Variable input
28	REAL paprs(klon,klev+1)
29	REAL tr_seri(klon,klev,nbtr)
30	REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
31	LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
32
33	! Variables locales
34	INTEGER n_lat ! number of latitudes in the input data
35	INTEGER n_lon ! number of longitudes in the input data
36	INTEGER, SAVE :: n_lev ! number of levels in the input data
37	!$OMP THREADPRIVATE(n_lev)
38	INTEGER n_mth ! number of months in the input data
39	INTEGER, SAVE :: mth_pre
40	!$OMP THREADPRIVATE(mth_pre)
41
42	! Champs reconstitues
43	REAL paprs_glo(klon_glo,klev+1)
44
45	REAL, ALLOCATABLE :: latitude(:)
46	! (of input data sorted in strictly ascending order)
47
48	REAL, ALLOCATABLE :: longitude(:)
49	! (of input data sorted in strictly ascending order)
50
51	REAL, ALLOCATABLE :: time(:)
52	! (of input data sorted in strictly ascending order)
53
54	REAL, ALLOCATABLE :: lev(:)
55	! levels of input data
56
57	REAL, ALLOCATABLE :: OCS_clim_in(:, :, :, :)
58	REAL, ALLOCATABLE :: SO2_clim_in(:, :, :, :)
59	REAL, ALLOCATABLE :: OCS_clim_mth(:, :, :)
60	REAL, ALLOCATABLE :: SO2_clim_mth(:, :, :)
61	REAL, ALLOCATABLE :: OCS_clim_tmp(:, :)
62	REAL, ALLOCATABLE :: SO2_clim_tmp(:, :)
63	REAL OCS_clim_glo(klon_glo,klev)
64	REAL SO2_clim_glo(klon_glo,klev)
65
66	REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :)
67	REAL, ALLOCATABLE :: O3_clim_mth(:, :, :)
68	REAL, ALLOCATABLE :: O3_clim_tmp(:, :)
69	REAL O3_clim_glo(klon_glo,klev)
70
71	! fixed climatos
72	REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
73	REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
74	REAL, ALLOCATABLE :: OCS_lifetime_mth(:, :, :)
75	REAL, ALLOCATABLE :: SO2_lifetime_mth(:, :, :)
76	REAL, ALLOCATABLE :: OCS_lifetime_tmp(:, :)
77	REAL, ALLOCATABLE :: SO2_lifetime_tmp(:, :)
78	REAL OCS_lifetime_glo(klon_glo,klev)
79	REAL SO2_lifetime_glo(klon_glo,klev)
80
81	REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :)
82	REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :)
83	REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :)
84
85	REAL H2SO4_lifetime_glo(klon_glo,klev)
86	!
87	REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
88	REAL, ALLOCATABLE, SAVE :: SO2_clim(:,:)
89	!$OMP THREADPRIVATE(OCS_clim,SO2_clim)
90	!
91	INTEGER i, k, kk, j
92	REAL p_bound
93
94	! For NetCDF:
95	INTEGER :: ncid_in ! IDs for input files
96	INTEGER :: varid, ncerr
97	CHARACTER (LEN=3) :: nc_lat, nc_lon
98	CHARACTER (LEN=256) :: nc_fname
99
100	INTEGER, PARAMETER :: lev_input=17
101	!--pressure at interfaces of input data (in Pa)
102	REAL, DIMENSION(lev_input+1), PARAMETER :: &
103	paprs_input=(/ &
104	1.00000002e+05, 6.06530673e+04, 3.67879449e+04, &
105	2.23130165e+04, 1.35335286e+04, 8.20850004e+03, &
106	4.97870695e+03, 3.01973841e+03, 1.83156393e+03, &
107	1.11089968e+03, 6.73794715e+02, 4.08677153e+02, &
108	2.47875223e+02, 1.50343923e+02, 9.11881985e+01, &
109	5.53084382e+01, 3.35462635e+01, 0.0 /)
110	!
111	IF (.NOT.ALLOCATED(OCS_clim)) ALLOCATE(OCS_clim(klon,klev))
112	IF (.NOT.ALLOCATED(SO2_clim)) ALLOCATE(SO2_clim(klon,klev))
113
114	IF (debutphy.OR.mth_cur/=mth_pre) THEN
115
116	!--preparation of global fields
117	CALL gather(paprs, paprs_glo)
118
119	IF (is_mpi_root.AND.is_omp_root) THEN
120
121	!--init ncdf variables
122	IF(flag_newclim_file) THEN
123	nc_fname = "ocs_so2_h2so4_annual_lmdz.nc"
124	nc_lat = "LAT"
125	nc_lon = "LON"
126	ELSE
127	! old file for retro compatibility
128	nc_fname = "ocs_so2_annual_lmdz.nc"
129	nc_lat = "lat"
130	nc_lon = "lon"
131	ENDIF
132
133	!--reading emission files
134	CALL nf95_open(nc_fname, nf90_nowrite, ncid_in)
135	IF(flag_verbose_strataer) print *, "Open file ", nc_fname
136
137	CALL nf95_inq_varid(ncid_in, "LEV", varid)
138	CALL nf95_gw_var(ncid_in, varid, lev)
139	n_lev = size(lev)
140
141	CALL nf95_inq_varid(ncid_in, nc_lat, varid)
142	CALL nf95_gw_var(ncid_in, varid, latitude)
143	n_lat = size(latitude)
144
145	CALL nf95_inq_varid(ncid_in, nc_lon, varid)
146	CALL nf95_gw_var(ncid_in, varid, longitude)
147	n_lon = size(longitude)
148
149	CALL nf95_inq_varid(ncid_in, "TIME", varid)
150	CALL nf95_gw_var(ncid_in, varid, time)
151	n_mth = size(time)
152
153	IF (.NOT.ALLOCATED(OCS_clim_in)) ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
154	IF (.NOT.ALLOCATED(SO2_clim_in)) ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
155	IF (.NOT.ALLOCATED(O3_clim_in)) ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth))
156
157	IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
158	IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
159	IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth))
160
161	CALL nf95_inq_varid(ncid_in, "OCS", varid)
162	ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
163	IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid
164
165	CALL nf95_inq_varid(ncid_in, "SO2", varid)
166	ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
167	IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid
168
169	CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid)
170	ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in)
171	IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid
172
173	CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid)
174	ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in)
175	IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid
176
177	IF(flag_newclim_file) THEN
178	CALL nf95_inq_varid(ncid_in, "O3", varid)
179	ncerr = nf90_get_var(ncid_in, varid, O3_clim_in)
180	IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid
181
182	CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid)
183	ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in)
184	IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid
185	ENDIF
186
187	CALL nf95_close(ncid_in)
188
189	IF (.NOT.ALLOCATED(OCS_clim_mth)) ALLOCATE(OCS_clim_mth(n_lon, n_lat, n_lev))
190	IF (.NOT.ALLOCATED(SO2_clim_mth)) ALLOCATE(SO2_clim_mth(n_lon, n_lat, n_lev))
191	IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
192	IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
193	IF (.NOT.ALLOCATED(O3_clim_mth)) ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev))
194	IF (.NOT.ALLOCATED(O3_clim_tmp)) ALLOCATE(O3_clim_tmp(klon_glo, n_lev))
195
196	IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
197	IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
198	IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
199	IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
200	IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev))
201	IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev))
202
203	!---select the correct month, undo multiplication with 1.e12 (precision reasons)
204	!---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
205	DO j=1,n_lat
206	SO2_clim_mth(:,j,:) = 1.e-12SO2_clim_in(:,n_lat+1-j,:,mth_cur)mSO2mol/mAIRmol
207	OCS_clim_mth(:,j,:) = 1.e-12OCS_clim_in(:,n_lat+1-j,:,mth_cur)mOCSmol/mAIRmol
208	SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
209	OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
210
211	! O3 from 2d model is not tracer, in VMR
212	IF(flag_newclim_file) THEN
213	H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur)
214	! new input files
215	O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur)
216	ELSE
217	H2SO4_lifetime_mth(:,j,:) = 1.e-6
218	O3_clim_mth(:,j,:) = 1.e-6
219	ENDIF
220	ENDDO
221
222	!---reduce to a klon_glo grid but keep the levels
223	CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
224	CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
225	CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp)
226
227	CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
228	CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
229	CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp)
230
231	!--set lifetime to very high value in uninsolated areas
232	DO i=1, klon_glo
233	DO kk=1, n_lev
234	!--added tests on VMR of species for realism
235	IF (OCS_clim_tmp(i,kk)<1.e-20) THEN
236	OCS_clim_tmp(i,kk)=1.0e-20
237	ENDIF
238	IF (SO2_clim_tmp(i,kk)<1.e-20) THEN
239	SO2_clim_tmp(i,kk)=1.0e-20
240	ENDIF
241	! 50 ppbv min for O3
242	IF (O3_clim_tmp(i,kk)<50.0e-9) THEN
243	O3_clim_tmp(i,kk)=50.0e-9
244	ENDIF
245
246	IF (OCS_lifetime_tmp(i,kk)==0.0) THEN
247	OCS_lifetime_tmp(i,kk)=1.0e12
248	ENDIF
249	IF (SO2_lifetime_tmp(i,kk)==0.0) THEN
250	SO2_lifetime_tmp(i,kk)=1.0e12
251	ENDIF
252	IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN
253	H2SO4_lifetime_tmp(i,kk)=1.0e12
254	ENDIF
255	ENDDO
256	ENDDO
257
258	!---regrid weighted lifetime and climatologies
259	DO i=1, klon_glo
260	DO k=1, klev
261	OCS_lifetime_glo(i,k)=0.0
262	SO2_lifetime_glo(i,k)=0.0
263	O3_clim_glo(i,k)=0.0
264
265	OCS_clim_glo(i,k)=0.0
266	SO2_clim_glo(i,k)=0.0
267	H2SO4_lifetime_glo(i,k)=0.0
268
269	DO kk=1, n_lev
270	OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
271	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
272	*OCS_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
273	SO2_lifetime_glo(i,k)=SO2_lifetime_glo(i,k)+ &
274	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
275	*SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
276	IF(flag_newclim_file) THEN
277	H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ &
278	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
279	-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
280	*H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
281	ENDIF
282
283	OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
284	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
285	*OCS_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
286	SO2_clim_glo(i,k)=SO2_clim_glo(i,k)+ &
287	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
288	*SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
289	IF(flag_newclim_file) THEN
290	O3_clim_glo(i,k)=O3_clim_glo(i,k)+ &
291	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
292	-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
293	*O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
294	ENDIF
295	ENDDO
296	ENDDO
297	ENDDO
298
299	ENDIF !--is_mpi_root and is_omp_root
300
301	!--keep memory of previous month
302	mth_pre=mth_cur
303
304	!--scatter global fields around
305	CALL scatter(OCS_clim_glo, OCS_clim)
306	CALL scatter(SO2_clim_glo, SO2_clim)
307	CALL scatter(O3_clim_glo, O3_clim)
308	CALL scatter(OCS_lifetime_glo, OCS_lifetime)
309	CALL scatter(SO2_lifetime_glo, SO2_lifetime)
310	CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime)
311
312	IF (is_mpi_root.AND.is_omp_root) THEN
313	!
314	DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in)
315	DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth)
316	DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp)
317	DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in)
318	DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth)
319	DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp)
320	!
321	ENDIF !--is_mpi_root and is_omp_root
322
323	ENDIF ! debutphy.OR.new month
324
325	!--set to background value everywhere in the very beginning, later only in the troposphere
326	!--a little dangerous as the MAXVAL is not computed on the global field
327	IF (debutphy.AND.MAXVAL(tr_seri)<1.e-30) THEN
328	p_bound=0.0
329	ELSE
330	p_bound=50000.
331	ENDIF
332
333	!--regridding tracer concentration on the vertical
334	DO i=1, klon
335	DO k=1, klev
336	!
337	!--OCS and SO2 prescribed back to their clim values below p_bound
338	IF (paprs(i,k)>p_bound) THEN
339	budg_3D_backgr_ocs(i,k)=OCS_clim(i,k)-tr_seri(i,k,id_OCS_strat)
340	budg_3D_backgr_so2(i,k)=SO2_clim(i,k)-tr_seri(i,k,id_SO2_strat)
341	tr_seri(i,k,id_OCS_strat)=OCS_clim(i,k)
342	tr_seri(i,k,id_SO2_strat)=SO2_clim(i,k)
343	ENDIF
344	ENDDO
345	ENDDO
346
347	!convert SO2_backgr_tend from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
348	DO i=1, klon
349	DO k=1, klev
350	budg_3D_backgr_ocs(i,k)=budg_3D_backgr_ocs(i,k)mSatom/mOCSmol(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
351	budg_3D_backgr_so2(i,k)=budg_3D_backgr_so2(i,k)mSatom/mSO2mol(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
352	ENDDO
353	ENDDO
354
355	RETURN
356
357	END SUBROUTINE interp_sulf_input

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

Download in other formats:

Original Format