/[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 61 - (hide annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years, 2 months ago) by guez
File size: 11873 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

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 61 use numer_rec_95, 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