Context Navigation

regr_horiz_time_climoz_m.F90 @ 5006

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

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

Download in other formats:

Original Format