Context Navigation

interp_sulf_input.F90 @ 3158

Last change on this file since 3158 was 3097, checked in by oboucher, 7 years ago
minor change to diag to avoid copy of arrays
File size: 9.4 KB

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

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

Download in other formats:

Original Format