/[lmdze]/trunk/libf/phylmd/Mobidic/regr_lat_time_coefoz.f90
ViewVC logotype

Annotation of /trunk/libf/phylmd/Mobidic/regr_lat_time_coefoz.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 52 - (hide annotations)
Fri Sep 23 12:28:01 2011 UTC (12 years, 8 months ago) by guez
File size: 11870 byte(s)
Split "conflx.f" into single-procedure files in directory "Conflx".

Split "cv_routines.f" into single-procedure files in directory
"CV_routines". Made module "cvparam" from included file
"cvparam.h". No included file other than "netcdf.inc" left in LMDZE.

1 guez 7 module regr_lat_time_coefoz_m
2 guez 3
3 guez 4 ! Author: Lionel GUEZ
4 guez 3
5     implicit none
6    
7     private
8 guez 9 public regr_lat_time_coefoz
9 guez 3
10     contains
11    
12 guez 7 subroutine regr_lat_time_coefoz
13 guez 3
14 guez 7 ! "regr_lat_time_coefoz" stands for "regrid latitude time
15     ! coefficients ozone".
16 guez 3
17 guez 7 ! This procedure reads from a NetCDF file parameters for ozone
18     ! chemistry, regrids them in latitude and time, and writes the
19     ! regridded fields to a new NetCDF file.
20    
21 guez 3 ! The input fields depend on time, pressure level and
22     ! latitude.
23     ! We assume that the input fields are step functions
24     ! of latitude.
25 guez 7 ! Regridding in latitude is made by averaging, with a cosine of
26     ! latitude factor.
27     ! The target LMDZ latitude grid is the "scalar" grid: "rlatu".
28 guez 3 ! The values of "rlatu" are taken to be the centers of intervals.
29 guez 7 ! Regridding in time is by linear interpolation.
30     ! Monthly values are processed to get daily values, on the basis
31     ! of a 360-day calendar.
32 guez 3
33 guez 7 ! We assume that in the input file:
34     ! -- the latitude is in degrees and strictly monotonic (as all
35     ! NetCDF coordinate variables should be);
36     ! -- time increases from January to December (even though we do
37     ! not use values of the input time coordinate).
38 guez 3
39 guez 7 use dimens_m, only: jjm
40 guez 3 use comgeom, only: rlatv
41 guez 39 use nr_util, only: pi
42 guez 52 use numer_rec, only: regr1_step_av, regr3_lint
43 guez 22 use netcdf95, only: nf95_open, nf95_gw_var, nf95_close, &
44 guez 7 nf95_inq_varid, handle_err, nf95_put_var
45     use netcdf, only: nf90_nowrite, nf90_get_var
46 guez 3
47     ! Variables local to the procedure:
48    
49     integer ncid_in, ncid_out ! NetCDF IDs for input and output files
50     integer n_plev ! number of pressure levels in the input data
51     integer n_lat! number of latitudes in the input data
52    
53     real, pointer:: latitude(:)
54     ! (of input data, converted to rad, sorted in strictly increasing order)
55    
56     real, allocatable:: lat_in_edg(:)
57     ! (edges of latitude intervals for input data, in rad, in strictly
58     ! increasing order)
59    
60 guez 7 real, pointer:: plev(:) ! pressure level of input data
61 guez 3 logical decr_lat ! decreasing latitude in the input file
62    
63 guez 7 real, allocatable:: o3_par_in(:, :, :) ! (n_lat, n_plev, 12)
64 guez 3 ! (ozone parameter from the input file)
65     ! ("o3_par_in(j, l, month)" is at latitude "latitude(j)" and pressure
66     ! level "plev(l)". "month" is between 1 and 12.)
67    
68 guez 7 real, allocatable:: v_regr_lat(:, :, :) ! (jjm + 1, n_plev, 0:13)
69     ! (mean of a variable "v" over a latitude interval)
70     ! (First dimension is latitude interval.
71     ! The latitude interval for "v_regr_lat(j,:, :)" contains "rlatu(j)".
72     ! If "j" is between 2 and "jjm" then the interval is:
73     ! [rlatv(j), rlatv(j-1)]
74     ! If "j" is 1 or "jjm + 1" then the interval is:
75     ! [rlatv(1), pi / 2]
76     ! or:
77     ! [- pi / 2, rlatv(jjm)]
78     ! respectively.
79     ! "v_regr_lat(:, l, :)" is for pressure level "plev(l)".
80     ! Last dimension is month number.)
81 guez 3
82 guez 7 real, allocatable:: o3_par_out(:, :, :) ! (jjm + 1, n_plev, 360)
83     ! (regridded ozone parameter)
84     ! ("o3_par_out(j, l, day)" is at latitude "rlatu(j)", pressure
85     ! level "plev(l)" and date "January 1st 0h" + "tmidday(day)", in a
86     ! 360-day calendar.)
87    
88     integer j
89 guez 3 integer i_v ! index of ozone parameter
90     integer, parameter:: n_o3_param = 8 ! number of ozone parameters
91    
92     character(len=11) name_in(n_o3_param)
93 guez 7 ! (name of NetCDF primary variable in the input file)
94 guez 4
95 guez 3 character(len=9) name_out(n_o3_param)
96 guez 7 ! (name of NetCDF primary variable in the output file)
97 guez 3
98 guez 7 integer varid_in(n_o3_param), varid_out(n_o3_param), varid_plev, varid_time
99 guez 22 integer ncerr, varid
100 guez 7 ! (for NetCDF)
101 guez 3
102 guez 7 real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * j, j = 0, 13)/)
103     ! (time to middle of month, in days since January 1st 0h, in a
104     ! 360-day calendar)
105     ! (We add values -15 and 375 so that, for example, day 3 of the year is
106     ! interpolated between the December and the January value.)
107 guez 3
108 guez 7 real, parameter:: tmidday(360) = (/(j + 0.5, j = 0, 359)/)
109     ! (time to middle of day, in days since January 1st 0h, in a
110     ! 360-day calendar)
111 guez 3
112     !---------------------------------
113    
114 guez 10 print *, "Call sequence information: regr_lat_time_coefoz"
115 guez 3
116 guez 7 ! Names of ozone parameters:
117 guez 3 i_v = 0
118    
119     i_v = i_v + 1
120     name_in(i_v) = "P_net"
121     name_out(i_v) = "P_net_Mob"
122    
123     i_v = i_v + 1
124     name_in(i_v) = "a2"
125     name_out(i_v) = "a2"
126    
127     i_v = i_v + 1
128 guez 23 name_in(i_v) = "tro3"
129 guez 3 name_out(i_v) = "r_Mob"
130    
131     i_v = i_v + 1
132     name_in(i_v) = "a4"
133     name_out(i_v) = "a4"
134    
135     i_v = i_v + 1
136     name_in(i_v) = "temperature"
137     name_out(i_v) = "temp_Mob"
138    
139     i_v = i_v + 1
140     name_in(i_v) = "a6"
141     name_out(i_v) = "a6"
142    
143     i_v = i_v + 1
144     name_in(i_v) = "Sigma"
145     name_out(i_v) = "Sigma_Mob"
146    
147     i_v = i_v + 1
148     name_in(i_v) = "R_Het"
149     name_out(i_v) = "R_Het"
150    
151 guez 22 call nf95_open("coefoz.nc", nf90_nowrite, ncid_in)
152 guez 3
153     ! Get coordinates from the input file:
154    
155 guez 22 call nf95_inq_varid(ncid_in, "latitude", varid)
156     call nf95_gw_var(ncid_in, varid, latitude)
157 guez 3 ! Convert from degrees to rad, because "rlatv" is in rad:
158 guez 4 latitude = latitude / 180. * pi
159 guez 3 n_lat = size(latitude)
160 guez 7 ! We need to supply the latitudes to "regr1_step_av" in
161 guez 3 ! increasing order, so invert order if necessary:
162     decr_lat = latitude(1) > latitude(n_lat)
163 guez 4 if (decr_lat) latitude = latitude(n_lat:1:-1)
164 guez 3
165     ! Compute edges of latitude intervals:
166     allocate(lat_in_edg(n_lat + 1))
167     lat_in_edg(1) = - pi / 2
168     forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
169     lat_in_edg(n_lat + 1) = pi / 2
170     deallocate(latitude) ! pointer
171    
172 guez 22 call nf95_inq_varid(ncid_in, "plev", varid)
173     call nf95_gw_var(ncid_in, varid, plev)
174 guez 3 n_plev = size(plev)
175 guez 7 ! (We only need the pressure coordinate to copy it to the output file.)
176 guez 3
177     ! Get the IDs of ozone parameters in the input file:
178     do i_v = 1, n_o3_param
179     call nf95_inq_varid(ncid_in, trim(name_in(i_v)), varid_in(i_v))
180     end do
181    
182 guez 7 ! Create the output file and get the variable IDs:
183     call prepare_out(ncid_in, varid_in, n_plev, name_out, ncid_out, &
184     varid_out, varid_plev, varid_time)
185 guez 3
186 guez 7 ! Write remaining coordinate variables:
187     call nf95_put_var(ncid_out, varid_time, tmidday)
188     call nf95_put_var(ncid_out, varid_plev, plev)
189    
190     deallocate(plev) ! pointer
191    
192     allocate(o3_par_in(n_lat, n_plev, 12))
193     allocate(v_regr_lat(jjm + 1, n_plev, 0:13))
194     allocate(o3_par_out(jjm + 1, n_plev, 360))
195    
196 guez 3 do i_v = 1, n_o3_param
197     ! Process ozone parameter "name_in(i_v)"
198    
199     ncerr = nf90_get_var(ncid_in, varid_in(i_v), o3_par_in)
200     call handle_err("nf90_get_var", ncerr, ncid_in)
201    
202     if (decr_lat) o3_par_in = o3_par_in(n_lat:1:-1, :, :)
203    
204     ! Regrid in latitude:
205     ! We average with respect to sine of latitude, which is
206     ! equivalent to weighting by cosine of latitude:
207 guez 7 v_regr_lat(jjm+1:1:-1, :, 1:12) = regr1_step_av(o3_par_in, &
208 guez 3 xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
209 guez 4 ! (invert order of indices in "v_regr_lat" because "rlatu" is
210     ! decreasing)
211 guez 3
212 guez 7 ! Duplicate January and December values, in preparation of time
213     ! interpolation:
214     v_regr_lat(:, :, 0) = v_regr_lat(:, :, 12)
215     v_regr_lat(:, :, 13) = v_regr_lat(:, :, 1)
216 guez 3
217 guez 7 ! Regrid in time by linear interpolation:
218     o3_par_out = regr3_lint(v_regr_lat, tmidmonth, tmidday)
219 guez 3
220     ! Write to file:
221 guez 7 call nf95_put_var(ncid_out, varid_out(i_v), &
222     o3_par_out(jjm+1:1:-1, :, :))
223 guez 3 ! (The order of "rlatu" is inverted in the output file)
224     end do
225    
226     call nf95_close(ncid_out)
227     call nf95_close(ncid_in)
228    
229 guez 7 end subroutine regr_lat_time_coefoz
230 guez 3
231     !********************************************
232    
233 guez 7 subroutine prepare_out(ncid_in, varid_in, n_plev, name_out, ncid_out, &
234     varid_out, varid_plev, varid_time)
235 guez 3
236     ! This subroutine creates the NetCDF output file, defines
237 guez 7 ! dimensions and variables, and writes one of the coordinate variables.
238 guez 3
239 guez 7 use dimens_m, only: jjm
240     use comgeom, only: rlatu
241 guez 39 use nr_util, only: assert_eq, pi
242 guez 3
243 guez 7 use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
244     nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
245     use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global
246 guez 3
247 guez 7 integer, intent(in):: ncid_in, varid_in(:), n_plev
248 guez 3 character(len=*), intent(in):: name_out(:) ! of NetCDF variables
249 guez 7 integer, intent(out):: ncid_out, varid_out(:), varid_plev, varid_time
250 guez 3
251     ! Variables local to the procedure:
252     integer ncerr
253 guez 7 integer dimid_rlatu, dimid_plev, dimid_time
254     integer varid_rlatu
255 guez 3 integer i, n_o3_param
256    
257     !---------------------------
258    
259     print *, "Call sequence information: prepare_out"
260     n_o3_param = assert_eq(size(varid_in), size(varid_out), &
261     size(name_out), "prepare_out")
262    
263     call nf95_create("coefoz_LMDZ.nc", nf90_clobber, ncid_out)
264    
265     ! Dimensions:
266 guez 7 call nf95_def_dim(ncid_out, "time", 360, dimid_time)
267     call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
268 guez 3 call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu)
269    
270 guez 7 ! Define coordinate variables:
271 guez 3
272 guez 7 call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
273     call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
274     call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
275     call nf95_put_att(ncid_out, varid_time, "standard_name", "time")
276 guez 3
277 guez 7 call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
278     call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
279     call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
280     call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")
281 guez 3
282     call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
283     call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
284 guez 7 call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
285 guez 3
286 guez 7 ! Define primary variables:
287 guez 3
288     do i = 1, n_o3_param
289     call nf95_def_var(ncid_out, name_out(i), nf90_float, &
290 guez 7 (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(i))
291    
292 guez 3 ! The following commands may fail. That is OK. It should just
293     ! mean that the attribute is not defined in the input file.
294 guez 7
295 guez 3 ncerr = nf90_copy_att(ncid_in, varid_in(i), "long_name",&
296     & ncid_out, varid_out(i))
297     call handle_err_copy_att("long_name")
298 guez 7
299 guez 3 ncerr = nf90_copy_att(ncid_in, varid_in(i), "units", ncid_out,&
300     & varid_out(i))
301     call handle_err_copy_att("units")
302 guez 7
303 guez 3 ncerr = nf90_copy_att(ncid_in, varid_in(i), "standard_name", ncid_out,&
304     & varid_out(i))
305     call handle_err_copy_att("standard_name")
306     end do
307    
308     ! Global attributes:
309 guez 7 call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, &
310     nf90_global)
311 guez 3 call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global)
312     call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global)
313     call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
314    
315     call nf95_enddef(ncid_out)
316    
317 guez 7 ! Write one of the coordinate variables:
318     call nf95_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.)
319 guez 3 ! (convert from rad to degrees and sort in increasing order)
320    
321     contains
322    
323     subroutine handle_err_copy_att(att_name)
324    
325 guez 7 use netcdf, only: nf90_noerr, nf90_strerror
326 guez 3
327     character(len=*), intent(in):: att_name
328    
329     !----------------------------------------
330    
331     if (ncerr /= nf90_noerr) then
332     print *, "prepare_out " // trim(name_out(i)) &
333     // " nf90_copy_att " // att_name // " -- " &
334     // trim(nf90_strerror(ncerr))
335     end if
336    
337     end subroutine handle_err_copy_att
338    
339     end subroutine prepare_out
340    
341 guez 7 end module regr_lat_time_coefoz_m

  ViewVC Help
Powered by ViewVC 1.1.21