Context Navigation

interp_sulf_input.F90 @ 3664

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

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

Download in other formats:

Original Format