Context Navigation

regr_horiz_time_climoz_m.f90 @ 5452

Last change on this file since 5452 was 5353, checked in by yann meurdesoif, 4 weeks ago
Nvidia compiler has some difficulties to compile correctly some complex array constructor. This commit decompose it the several phases in order to achieve the compilation. Please Lionel and David, have a look to this in order to validate. Probably, in future, when compiler heuristic will be improved, this commit can be reversed. YM
File size: 32.4 KB

Line
1	MODULE regr_horiz_time_climoz_m
2
3	USE interpolation, ONLY: locate
4	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
5	USE nrtype, ONLY: pi
6	USE netcdf, ONLY: NF90_CLOBBER, NF90_FLOAT, NF90_OPEN, &
7	NF90_NOWRITE, NF90_NOERR, NF90_GET_ATT, NF90_GLOBAL
8	USE netcdf95, ONLY: NF95_DEF_DIM, NF95_INQ_DIMID, NF95_INQUIRE_DIMENSION, &
9	NF95_DEF_VAR, NF95_INQ_VARID, NF95_INQUIRE_VARIABLE, &
10	NF95_OPEN, NF95_CREATE, NF95_GET_ATT, NF95_GW_VAR, nf95_get_var, &
11	NF95_CLOSE, NF95_ENDDEF, NF95_PUT_ATT, NF95_PUT_VAR, NF95_COPY_ATT
12	USE print_control_mod, ONLY: lunout
13	USE dimphy
14	IMPLICIT NONE
15	PRIVATE
16	PUBLIC :: regr_horiz_time_climoz
17	REAL, PARAMETER :: deg2rad=pi/180.
18	CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3 ','tro3_daylight']
19
20	INTEGER :: nlat_ou, nlon_ou
21	REAL, ALLOCATABLE :: latitude_glo(:)
22	!$OMP THREADPRIVATE(latitude_glo)
23	INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
24	!$OMP THREADPRIVATE(ind_cell_glo_glo)
25
26	CONTAINS
27
28	!-------------------------------------------------------------------------------
29	!
30	SUBROUTINE regr_horiz_time_climoz(read_climoz,interpt)
31	!
32	!-------------------------------------------------------------------------------
33	! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
34	! * Read ozone climatology from netcdf file
35	! * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
36	! * If interpt=T, interpolate linearly in time (one record each day)
37	! If interpt=F, keep original time sampling (14 months).
38	! * Save it to a new netcdf file.
39	!-------------------------------------------------------------------------------
40	! Remarks:
41	! * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
42	! * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
43	! * Output grid cells centers coordinates given by [rlonv,] rlatu.
44	! * Output grid cells edges coordinates given by [rlonu,] rlatv.
45	! * Input file [longitudes and] latitudes given in degrees.
46	! * Input file pressure levels are given in Pa or hPa.
47	! * All coordinates variables are stricly monotonic.
48	! * Monthly fields are interpolated linearly in time to get daily values.
49	! * Fields are known at the middle of the months, so interpolation requires an
50	! additional record both for 1st half of january and 2nd half of december:
51	! - For a 14-records "climoz.nc": records 1 and 14.
52	! - For 12-records files:
53	! record 12 of "climoz_m.nc" if available, or record 1 of "climoz.nc".
54	! record 1 of "climoz_p.nc" if available, or record 12 of "climoz.nc".
55	! * Calendar is taken into account to get one record each day (not 360 always).
56	! * Missing values are filled in from sky to ground by copying lowest valid one.
57	! Attribute "missing_value" or "_FillValue" must be present in input file.
58	!-------------------------------------------------------------------------------
59	USE assert_m, ONLY: assert
60	USE cal_tools_m, ONLY: year_len, mid_month
61	!! USE control_mod, ONLY: anneeref
62	USE time_phylmdz_mod, ONLY: annee_ref
63	USE ioipsl, ONLY: ioget_year_len, ioget_calendar
64	USE regr_conserv_m, ONLY: regr_conserv
65	USE regr_lint_m, ONLY: regr_lint
66	USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
67	USE slopes_m, ONLY: slopes
68	USE lmdz_xios
69	USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
70	USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
71	USE mod_grid_phy_lmdz, ONLY: klon_glo
72
73	!-------------------------------------------------------------------------------
74	! Arguments:
75	INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
76	! 1: read a single ozone climatology used day and night
77	! 2: same + read also a daylight climatology
78	LOGICAL, INTENT(IN) :: interpt ! TRUE => daily interpolation
79	! FALSE => no interpolation (14 months)
80	!-------------------------------------------------------------------------------
81	! Local variables:
82
83	!--- Input files variables
84	INTEGER :: nlon_in ! Number of longitudes
85	INTEGER :: nlat_in ! Number of latitudes
86	INTEGER :: nlev_in ! Number of pressure levels
87	INTEGER :: nmth_in ! Number of months
88	REAL, ALLOCATABLE:: lon_in(:) ! Longitudes (ascending order, rad)
89	REAL, ALLOCATABLE:: lat_in(:) ! Latitudes (ascending order, rad)
90	REAL, ALLOCATABLE:: lev_in(:) ! Pressure levels (ascen. order, hPa)
91	REAL, ALLOCATABLE :: lon_in_edge(:) ! Longitude intervals edges
92	! (ascending order, / )
93	REAL, ALLOCATABLE :: sinlat_in_edge(:) ! Sinus of latitude intervals edges
94	! (ascending order, / )
95	LOGICAL :: ldec_lon, ldec_lat, ldec_lev ! Decreasing order in input file
96	CHARACTER(LEN=20) :: cal_in ! Calendar
97	REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:) ! Ozone climatologies
98	REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:) ! Ozone climatologies
99	REAL, ALLOCATABLE :: o3_in2 (:,:,:,:) ! Ozone climatologies
100	REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:) ! Ozone climatologies
101	! last index: 1 for the day-night average, 2 for the daylight field.
102	REAL :: NaN
103
104	!--- Partially or totally regridded variables (:,:,nlev_in,:,read_climoz)
105	REAL, ALLOCATABLE :: o3_regr_lon (:,:,:,:,:) ! (nlon_ou,nlat_in,:,0:13 ,:)
106	REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13 ,:)
107	REAL, ALLOCATABLE :: o3_out3 (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
108	REAL, ALLOCATABLE :: o3_out3_glo (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
109	REAL, ALLOCATABLE :: o3_regr_lat (:,:,:,:) ! (nlat_in,:,0:13 ,:)
110	REAL, ALLOCATABLE :: o3_out2 (:,:,:,:) ! (nlat_ou,:,ntim_ou,:)
111	REAL, ALLOCATABLE :: o3_out2_glo (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
112	REAL, ALLOCATABLE :: o3_out (:,:,:,:) ! (nbp_lat,:,ntim_ou,:)
113	! Dimension number \| Interval \| Contains \| For variables:
114	! 1 (longitude) \| [rlonu(i-1), rlonu(i)] \| rlonv(i) \| all
115	! 2 (latitude) \| [rlatv(j), rlatv(j-1)] \| rlatu(j) \| all but o3_regr_lon
116	! 3 (press level) \| \| lev(k) \| all
117	! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
118	! Dimension 4 is: month number (all vars but o3_out)
119	! days elapsed since Jan. 1st 0h at mid-day (o3_out only)
120	REAL, ALLOCATABLE :: v1(:)
121
122	!--- For NetCDF:
123	INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
124	INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
125	INTEGER, ALLOCATABLE :: dIDs(:)
126	CHARACTER(LEN=20) :: cal_ou !--- Calendar; no time inter => same as input
127	CHARACTER(LEN=80) :: press_unit !--- Pressure unit
128	REAL :: tmidmonth(0:13) !--- Elapsed days since Jan-1 0h at mid-months
129	! Additional records 0, 13 for interpolation
130	REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
131	LOGICAL :: lprev, lnext !--- Flags: previous/next files are present
132	LOGICAL :: l3D, l2D !--- Flag: input fields are 3D or zonal
133	INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
134	INTEGER, ALLOCATABLE :: sta(:), cnt(:)
135	CHARACTER(LEN=80) :: sub, dim_nam, msg
136	REAL :: null_array(0)
137	LOGICAL,SAVE :: first=.TRUE.
138	!$OMP THREADPRIVATE(first)
139	REAL, ALLOCATABLE :: test_o3_in(:,:)
140	REAL, ALLOCATABLE :: test_o3_out(:)
141	REAL,ALLOCATABLE :: tmp(:)
142
143
144	IF (grid_type==unstructured) THEN
145	IF (first) THEN
146	IF (is_master) THEN
147	ALLOCATE(latitude_glo(klon_glo))
148	ALLOCATE(ind_cell_glo_glo(klon_glo))
149	ELSE
150	ALLOCATE(latitude_glo(0))
151	ALLOCATE(ind_cell_glo_glo(0))
152	ENDIF
153	CALL gather(latitude_deg, latitude_glo)
154	CALL gather(ind_cell_glo, ind_cell_glo_glo)
155	ENDIF
156	ENDIF
157
158	IF (is_omp_master) THEN
159	nlat_ou=nbp_lat
160	nlon_ou=nbp_lon
161
162	!-------------------------------------------------------------------------------
163	IF (is_mpi_root) THEN
164	sub="regr_horiz_time_climoz"
165	WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
166	CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
167
168	CALL NF95_OPEN("climoz.nc" , NF90_NOWRITE, fID_in)
169	lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
170	lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
171
172	!--- Get coordinates from the input file. Converts lon/lat in radians.
173	! Few inversions because "regr_conserv" and gcm need ascending vectors.
174	CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
175	CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
176	l3D=ndims==4; l2D=ndims==3
177	IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
178	IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
179	DO i=1,ndims
180	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
181	CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
182	ii=i; IF(l2D) ii=i+1 !--- ndims==3:NO LONGITUDE
183	SELECT CASE(ii)
184	CASE(1) !--- LONGITUDE
185	CALL NF95_GW_VAR(fID_in, varid, lon_in)
186	ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
187	nlon_in=dln; lon_in=lon_in*deg2rad
188	CASE(2) !--- LATITUDE
189	CALL NF95_GW_VAR(fID_in, varid, lat_in)
190	ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
191	nlat_in=dln; lat_in=lat_in*deg2rad
192	CASE(3) !--- PRESSURE LEVELS
193	CALL NF95_GW_VAR(fID_in, varid, lev_in)
194	ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
195	nlev_in=dln
196	CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
197	k=LEN_TRIM(press_unit)
198	DO WHILE(ICHAR(press_unit(k:k))==0)
199	press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
200	END DO
201	IF(press_unit == "Pa") THEN
202	lev_in = lev_in/100. !--- CONVERT TO hPa
203	ELSE IF(press_unit /= "hPa") THEN
204	CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
205	END IF
206	CASE(4) !--- TIME
207	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
208	cal_in='gregorian'
209	IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR) &
210	WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'// &
211	TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
212	k=LEN_TRIM(cal_in)
213	DO WHILE(ICHAR(cal_in(k:k))==0)
214	cal_in(k:k)=' '; k=LEN_TRIM(cal_in) !--- REMOVE NULL END CHAR
215	END DO
216	END SELECT
217	END DO
218
219	!--- Prepare quantities for time interpolation
220	tmidmonth=mid_month(annee_ref, cal_in)
221	IF(interpt) THEN
222	ntim_ou=ioget_year_len(annee_ref)
223	ALLOCATE(tmidday(ntim_ou))
224	tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
225	CALL ioget_calendar(cal_ou)
226	ELSE
227	ntim_ou=14
228	cal_ou=cal_in
229	END IF
230	ENDIF
231
232	IF (grid_type==unstructured) THEN
233	CALL bcast_mpi(nlon_in)
234	CALL bcast_mpi(nlat_in)
235	CALL bcast_mpi(nlev_in)
236	CALL bcast_mpi(l3d)
237	CALL bcast_mpi(tmidmonth)
238	IF(interpt) CALL bcast_mpi(tmidday)
239	CALL bcast_mpi(ntim_ou)
240
241	IF (is_mpi_root) THEN
242	CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
243	IF (l3D) THEN
244	CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
245	ELSE
246	CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
247	ENDIF
248	ELSE
249	CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
250	IF (l3D) THEN
251	CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
252	ELSE
253	CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
254	ENDIF
255	ENDIF
256	CALL xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
257	CALL xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
258	CALL xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
259	CALL xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
260	CALL xios_set_field_attr("tro3_out", enabled=.TRUE.)
261	CALL xios_set_field_attr("tro3_out", enabled=.TRUE.)
262
263	IF (first) THEN
264	first=.FALSE.
265	RETURN
266	ENDIF
267	ENDIF
268
269
270	IF (is_mpi_root) THEN
271	!--- Longitudes management:
272	! * Need to shift data if the origin of input file longitudes /= -pi
273	! * Need to add some margin in longitude to ensure input interval contains
274	! all the output intervals => at least one longitudes slice has to be
275	! duplicated, possibly more for undersampling.
276	IF(l3D) THEN
277	IF (grid_type==unstructured) THEN
278	dx2=0
279	ELSE
280	!--- Compute input edges longitudes vector (no end point yet)
281	ALLOCATE(v1(nlon_in+1))
282	v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
283	FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
284	v1(nlon_in+1)=v1(1)+2.*pi
285	DEALLOCATE(lon_in)
286
287	!--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
288	v1=v1+2piREAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
289
290	!--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
291	dx1=locate(v1,boundslon_reg(1,west))-1
292	v1=CSHIFT(v1,SHIFT=dx1,DIM=1)
293	v1(nlon_in-dx1+2:)=v1(nlon_in-dx1+2:)+2.*pi
294
295	!--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
296	dx2=0; DO WHILE(v1(1+dx2)+2.*pi<=boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
297
298	!--- Final edges longitudes vector (with margin and end point)
299	ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
300	DEALLOCATE(v1)
301	ENDIF
302	END IF
303
304	!--- Compute sinus of intervals edges latitudes:
305	ALLOCATE(sinlat_in_edge(nlat_in+1))
306	sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
307	FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
308	DEALLOCATE(lat_in)
309
310
311
312	!--- Check for contiguous years:
313	ib=0; ie=13
314	IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
315	WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
316	ELSE
317	IF( lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
318	IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
319	IF( lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
320	IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
321	IF(.NOT.lprev) ib=1
322	IF(.NOT.lnext) ie=12
323	END IF
324	ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
325	IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
326	IF(l2D) cnt=[ nlat_in,nlev_in,1]
327	IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
328	IF(l2D) ALLOCATE(o3_in2( nlat_in,nlev_in,ib:ie,read_climoz))
329
330	!--- Read full current file and one record each available contiguous file
331	DO iv=1,read_climoz
332	CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
333	IF(l3D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
334	IF(l2D) call NF95_GET_VAR(fID_in, vID_in(iv), o3_in2( :,:,1:12,iv))
335	IF(lprev) THEN; sta(ndims)=12
336	CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
337	IF(l3D) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
338	IF(l2d) call NF95_GET_VAR(fID_in_m,vID_in(iv),o3_in2( :,:, 0,iv),sta,cnt)
339	END IF
340	IF(lnext) THEN; sta(ndims)=1
341	CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
342	IF(l3D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
343	IF(l2D) call NF95_GET_VAR(fID_in_p,vID_in(iv),o3_in2( :,:,13,iv),sta,cnt)
344	END IF
345	END DO
346	IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
347	IF(lprev) CALL NF95_CLOSE(fID_in_m)
348	IF(lnext) CALL NF95_CLOSE(fID_in_p)
349
350	!--- Revert decreasing coordinates vector
351	IF(l3D) THEN
352	IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
353	IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
354	IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
355
356	IF (grid_type /= unstructured) THEN
357	!--- Shift values for longitude and duplicate some longitudes slices
358	o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
359	o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
360	ENDIF
361	ELSE
362	IF(ldec_lat) o3_in2 = o3_in2( nlat_in:1:-1,:,:,:)
363	IF(ldec_lev) o3_in2 = o3_in2( :,nlev_in:1:-1,:,:)
364	END IF
365
366	!--- Deal with missing values
367	DO m=1, read_climoz
368	WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
369	IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
370	IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
371	WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
372	END IF
373	END IF
374	WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
375	WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
376
377	!--- Check top layer contains no NaNs & search NaNs from top to ground
378	msg=TRIM(sub)//": NaNs in top layer !"
379	IF(l3D) THEN
380	IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
381	DO k = 2,nlev_in
382	WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
383	END DO
384	ELSE
385	IF(ANY(o3_in2( :,1,:,m)==NaN)) THEN
386	WRITE(lunout,*)msg
387	!--- Fill in latitudes where all values are missing
388	DO l=1,nmth_in
389	!--- Next to south pole
390	j=1; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
391	IF(j>1) &
392	o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
393	!--- Next to north pole
394	j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
395	IF(j<nlat_in) &
396	o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
397	END DO
398	END IF
399
400	!--- Fill in high latitudes missing values
401	!--- Highest level been filled-in, so has always valid values.
402	DO k = 2,nlev_in
403	WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
404	END DO
405	END IF
406	END DO
407
408	ENDIF
409
410	!=============================================================================
411	IF(l3D) THEN !=== 3D FIELDS
412	!=============================================================================
413	IF (grid_type==unstructured) THEN
414	nlat_ou=klon_mpi
415
416	IF (is_mpi_root) THEN
417	ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz))
418	o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
419	ELSE
420	ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
421	ENDIF
422	ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
423
424	CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
425	CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
426	ELSE
427
428	!--- Regrid in longitude
429	ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
430	tmp = [boundslon_reg(1,west),boundslon_reg(:,east)]
431	CALL regr_conserv(1, o3_in3, xs = lon_in_edge, &
432	xt = tmp , &
433	vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
434	DEALLOCATE(o3_in3)
435
436	!--- Regrid in latitude: averaging with respect to SIN(lat) is
437	! equivalent to weighting by COS(lat)
438	!--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
439	ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
440	tmp = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.]
441	CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge, &
442	xt = tmp, &
443	vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:), &
444	slope = slopes(2,o3_regr_lon, sinlat_in_edge))
445	DEALLOCATE(o3_regr_lon)
446
447	ENDIF
448
449	!--- Duplicate previous/next record(s) if they are not available
450	IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
451	IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
452
453	!--- Regrid in time by linear interpolation:
454	ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
455	IF( interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
456	IF(.NOT.interpt) o3_out3=o3_regr_lonlat
457	DEALLOCATE(o3_regr_lonlat)
458
459	nlat_ou=nbp_lat
460	IF (grid_type==unstructured) THEN
461	CALL xios_send_field('o3_out',o3_out3)
462	ndims=3
463	ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
464	CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
465	ENDIF
466
467	!--- Create the output file and get the variable IDs:
468	CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
469	ndims, cal_ou)
470
471	IF (is_mpi_root) THEN
472	!--- Write remaining coordinate variables:
473	CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
474	IF( interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
475	IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
476
477	!--- Write to file (the order of "rlatu" is inverted in the output file):
478	IF (grid_type==unstructured) THEN
479
480	ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
481	DO i=1,klon_glo
482	o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
483	ENDDO
484
485	DO m = 1, read_climoz
486	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
487	END DO
488
489	ELSE
490	DO m = 1, read_climoz
491	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
492	END DO
493	ENDIF
494	CALL NF95_CLOSE(fID_ou)
495
496
497	ENDIF
498
499
500	!=============================================================================
501	ELSE !=== ZONAL FIELDS
502	!=============================================================================
503
504	IF (grid_type==unstructured) THEN
505	nlat_ou=klon_mpi
506
507	IF (is_mpi_root) THEN
508	ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz))
509	o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
510	ELSE
511	ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
512	ENDIF
513	ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
514	CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
515	CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
516	IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
517	IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
518
519	ELSE
520	!--- Regrid in latitude: averaging with respect to SIN(lat) is
521	! equivalent to weighting by COS(lat)
522	!--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
523	ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
524	tmp = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.]
525	CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge, &
526	xt = tmp, &
527	vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:), &
528	slope = slopes(1,o3_in2, sinlat_in_edge))
529	DEALLOCATE(o3_in2)
530
531	!--- Duplicate previous/next record(s) if they are not available
532	IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
533	IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
534
535	ENDIF
536
537	!--- Regrid in time by linear interpolation:
538	ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
539	IF( interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
540	IF(.NOT.interpt) o3_out2=o3_regr_lat
541	DEALLOCATE(o3_regr_lat)
542
543	nlat_ou=nbp_lat
544
545	IF (grid_type==unstructured) THEN
546	ndims=3
547	ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
548	CALL gather_mpi(o3_out2, o3_out2_glo)
549	ENDIF
550
551	!--- Create the output file and get the variable IDs:
552	CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
553	ndims, cal_ou)
554
555	IF (is_mpi_root) THEN
556
557	!--- Write remaining coordinate variables:
558	CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
559	IF( interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
560	IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
561
562	IF (grid_type==unstructured) THEN
563
564	ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
565	DO i=1,klon_glo
566	o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
567	ENDDO
568
569
570	DO m = 1, read_climoz
571	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
572	END DO
573	ELSE
574	!--- Write to file (the order of "rlatu" is inverted in the output file):
575	DO m = 1, read_climoz
576	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
577	END DO
578	ENDIF
579
580	CALL NF95_CLOSE(fID_ou)
581
582	ENDIF
583
584	!=============================================================================
585	END IF
586	!=============================================================================
587
588	IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
589
590	ENDIF ! is_omp_master
591
592	first=.FALSE.
593	END SUBROUTINE regr_horiz_time_climoz
594	!
595	!-------------------------------------------------------------------------------
596
597
598	!-------------------------------------------------------------------------------
599	!
600	SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
601	vID_ou, ndims, cal_ou)
602	!-------------------------------------------------------------------------------
603	! Purpose: This subroutine creates the NetCDF output file, defines
604	! dimensions and variables, and writes some of the coordinate variables.
605	!-------------------------------------------------------------------------------
606	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
607	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
608	USE mod_phys_lmdz_para, ONLY: is_mpi_root
609	USE mod_grid_phy_lmdz, ONLY: klon_glo
610	!
611	!-------------------------------------------------------------------------------
612	! Arguments:
613	INTEGER, INTENT(IN) :: fID_in, nlev_in, ntim_ou
614	INTEGER, INTENT(OUT) :: fID_ou, vlevID, vtimID
615	INTEGER, INTENT(OUT) :: vID_ou(:) ! dim(1/2) 1: O3day&night 2: O3daylight
616	INTEGER, INTENT(IN) :: ndims ! fields rank (3 or 4)
617	CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
618	!-------------------------------------------------------------------------------
619	! Local variables:
620	INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
621	INTEGER :: vlonID, vlatID, ncerr, is
622	REAL,ALLOCATABLE :: latitude_glo_(:)
623	CHARACTER(LEN=80) :: sub
624	INTEGER :: i
625
626
627	!-------------------------------------------------------------------------------
628
629	IF (is_mpi_root) THEN
630	sub="prepare_out"
631	WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
632	CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
633
634	!--- Dimensions:
635	IF(ndims==4) &
636	CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
637	CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
638	CALL NF95_DEF_DIM(fID_ou, "plev", nlev_in, dlevID)
639	CALL NF95_DEF_DIM(fID_ou, "time", ntim_ou, dtimID)
640
641	!--- Define coordinate variables:
642	IF(ndims==4) &
643	CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
644	CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
645	CALL NF95_DEF_VAR(fID_ou, "plev", NF90_FLOAT, dlevID, vlevID)
646	CALL NF95_DEF_VAR(fID_ou, "time", NF90_FLOAT, dtimID, vtimID)
647	IF(ndims==4) &
648	CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
649	CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
650	CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
651	CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
652	IF(ndims==4) &
653	CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
654	CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
655	CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
656	CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
657	CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name", "air pressure")
658	CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar", cal_ou)
659
660	!--- Define the main variables:
661	IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
662	IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
663	CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
664	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
665	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
666	&_in_air")
667	IF(SIZE(vID_ou) == 2) THEN
668	CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
669	CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
670	&ylight")
671	END IF
672
673	!--- Global attributes:
674	! The following commands, copying attributes, may fail. That is OK.
675	! It should just mean that the attribute is not defined in the input file.
676	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
677	CALL handle_err_copy_att("Conventions")
678	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title", fID_ou,NF90_GLOBAL, ncerr)
679	CALL handle_err_copy_att("title")
680	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
681	CALL handle_err_copy_att("institution")
682	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source", fID_ou,NF90_GLOBAL, ncerr)
683	CALL handle_err_copy_att("source")
684	CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
685	CALL NF95_ENDDEF(fID_ou)
686
687	IF (grid_type==unstructured) THEN
688	ALLOCATE(latitude_glo_(klon_glo))
689	DO i=1,klon_glo
690	latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
691	ENDDO
692	CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
693	ELSE
694	!--- Write one of the coordinate variables:
695	IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
696	CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
697	! (convert from rad to degrees and sort in ascending order)
698	ENDIF
699	ENDIF
700
701	CONTAINS
702
703	!-------------------------------------------------------------------------------
704	!
705	SUBROUTINE handle_err_copy_att(att_name)
706	!
707	!-------------------------------------------------------------------------------
708	USE netcdf, ONLY: NF90_NOERR, NF90_strerror
709	!-------------------------------------------------------------------------------
710	! Arguments:
711	CHARACTER(LEN=*), INTENT(IN) :: att_name
712	!-------------------------------------------------------------------------------
713	IF(ncerr /= NF90_NOERR) &
714	WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)// &
715	" -- "//TRIM(NF90_strerror(ncerr))
716
717	END SUBROUTINE handle_err_copy_att
718	!
719	!-------------------------------------------------------------------------------
720
721	END SUBROUTINE prepare_out
722	!
723	!-------------------------------------------------------------------------------
724
725	END MODULE regr_horiz_time_climoz_m
726	!
727	!-------------------------------------------------------------------------------

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

Download in other formats:

Original Format