Context Navigation

regr_horiz_time_climoz_m.F90 @ 3776

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

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

Download in other formats:

Original Format