Context Navigation

interp_sulf_input.F90 @ 3526

Last change on this file since 3526 was 3526, checked in by Laurent Fairhead, 5 years ago
Olivier's modifications + addition of Id keyword
Property svn:keywords set to `Id`
File size: 9.4 KB

Line
1	!
2	! $Id: interp_sulf_input.F90 3526 2019-05-28 13:00:44Z fairhead $
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 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
15	USE mod_phys_lmdz_para
16	USE dimphy
17	USE phys_cal_mod
18	USE infotrac
19	USE aerophys
20	USE YOMCST
21
22	IMPLICIT NONE
23
24	include "dimensions.h"
25
26	! Variable input
27	REAL paprs(klon,klev+1)
28	REAL tr_seri(klon,klev,nbtr)
29	REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
30	LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
31
32	! Variables locales
33	INTEGER n_lat ! number of latitudes in the input data
34	INTEGER n_lon ! number of longitudes in the input data
35	INTEGER, SAVE :: n_lev ! number of levels in the input data
36	INTEGER n_mth ! number of months in the input data
37	INTEGER, SAVE :: mth_pre
38	!$OMP THREADPRIVATE(mth_pre)
39
40	! Champs reconstitues
41	REAL paprs_glo(klon_glo,klev+1)
42
43	REAL, POINTER:: latitude(:)
44	! (of input data sorted in strictly ascending order)
45
46	REAL, POINTER:: longitude(:)
47	! (of input data sorted in strictly ascending order)
48
49	REAL, POINTER:: time(:)
50	! (of input data sorted in strictly ascending order)
51
52	REAL, POINTER:: lev(:)
53	! levels of input data
54
55	REAL, ALLOCATABLE :: OCS_clim_in(:, :, :, :)
56	REAL, ALLOCATABLE :: SO2_clim_in(:, :, :, :)
57	REAL, ALLOCATABLE :: OCS_clim_mth(:, :, :)
58	REAL, ALLOCATABLE :: SO2_clim_mth(:, :, :)
59	REAL, ALLOCATABLE :: OCS_clim_tmp(:, :)
60	REAL, ALLOCATABLE :: SO2_clim_tmp(:, :)
61	REAL OCS_clim_glo(klon_glo,klev)
62	REAL SO2_clim_glo(klon_glo,klev)
63	REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
64	REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
65	REAL, ALLOCATABLE :: OCS_lifetime_mth(:, :, :)
66	REAL, ALLOCATABLE :: SO2_lifetime_mth(:, :, :)
67	REAL, ALLOCATABLE :: OCS_lifetime_tmp(:, :)
68	REAL, ALLOCATABLE :: SO2_lifetime_tmp(:, :)
69	REAL OCS_lifetime_glo(klon_glo,klev)
70	REAL SO2_lifetime_glo(klon_glo,klev)
71	!
72	REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
73	REAL, ALLOCATABLE, SAVE :: SO2_clim(:,:)
74	!$OMP THREADPRIVATE(OCS_clim,SO2_clim)
75	!
76	INTEGER i, k, kk, j
77	REAL p_bound
78
79	! For NetCDF:
80	INTEGER ncid_in ! IDs for input files
81	INTEGER varid, ncerr
82
83	INTEGER, PARAMETER :: lev_input=17
84	!--pressure at interfaces of input data (in Pa)
85	REAL, DIMENSION(lev_input+1), PARAMETER :: &
86	paprs_input=(/ &
87	1.00000002e+05, 6.06530673e+04, 3.67879449e+04, &
88	2.23130165e+04, 1.35335286e+04, 8.20850004e+03, &
89	4.97870695e+03, 3.01973841e+03, 1.83156393e+03, &
90	1.11089968e+03, 6.73794715e+02, 4.08677153e+02, &
91	2.47875223e+02, 1.50343923e+02, 9.11881985e+01, &
92	5.53084382e+01, 3.35462635e+01, 0.0 /)
93	!
94	IF (.NOT.ALLOCATED(OCS_clim)) ALLOCATE(OCS_clim(klon,klev))
95	IF (.NOT.ALLOCATED(SO2_clim)) ALLOCATE(SO2_clim(klon,klev))
96
97	IF (debutphy.OR.mth_cur.NE.mth_pre) THEN
98
99	!--preparation of global fields
100	CALL gather(paprs, paprs_glo)
101
102	IF (is_mpi_root.AND.is_omp_root) THEN
103
104	!--reading emission files
105	CALL nf95_open("ocs_so2_annual_lmdz.nc", nf90_nowrite, ncid_in)
106
107	CALL nf95_inq_varid(ncid_in, "LEV", varid)
108	CALL nf95_gw_var(ncid_in, varid, lev)
109	n_lev = size(lev)
110
111	CALL nf95_inq_varid(ncid_in, "lat", varid)
112	CALL nf95_gw_var(ncid_in, varid, latitude)
113	n_lat = size(latitude)
114
115	CALL nf95_inq_varid(ncid_in, "lon", varid)
116	CALL nf95_gw_var(ncid_in, varid, longitude)
117	n_lon = size(longitude)
118
119	CALL nf95_inq_varid(ncid_in, "TIME", varid)
120	CALL nf95_gw_var(ncid_in, varid, time)
121	n_mth = size(time)
122
123	IF (.NOT.ALLOCATED(OCS_clim_in)) ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
124	IF (.NOT.ALLOCATED(SO2_clim_in)) ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
125	IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
126	IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
127
128	CALL nf95_inq_varid(ncid_in, "OCS", varid)
129	ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
130	print *,'code erreur OCS=', ncerr, varid
131
132	CALL nf95_inq_varid(ncid_in, "SO2", varid)
133	ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
134	print *,'code erreur SO2=', ncerr, varid
135
136	CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid)
137	ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in)
138	print *,'code erreur OCS_lifetime_in=', ncerr, varid
139
140	CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid)
141	ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in)
142	print *,'code erreur SO2_lifetime_in=', ncerr, varid
143
144	CALL nf95_close(ncid_in)
145
146	IF (.NOT.ALLOCATED(OCS_clim_mth)) ALLOCATE(OCS_clim_mth(n_lon, n_lat, n_lev))
147	IF (.NOT.ALLOCATED(SO2_clim_mth)) ALLOCATE(SO2_clim_mth(n_lon, n_lat, n_lev))
148	IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
149	IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
150	IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
151	IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
152	IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
153	IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
154
155	!---select the correct month, undo multiplication with 1.e12 (precision reasons)
156	!---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
157	DO j=1,n_lat
158	SO2_clim_mth(:,j,:) = 1.e-12SO2_clim_in(:,n_lat+1-j,:,mth_cur)mSO2mol/mAIRmol
159	OCS_clim_mth(:,j,:) = 1.e-12OCS_clim_in(:,n_lat+1-j,:,mth_cur)mOCSmol/mAIRmol
160	SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
161	OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
162	ENDDO
163
164	!---reduce to a klon_glo grid but keep the levels
165	CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
166	CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
167	CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
168	CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
169
170	!--set lifetime to very high value in uninsolated areas
171	DO i=1, klon_glo
172	DO kk=1, n_lev
173	IF (OCS_lifetime_tmp(i,kk)==0.0) THEN
174	OCS_lifetime_tmp(i,kk)=1.0e12
175	ENDIF
176	IF (SO2_lifetime_tmp(i,kk)==0.0) THEN
177	SO2_lifetime_tmp(i,kk)=1.0e12
178	ENDIF
179	ENDDO
180	ENDDO
181
182	!---regrid weighted lifetime and climatologies
183	DO i=1, klon_glo
184	DO k=1, klev
185	OCS_lifetime_glo(i,k)=0.0
186	SO2_lifetime_glo(i,k)=0.0
187	OCS_clim_glo(i,k)=0.0
188	SO2_clim_glo(i,k)=0.0
189	DO kk=1, n_lev
190	OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
191	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
192	*OCS_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
193	SO2_lifetime_glo(i,k)=SO2_lifetime_glo(i,k)+ &
194	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
195	*SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
196	OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
197	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
198	*OCS_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
199	SO2_clim_glo(i,k)=SO2_clim_glo(i,k)+ &
200	MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
201	*SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
202	ENDDO
203	ENDDO
204	ENDDO
205
206	ENDIF !--is_mpi_root and is_omp_root
207
208	!--keep memory of previous month
209	mth_pre=mth_cur
210
211	!--scatter global fields around
212	CALL scatter(OCS_clim_glo, OCS_clim)
213	CALL scatter(SO2_clim_glo, SO2_clim)
214	CALL scatter(OCS_lifetime_glo, OCS_lifetime)
215	CALL scatter(SO2_lifetime_glo, SO2_lifetime)
216
217	IF (is_mpi_root.AND.is_omp_root) THEN
218	!
219	DEALLOCATE(OCS_clim_in,SO2_clim_in)
220	DEALLOCATE(OCS_clim_mth,SO2_clim_mth)
221	DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp)
222	DEALLOCATE(OCS_lifetime_in,SO2_lifetime_in)
223	DEALLOCATE(OCS_lifetime_mth,SO2_lifetime_mth)
224	DEALLOCATE(OCS_lifetime_tmp,SO2_lifetime_tmp)
225	!
226	ENDIF !--is_mpi_root and is_omp_root
227
228	ENDIF ! debutphy.OR.new month
229
230	!--set to background value everywhere in the very beginning, later only in the troposphere
231	!--a little dangerous as the MAXVAL is not computed on the global field
232	IF (debutphy.AND.MAXVAL(tr_seri).LT.1.e-30) THEN
233	p_bound=0.0
234	ELSE
235	p_bound=50000.
236	ENDIF
237
238	!--regridding tracer concentration on the vertical
239	DO i=1, klon
240	DO k=1, klev
241	!
242	!--OCS and SO2 prescribed back to their clim values below p_bound
243	IF (paprs(i,k).GT.p_bound) THEN
244	budg_3D_backgr_ocs(i,k)=OCS_clim(i,k)-tr_seri(i,k,id_OCS_strat)
245	budg_3D_backgr_so2(i,k)=SO2_clim(i,k)-tr_seri(i,k,id_SO2_strat)
246	tr_seri(i,k,id_OCS_strat)=OCS_clim(i,k)
247	tr_seri(i,k,id_SO2_strat)=SO2_clim(i,k)
248	ENDIF
249	ENDDO
250	ENDDO
251
252	!convert SO2_backgr_tend from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
253	DO i=1, klon
254	DO k=1, klev
255	budg_3D_backgr_ocs(i,k)=budg_3D_backgr_ocs(i,k)mSatom/mOCSmol(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
256	budg_3D_backgr_so2(i,k)=budg_3D_backgr_so2(i,k)mSatom/mSO2mol(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
257	ENDDO
258	ENDDO
259
260	RETURN
261
262	END SUBROUTINE interp_sulf_input

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

Download in other formats:

Original Format