Context Navigation

regr_horiz_time_climoz_m.F90 @ 2908

Last change on this file since 2908 was 2788, checked in by dcugnet, 8 years ago

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

records 2-13 are obtained with records 1-12 of "climoz.nc".
records 1 and 14 are obtained respectively with:
- record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
- record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.

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