Context Navigation

regr_horiz_time_climoz_m.F90 @ 3605

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

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

Download in other formats:

Original Format