Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

regr_horiz_time_climoz_m.F90 @ 2839

Last change on this file since 2839 was 2839, checked in by Laurent Fairhead, 7 years ago
Merged trunk changes r2785:2838 into testing branch
File size: 24.7 KB

Line
1	MODULE regr_horiz_time_climoz_m
2
3	USE interpolation, ONLY: locate
4	USE mod_grid_phy_lmdz, ONLY: nlon_ou => nbp_lon, nlat_ou => nbp_lat
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	IMPLICIT NONE
14	PRIVATE
15	PUBLIC :: regr_horiz_time_climoz
16	REAL, PARAMETER :: deg2rad=pi/180.
17	CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3 ','tro3_daylight']
18
19	CONTAINS
20
21	!-------------------------------------------------------------------------------
22	!
23	SUBROUTINE regr_horiz_time_climoz(read_climoz,interpt)
24	!
25	!-------------------------------------------------------------------------------
26	! Purpose: Regrid horizontally and in time zonal or 3D ozone climatologies.
27	! * Read ozone climatology from netcdf file
28	! * Regrid it horizontaly to LMDZ grid (quasi-conservative method)
29	! * If interpt=T, interpolate linearly in time (one record each day)
30	! If interpt=F, keep original time sampling (14 months).
31	! * Save it to a new netcdf file.
32	!-------------------------------------------------------------------------------
33	! Remarks:
34	! * Up to 2 variables treated: "tro3" and "tro3_daylight" (if read_climoz=2)
35	! * Input fields coordinates: (longitudes, latitudes, pressure_levels, time)
36	! * Output grid cells centers coordinates given by [rlonv,] rlatu.
37	! * Output grid cells edges coordinates given by [rlonu,] rlatv.
38	! * Input file [longitudes and] latitudes given in degrees.
39	! * Input file pressure levels are given in Pa or hPa.
40	! * All coordinates variables are stricly monotonic.
41	! * Monthly fields are interpolated linearly in time to get daily values.
42	! * Fields are known at the middle of the months, so interpolation requires an
43	! additional record both for 1st half of january and 2nd half of december:
44	! - For a 14-records "climoz.nc": records 1 and 14.
45	! - For 12-records files:
46	! record 12 of "climoz_m.nc" if available, or record 1 of "climoz.nc".
47	! record 1 of "climoz_p.nc" if available, or record 12 of "climoz.nc".
48	! * Calendar is taken into account to get one record each day (not 360 always).
49	! * Missing values are filled in from sky to ground by copying lowest valid one.
50	! Attribute "missing_value" or "_FillValue" must be present in input file.
51	!-------------------------------------------------------------------------------
52	USE assert_m, ONLY: assert
53	USE cal_tools_m, ONLY: year_len, mid_month
54	USE control_mod, ONLY: anneeref
55	USE ioipsl, ONLY: ioget_year_len, ioget_calendar
56	USE regr_conserv_m, ONLY: regr_conserv
57	USE regr_lint_m, ONLY: regr_lint
58	USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
59	USE slopes_m, ONLY: slopes
60	!-------------------------------------------------------------------------------
61	! Arguments:
62	INTEGER, INTENT(IN) :: read_climoz ! read ozone climatology, 1 or 2
63	! 1: read a single ozone climatology used day and night
64	! 2: same + read also a daylight climatology
65	LOGICAL, INTENT(IN) :: interpt ! TRUE => daily interpolation
66	! FALSE => no interpolation (14 months)
67	!-------------------------------------------------------------------------------
68	! Local variables:
69
70	!--- Input files variables
71	INTEGER :: nlon_in ! Number of longitudes
72	INTEGER :: nlat_in ! Number of latitudes
73	INTEGER :: nlev_in ! Number of pressure levels
74	INTEGER :: nmth_in ! Number of months
75	REAL, POINTER :: lon_in(:) ! Longitudes (ascending order, rad)
76	REAL, POINTER :: lat_in(:) ! Latitudes (ascending order, rad)
77	REAL, POINTER :: lev_in(:) ! Pressure levels (ascen. order, hPa)
78	REAL, ALLOCATABLE :: lon_in_edge(:) ! Longitude intervals edges
79	! (ascending order, / )
80	REAL, ALLOCATABLE :: sinlat_in_edge(:) ! Sinus of latitude intervals edges
81	! (ascending order, / )
82	LOGICAL :: ldec_lon, ldec_lat, ldec_lev ! Decreasing order in input file
83	CHARACTER(LEN=20) :: cal_in ! Calendar
84	REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:) ! Ozone climatologies
85	REAL, ALLOCATABLE :: o3_in2 (:,:,:,:) ! Ozone climatologies
86	! last index: 1 for the day-night average, 2 for the daylight field.
87	REAL :: NaN
88
89	!--- Partially or totally regridded variables (:,:,nlev_in,:,read_climoz)
90	REAL, ALLOCATABLE :: o3_regr_lon (:,:,:,:,:) ! (nlon_ou,nlat_in,:,0:13 ,:)
91	REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13 ,:)
92	REAL, ALLOCATABLE :: o3_out3 (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
93	REAL, ALLOCATABLE :: o3_regr_lat (:,:,:,:) ! (nlat_in,:,0:13 ,:)
94	REAL, ALLOCATABLE :: o3_out2 (:,:,:,:) ! (nlat_ou,:,ntim_ou,:)
95	! Dimension number \| Interval \| Contains \| For variables:
96	! 1 (longitude) \| [rlonu(i-1), rlonu(i)] \| rlonv(i) \| all
97	! 2 (latitude) \| [rlatv(j), rlatv(j-1)] \| rlatu(j) \| all but o3_regr_lon
98	! 3 (press level) \| \| lev(k) \| all
99	! Note that rlatv(0)=pi/2 and rlatv(nlat_ou)=-pi/2.
100	! Dimension 4 is: month number (all vars but o3_out)
101	! days elapsed since Jan. 1st 0h at mid-day (o3_out only)
102	REAL, ALLOCATABLE :: v1(:)
103
104	!--- For NetCDF:
105	INTEGER :: fID_in_m, fID_in, levID_ou, dimid, vID_in(read_climoz), ntim_ou
106	INTEGER :: fID_in_p, fID_ou, timID_ou, varid, vID_ou(read_climoz), ndims, ncerr
107	INTEGER, POINTER :: dIDs(:)
108	CHARACTER(LEN=20) :: cal_ou !--- Calendar; no time inter => same as input
109	CHARACTER(LEN=80) :: press_unit !--- Pressure unit
110	REAL :: tmidmonth(0:13) !--- Elapsed days since Jan-1 0h at mid-months
111	! Additional records 0, 13 for interpolation
112	REAL, ALLOCATABLE :: tmidday(:) !--- Output times (mid-days since Jan 1st 0h)
113	LOGICAL :: lprev, lnext !--- Flags: previous/next files are present
114	LOGICAL :: l3D, l2D !--- Flag: input fields are 3D or zonal
115	INTEGER :: ii, i, j, k, l, m, dln, ib, ie, iv, dx1, dx2
116	INTEGER, ALLOCATABLE :: sta(:), cnt(:)
117	CHARACTER(LEN=80) :: sub, dim_nam, msg
118	!-------------------------------------------------------------------------------
119	sub="regr_horiz_time_climoz"
120	WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
121	CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
122
123	CALL NF95_OPEN("climoz.nc" , NF90_NOWRITE, fID_in)
124	lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
125	lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
126
127	!--- Get coordinates from the input file. Converts lon/lat in radians.
128	! Few inversions because "regr_conserv" and gcm need ascending vectors.
129	CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
130	CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
131	l3D=ndims==4; l2D=ndims==3
132	IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
133	IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
134	DO i=1,ndims
135	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
136	CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
137	ii=i; IF(l2D) ii=i+1 !--- ndims==3:NO LONGITUDE
138	SELECT CASE(ii)
139	CASE(1) !--- LONGITUDE
140	CALL NF95_GW_VAR(fID_in, varid, lon_in)
141	ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
142	nlon_in=dln; lon_in=lon_in*deg2rad
143	CASE(2) !--- LATITUDE
144	CALL NF95_GW_VAR(fID_in, varid, lat_in)
145	ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
146	nlat_in=dln; lat_in=lat_in*deg2rad
147	CASE(3) !--- PRESSURE LEVELS
148	CALL NF95_GW_VAR(fID_in, varid, lev_in)
149	ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
150	nlev_in=dln
151	CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
152	k=LEN_TRIM(press_unit)
153	DO WHILE(ICHAR(press_unit(k:k))==0)
154	press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
155	END DO
156	IF(press_unit == "Pa") THEN
157	lev_in = lev_in/100. !--- CONVERT TO hPa
158	ELSE IF(press_unit /= "hPa") THEN
159	CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
160	END IF
161	CASE(4) !--- TIME
162	CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
163	cal_in='gregorian'
164	IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR) &
165	WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'// &
166	TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
167	k=LEN_TRIM(cal_in)
168	DO WHILE(ICHAR(cal_in(k:k))==0)
169	cal_in(k:k)=' '; k=LEN_TRIM(cal_in) !--- REMOVE NULL END CHAR
170	END DO
171	END SELECT
172	END DO
173
174	!--- Longitudes management:
175	! * Need to shift data if the origin of input file longitudes /= -pi
176	! * Need to add some margin in longitude to ensure input interval contains
177	! all the output intervals => at least one longitudes slice has to be
178	! duplicated, possibly more for undersampling.
179	IF(l3D) THEN
180	!--- Compute input edges longitudes vector (no end point yet)
181	ALLOCATE(v1(nlon_in))
182	v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
183	FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
184	DEALLOCATE(lon_in)
185
186	!--- Shift input longitudes vector until it contains first output point boundslon_reg(1,east)
187	v1=v1+2piREAL(FLOOR((boundslon_reg(1,east)-v1(1))/(2.*pi)))
188
189	!--- Ensure first input longitudes interval contains first output point boundslon_reg(1,east)
190	dx1=locate(v1,boundslon_reg(1,east))-1
191	v1=CSHIFT(v1,SHIFT=dx1,DIM=1); v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
192
193	!--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,west)
194	dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlat_ou,west)); dx2=dx2+1; END DO
195
196	!--- Final edges longitudes vector (with margin and end point)
197	ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(1:1+dx2)+2.*pi]
198	DEALLOCATE(v1)
199	END IF
200
201	!--- Compute sinus of intervals edges latitudes:
202	ALLOCATE(sinlat_in_edge(nlat_in+1))
203	sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
204	FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
205	DEALLOCATE(lat_in)
206
207	!--- Prepare quantities for time interpolation
208	tmidmonth=mid_month(anneeref, cal_in)
209	IF(interpt) THEN
210	ntim_ou=ioget_year_len(anneeref)
211	ALLOCATE(tmidday(ntim_ou))
212	tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
213	CALL ioget_calendar(cal_ou)
214	ELSE
215	ntim_ou=14
216	cal_ou=cal_in
217	END IF
218
219	!--- Create the output file and get the variable IDs:
220	CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
221	ndims, cal_ou)
222
223	!--- Write remaining coordinate variables:
224	CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
225	IF( interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
226	IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
227
228	!--- Check for contiguous years:
229	ib=0; ie=13
230	IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
231	WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
232	ELSE
233	IF( lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
234	IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
235	IF( lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
236	IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
237	IF(.NOT.lprev) ib=1
238	IF(.NOT.lnext) ie=12
239	END IF
240	ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
241	IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
242	IF(l2D) cnt=[ nlat_in,nlev_in,1]
243	IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
244	IF(l2D) ALLOCATE(o3_in2( nlat_in,nlev_in,ib:ie,read_climoz))
245
246	!--- Read full current file and one record each available contiguous file
247	DO iv=1,read_climoz
248	msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
249	CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
250	IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
251	IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2( :,:,1:12,iv))
252	CALL handle_err(TRIM(msg), ncerr, fID_in)
253	IF(lprev) THEN; sta(ndims)=12
254	CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
255	IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
256	IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2( :,:, 0,iv),sta,cnt)
257	CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
258	END IF
259	IF(lnext) THEN; sta(ndims)=1
260	CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
261	IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
262	IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2( :,:,13,iv),sta,cnt)
263	CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
264	END IF
265	END DO
266	IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
267	IF(lprev) CALL NF95_CLOSE(fID_in_m)
268	IF(lnext) CALL NF95_CLOSE(fID_in_p)
269
270	!--- Revert decreasing coordinates vector
271	IF(l3D) THEN
272	IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
273	IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
274	IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
275	!--- Shift values for longitude and duplicate some longitudes slices
276	o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
277	o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
278	ELSE
279	IF(ldec_lat) o3_in2 = o3_in2( nlat_in:1:-1,:,:,:)
280	IF(ldec_lev) o3_in2 = o3_in2( :,nlev_in:1:-1,:,:)
281	END IF
282
283	!--- Deal with missing values
284	DO m=1, read_climoz
285	WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
286	IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
287	IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
288	WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
289	END IF
290	END IF
291	WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
292	WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
293
294	!--- Check top layer contains no NaNs & search NaNs from top to ground
295	msg=TRIM(sub)//": NaNs in top layer !"
296	IF(l3D) THEN
297	IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
298	DO k = 2,nlev_in
299	WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
300	END DO
301	ELSE
302	IF(ANY(o3_in2( :,1,:,m)==NaN)) THEN
303	WRITE(lunout,*)msg
304	!--- Fill in latitudes where all values are missing
305	DO l=1,nmth_in
306	!--- Next to south pole
307	j=1; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
308	IF(j>1) &
309	o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
310	!--- Next to north pole
311	j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
312	IF(j<nlat_in) &
313	o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
314	END DO
315	END IF
316
317	!--- Fill in high latitudes missing values
318	!--- Highest level been filled-in, so has always valid values.
319	DO k = 2,nlev_in
320	WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
321	END DO
322	END IF
323	END DO
324	CALL NF95_CLOSE(fID_in)
325
326	!=============================================================================
327	IF(l3D) THEN !=== 3D FIELDS
328	!=============================================================================
329	!--- Regrid in longitude
330	ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
331	CALL regr_conserv(1, o3_in3, xs = lon_in_edge, &
332	xt = [boundslon_reg(1,east),boundslon_reg(:,west)], &
333	vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
334	DEALLOCATE(o3_in3)
335
336	!--- Regrid in latitude: averaging with respect to SIN(lat) is
337	! equivalent to weighting by COS(lat)
338	!--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
339	ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
340	CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge, &
341	xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
342	vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:), &
343	slope = slopes(2,o3_regr_lon, sinlat_in_edge))
344	DEALLOCATE(o3_regr_lon)
345
346	!--- Duplicate previous/next record(s) if they are not available
347	IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
348	IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
349
350	!--- Regrid in time by linear interpolation:
351	ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
352	IF( interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
353	IF(.NOT.interpt) o3_out3=o3_regr_lonlat
354	DEALLOCATE(o3_regr_lonlat)
355
356	!--- Write to file (the order of "rlatu" is inverted in the output file):
357	DO m = 1, read_climoz
358	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
359	END DO
360
361	!=============================================================================
362	ELSE !=== ZONAL FIELDS
363	!=============================================================================
364	!--- Regrid in latitude: averaging with respect to SIN(lat) is
365	! equivalent to weighting by COS(lat)
366	!--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
367	ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
368	CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge, &
369	xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
370	vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:), &
371	slope = slopes(1,o3_in2, sinlat_in_edge))
372	DEALLOCATE(o3_in2)
373
374	!--- Duplicate previous/next record(s) if they are not available
375	IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
376	IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
377
378	!--- Regrid in time by linear interpolation:
379	ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
380	IF( interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
381	IF(.NOT.interpt) o3_out2=o3_regr_lat
382	DEALLOCATE(o3_regr_lat)
383
384	!--- Write to file (the order of "rlatu" is inverted in the output file):
385	DO m = 1, read_climoz
386	CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
387	END DO
388
389	!=============================================================================
390	END IF
391	!=============================================================================
392
393	CALL NF95_CLOSE(fID_ou)
394
395	END SUBROUTINE regr_horiz_time_climoz
396	!
397	!-------------------------------------------------------------------------------
398
399
400	!-------------------------------------------------------------------------------
401	!
402	SUBROUTINE prepare_out(fID_in, nlev_in, ntim_ou, fID_ou, vlevID, vtimID, &
403	vID_ou, ndims, cal_ou)
404	!-------------------------------------------------------------------------------
405	! Purpose: This subroutine creates the NetCDF output file, defines
406	! dimensions and variables, and writes some of the coordinate variables.
407	!-------------------------------------------------------------------------------
408	USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
409	!-------------------------------------------------------------------------------
410	! Arguments:
411	INTEGER, INTENT(IN) :: fID_in, nlev_in, ntim_ou
412	INTEGER, INTENT(OUT) :: fID_ou, vlevID, vtimID
413	INTEGER, INTENT(OUT) :: vID_ou(:) ! dim(1/2) 1: O3day&night 2: O3daylight
414	INTEGER, INTENT(IN) :: ndims ! fields rank (3 or 4)
415	CHARACTER(LEN=*), INTENT(IN) :: cal_ou ! calendar
416	!-------------------------------------------------------------------------------
417	! Local variables:
418	INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
419	INTEGER :: vlonID, vlatID, ncerr, is
420	CHARACTER(LEN=80) :: sub
421	!-------------------------------------------------------------------------------
422	sub="prepare_out"
423	WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
424	CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
425
426	!--- Dimensions:
427	IF(ndims==4) &
428	CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
429	CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
430	CALL NF95_DEF_DIM(fID_ou, "plev", nlev_in, dlevID)
431	CALL NF95_DEF_DIM(fID_ou, "time", ntim_ou, dtimID)
432
433	!--- Define coordinate variables:
434	IF(ndims==4) &
435	CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
436	CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
437	CALL NF95_DEF_VAR(fID_ou, "plev", NF90_FLOAT, dlevID, vlevID)
438	CALL NF95_DEF_VAR(fID_ou, "time", NF90_FLOAT, dtimID, vtimID)
439	IF(ndims==4) &
440	CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
441	CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
442	CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
443	CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
444	IF(ndims==4) &
445	CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
446	CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
447	CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
448	CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
449	CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name", "air pressure")
450	CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar", cal_ou)
451
452	!--- Define the main variables:
453	IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
454	IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
455	CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
456	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
457	CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
458	&_in_air")
459	IF(SIZE(vID_ou) == 2) THEN
460	CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
461	CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
462	&ylight")
463	END IF
464
465	!--- Global attributes:
466	! The following commands, copying attributes, may fail. That is OK.
467	! It should just mean that the attribute is not defined in the input file.
468	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
469	CALL handle_err_copy_att("Conventions")
470	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title", fID_ou,NF90_GLOBAL, ncerr)
471	CALL handle_err_copy_att("title")
472	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
473	CALL handle_err_copy_att("institution")
474	CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source", fID_ou,NF90_GLOBAL, ncerr)
475	CALL handle_err_copy_att("source")
476	CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
477	CALL NF95_ENDDEF(fID_ou)
478
479	!--- Write one of the coordinate variables:
480	IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
481	CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
482	! (convert from rad to degrees and sort in ascending order)
483
484	CONTAINS
485
486	!-------------------------------------------------------------------------------
487	!
488	SUBROUTINE handle_err_copy_att(att_name)
489	!
490	!-------------------------------------------------------------------------------
491	USE netcdf, ONLY: NF90_NOERR, NF90_strerror
492	!-------------------------------------------------------------------------------
493	! Arguments:
494	CHARACTER(LEN=*), INTENT(IN) :: att_name
495	!-------------------------------------------------------------------------------
496	IF(ncerr /= NF90_NOERR) &
497	WRITE(lunout,*)TRIM(sub)//" prepare_out NF95_COPY_ATT "//TRIM(att_name)// &
498	" -- "//TRIM(NF90_strerror(ncerr))
499
500	END SUBROUTINE handle_err_copy_att
501	!
502	!-------------------------------------------------------------------------------
503
504	END SUBROUTINE prepare_out
505	!
506	!-------------------------------------------------------------------------------
507
508	END MODULE regr_horiz_time_climoz_m
509	!
510	!-------------------------------------------------------------------------------

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

Download in other formats:

Original Format