/[lmdze]/trunk/libf/dyn3d/regr_coefoz.f90
ViewVC logotype

Annotation of /trunk/libf/dyn3d/regr_coefoz.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4 - (hide annotations)
Thu Feb 28 18:05:06 2008 UTC (16 years, 2 months ago) by guez
File size: 15577 byte(s)
Cosmetic changes
1 guez 3 module regr_coefoz_m
2    
3     ! This module is clean: no C preprocessor directive, no include line.
4 guez 4 ! Author: Lionel GUEZ
5 guez 3
6     implicit none
7    
8     private
9     public regr_coefoz
10    
11     contains
12    
13     subroutine regr_coefoz
14    
15     ! "regr_coefoz" stands for "regrid coefficients ozone".
16    
17     ! This procedure reads from a file parameters for ozone
18     ! chemistry and regrids them for LMDZ.
19     ! The input fields depend on time, pressure level and
20     ! latitude.
21     ! We assume that the input fields are step functions
22     ! of latitude.
23     ! Horizontal regridding is made by averaging on latitude, with a
24     ! cosine of latitude factor.
25     ! The target horizontal LMDZ grid is the "scalar" grid: "rlonv", "rlatu".
26     ! The values of "rlatu" are taken to be the centers of intervals.
27     ! The target vertical LMDZ grid is the grid of mid-layers.
28     ! The input latitude values are different from the latitude values
29     ! of the LMDZ "scalar" grid.
30     ! The input data does not depend on longitude, but the pressure
31     ! at LMDZ mid-layers does.
32     ! Therefore, the values on the LMDZ grid do depend on longitude.
33     ! The regridded fields are written to a file.
34    
35     ! We assume that in the input file the latitude is in degrees
36     ! and the pressure level is in hPa, and that both are strictly
37 guez 4 ! monotonic (as all NetCDF coordinate variables should be).
38 guez 3
39     use dimens_m, only: iim, jjm, llm
40     use comgeom, only: rlatv
41     use comconst, only: pi
42     use pressure_m, only: p3d, pls
43     use regr1_step_av_m, only: regr1_step_av
44     use regr1_lint_m, only: regr1_lint
45     use netcdf95, only: nf95_open, nf90_nowrite, coordin, nf95_close, &
46     nf95_inq_varid, handle_err, nf90_put_var, nf90_get_var
47    
48     ! Variables local to the procedure:
49    
50     integer ncid_in, ncid_out ! NetCDF IDs for input and output files
51    
52     integer n_plev ! number of pressure levels in the input data
53     integer n_lat! number of latitudes in the input data
54    
55     real, pointer:: latitude(:)
56     ! (of input data, converted to rad, sorted in strictly increasing order)
57    
58     real, allocatable:: lat_in_edg(:)
59     ! (edges of latitude intervals for input data, in rad, in strictly
60     ! increasing order)
61    
62     real, pointer:: plev(:)
63     ! (pressure level of input data, converted to Pa, sorted
64     ! in strictly increasing order)
65    
66     real, allocatable:: press_in_edg(:)
67     ! (edges of pressure intervals for input data, in Pa, in strictly
68     ! increasing order)
69    
70     logical decr_lat ! decreasing latitude in the input file
71     logical decr_plev ! decreasing pressure level in the input file
72    
73     real, allocatable:: o3_par_in(:, :, :)
74     ! (ozone parameter from the input file)
75     ! ("o3_par_in(j, l, month)" is at latitude "latitude(j)" and pressure
76     ! level "plev(l)". "month" is between 1 and 12.)
77    
78     real o3_par_out(iim + 1, jjm + 1, llm, 12)
79     ! (ozone parameter adapted to the LMDZ grid)
80     ! (Last dimension is month number.
81     ! "o3_par_out(i, j, l, month)" is at longitude "rlonv(i)", latitude
82     ! "rlatu(j)" and pressure level "pls(i, j, l)")
83    
84     integer i, j
85     integer i_v ! index of ozone parameter
86     integer, parameter:: n_o3_param = 8 ! number of ozone parameters
87    
88     character(len=11) name_in(n_o3_param)
89 guez 4 ! (name of NetCDF variable in the input file)
90    
91 guez 3 character(len=9) name_out(n_o3_param)
92 guez 4 ! (name of NetCDF variable in the output file)
93 guez 3
94     logical:: stepav_choice(n_o3_param) = .true.
95 guez 4 ! (vertical regridding by step average, otherwise linear interpolation)
96 guez 3
97     real top_value(n_o3_param)
98     ! (value at 0 pressure, only used for linear interpolation)
99    
100     integer varid_in(n_o3_param), varid_out(n_o3_param), ncerr ! for NetCDF
101    
102 guez 4 real, allocatable:: v_regr_lat(:, :, :) ! (jjm + 1, 0:n_plev, 12)
103     ! (mean of a variable "v" over a latitude interval)
104 guez 3 ! First dimension is latitude interval.
105     ! The latitude interval for "v_regr_lat(j,:, :)" contains "rlatu(j)".
106     ! If "j" is between 2 and "jjm" then the interval is:
107     ! [rlatv(j), rlatv(j-1)]
108     ! If "j" is 1 or "jjm + 1" then the interval is:
109     ! [rlatv(1), pi / 2]
110     ! or:
111     ! [- pi / 2, rlatv(jjm)]
112     ! respectively.
113     ! "v_regr_lat(:, l, :)" is for pressure interval
114     ! "[press_in_edg(l), press_in_edg(l+1)]" or pressure level "plev(l)",
115     ! depending on the type of vertical regridding, step average or linear
116     ! interpolation.
117     ! Last dimension is month number.)
118    
119     !---------------------------------
120    
121     print *, "Call sequence information: regr_coefoz"
122    
123     ! Details for each ozone parameter:
124     i_v = 0
125    
126     i_v = i_v + 1
127     name_in(i_v) = "P_net"
128     name_out(i_v) = "P_net_Mob"
129    
130     i_v = i_v + 1
131     name_in(i_v) = "a2"
132     name_out(i_v) = "a2"
133    
134     i_v = i_v + 1
135     name_in(i_v) = "r"
136     name_out(i_v) = "r_Mob"
137    
138     i_v = i_v + 1
139     name_in(i_v) = "a4"
140     name_out(i_v) = "a4"
141    
142     i_v = i_v + 1
143     name_in(i_v) = "temperature"
144     name_out(i_v) = "temp_Mob"
145    
146     i_v = i_v + 1
147     name_in(i_v) = "a6"
148     name_out(i_v) = "a6"
149    
150     i_v = i_v + 1
151     name_in(i_v) = "Sigma"
152     name_out(i_v) = "Sigma_Mob"
153     stepav_choice(i_v) = .false.
154     top_value(i_v) = 0.
155    
156     i_v = i_v + 1
157     name_in(i_v) = "R_Het"
158     name_out(i_v) = "R_Het"
159    
160     call nf95_open("coefoz_v2_3.nc", nf90_nowrite, ncid_in)
161    
162     ! Get coordinates from the input file:
163    
164     latitude => coordin(ncid_in, "latitude")
165     ! Convert from degrees to rad, because "rlatv" is in rad:
166 guez 4 latitude = latitude / 180. * pi
167 guez 3 n_lat = size(latitude)
168     ! We need to supply the latitudes to "stepav" in
169     ! increasing order, so invert order if necessary:
170     decr_lat = latitude(1) > latitude(n_lat)
171 guez 4 if (decr_lat) latitude = latitude(n_lat:1:-1)
172 guez 3
173     ! Compute edges of latitude intervals:
174     allocate(lat_in_edg(n_lat + 1))
175     lat_in_edg(1) = - pi / 2
176     forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
177     lat_in_edg(n_lat + 1) = pi / 2
178     deallocate(latitude) ! pointer
179    
180     plev => coordin(ncid_in, "plev")
181 guez 4 ! Convert from hPa to Pa because "p3d" and "pls" are in Pa:
182     plev = plev * 100.
183 guez 3 n_plev = size(plev)
184     ! We need to supply the pressure levels to "stepav" in
185     ! increasing order, so invert order if necessary:
186     decr_plev = plev(1) > plev(n_plev)
187 guez 4 if (decr_plev) plev = plev(n_plev:1:-1)
188 guez 3
189     ! Compute edges of pressure intervals:
190     allocate(press_in_edg(n_plev + 1))
191     press_in_edg(1) = 0.
192     ! We choose edges halfway in logarithm:
193     forall (j = 2:n_plev) press_in_edg(j) = sqrt(plev(j - 1) * plev(j))
194     press_in_edg(n_plev + 1) = huge(0.)
195     ! (infinity, but any value guaranteed to be greater than the
196     ! surface pressure would do)
197    
198     ! Get the IDs of ozone parameters in the input file:
199     do i_v = 1, n_o3_param
200     call nf95_inq_varid(ncid_in, trim(name_in(i_v)), varid_in(i_v))
201     end do
202    
203     call prepare_out(ncid_in, varid_in, name_out, ncid_out, varid_out)
204     allocate(o3_par_in(n_lat, n_plev, 12), v_regr_lat(jjm + 1, 0:n_plev, 12))
205    
206     do i_v = 1, n_o3_param
207     ! Process ozone parameter "name_in(i_v)"
208    
209     ncerr = nf90_get_var(ncid_in, varid_in(i_v), o3_par_in)
210     call handle_err("nf90_get_var", ncerr, ncid_in)
211    
212     if (decr_lat) o3_par_in = o3_par_in(n_lat:1:-1, :, :)
213     if (decr_plev) o3_par_in = o3_par_in(:, n_plev:1:-1, :)
214    
215     ! Regrid in latitude:
216     ! We average with respect to sine of latitude, which is
217     ! equivalent to weighting by cosine of latitude:
218     v_regr_lat(jjm+1:1:-1, 1:, :) = regr1_step_av(o3_par_in, &
219     xs=sin(lat_in_edg), xt=sin((/- pi / 2, rlatv(jjm:1:-1), pi / 2/)))
220 guez 4 ! (invert order of indices in "v_regr_lat" because "rlatu" is
221     ! decreasing)
222 guez 3
223     ! Regrid in pressure at each horizontal position:
224    
225     if (stepav_choice(i_v)) then
226     ! Regrid in pressure by averaging a step function of pressure
227    
228     ! Poles ("p3d" does not depend on longitude):
229     do j = 1, jjm + 1, jjm
230     ! Average on pressure, only for first longitude:
231     o3_par_out(1, j, llm:1:-1, :) &
232     = regr1_step_av(v_regr_lat(j, 1:, :), press_in_edg, &
233     p3d(1, j, llm+1:1:-1))
234     ! (invert order of indices because "p3d" is decreasing)
235     end do
236    
237     ! Latitudes other than poles ("p3d" depends on longitude):
238     do j = 2, jjm
239     ! Average on pressure at each longitude:
240     do i = 1, iim
241     o3_par_out(i, j, llm:1:-1, :) &
242     = regr1_step_av(v_regr_lat(j, 1:, :), press_in_edg, &
243     p3d(i, j, llm+1:1:-1))
244     ! (invert order of indices because "p3d" is decreasing)
245     end do
246     end do
247     else
248     ! Regrid in pressure by linear interpolation
249    
250     ! Complete "v_regr_lat" with the value at 0 pressure:
251     v_regr_lat(:, 0, :) = top_value(i_v)
252    
253     ! Poles ("pls" does not depend on longitude):
254     do j = 1, jjm + 1, jjm
255     ! Interpolate in pressure, only for first longitude:
256     o3_par_out(1, j, llm:1:-1, :) = regr1_lint(v_regr_lat(j, :, :), &
257     xs=(/0., plev/), xt=pls(1, j, llm:1:-1))
258     ! (invert order of indices because "pls" is decreasing)
259     end do
260    
261     ! Latitudes other than poles ("pls" depends on longitude):
262     do j = 2, jjm
263     ! Average on pressure at each longitude:
264     do i = 1, iim
265     o3_par_out(i, j, llm:1:-1, :) &
266     = regr1_lint(v_regr_lat(j, :, :), xs=(/0., plev/), &
267     xt=pls(i, j, llm:1:-1))
268     ! (invert order of indices because "pls" is decreasing)
269     end do
270     end do
271     end if
272    
273     ! Duplicate pole values on all longitudes:
274     o3_par_out(2:, 1, :, :) &
275     = spread(o3_par_out(1, 1, :, :), dim=1, ncopies=iim)
276     o3_par_out(2:, jjm + 1, :, :) &
277     = spread(o3_par_out(1, jjm + 1, :, :), dim=1, ncopies=iim)
278    
279     ! Duplicate first longitude to last longitude:
280     o3_par_out(iim + 1, 2:jjm, :, :) = o3_par_out(1, 2:jjm, :, :)
281    
282     ! Write to file:
283    
284     ncerr = nf90_put_var(ncid_out, varid_out(i_v), &
285     o3_par_out(:,jjm+1:1:-1, :, :))
286     ! (The order of "rlatu" is inverted in the output file)
287     call handle_err("nf90_put_var", ncerr, ncid_out)
288     end do
289    
290     deallocate(plev) ! pointer
291     call nf95_close(ncid_out)
292     call nf95_close(ncid_in)
293    
294     end subroutine regr_coefoz
295    
296     !********************************************
297    
298     subroutine prepare_out(ncid_in, varid_in, name_out, ncid_out, varid_out)
299    
300     ! This subroutine creates the NetCDF output file, defines
301 guez 4 ! dimensions and variables, and writes coordinate variables and "pls".
302 guez 3
303     use dimens_m, only: iim, jjm, llm
304     use comgeom, only: rlatu, rlonv
305     use pressure_m, only: pls
306     use comconst, only: pi
307     use nrutil, only: assert_eq
308    
309     use netcdf95, only: nf95_create, nf90_clobber, nf95_def_dim, &
310     nf95_def_var, nf90_float, nf90_int, nf95_put_att, nf95_enddef, &
311     nf90_put_var, handle_err, nf90_copy_att, nf95_copy_att, nf90_global
312    
313     integer, intent(in):: ncid_in, varid_in(:)
314     character(len=*), intent(in):: name_out(:) ! of NetCDF variables
315     integer, intent(out):: ncid_out, varid_out(:)
316    
317     ! Variables local to the procedure:
318     integer ncerr
319     integer dimid_rlonv, dimid_rlatu, dimid_sigs, dimid_month
320     integer varid_rlonv, varid_rlatu, varid_sigs, varid_month
321     integer varid_pls
322     integer i, n_o3_param
323    
324     !---------------------------
325    
326     print *, "Call sequence information: prepare_out"
327     n_o3_param = assert_eq(size(varid_in), size(varid_out), &
328     size(name_out), "prepare_out")
329    
330     call nf95_create("coefoz_LMDZ.nc", nf90_clobber, ncid_out)
331    
332     ! Dimensions:
333     call nf95_def_dim(ncid_out, "month", 12, dimid_month)
334     call nf95_def_dim(ncid_out, "sigs", llm, dimid_sigs)
335     call nf95_def_dim(ncid_out, "rlatu", jjm + 1, dimid_rlatu)
336     call nf95_def_dim(ncid_out, "rlonv", iim + 1, dimid_rlonv)
337    
338     ! Coordinate variables:
339    
340     call nf95_def_var(ncid_out, "month", nf90_float, dimid_month, varid_month)
341     call nf95_put_att(ncid_out, varid_month, "units", &
342     "calendar_month as %m.%f")
343     call nf95_put_att(ncid_out, varid_month, "long_name", "seasonal phase")
344    
345     call nf95_def_var(ncid_out, "sigs", nf90_int, dimid_sigs, varid_sigs)
346     call nf95_put_att(ncid_out, varid_sigs, "long_name", "s-layer index")
347    
348     call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
349     call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
350     call nf95_put_att(ncid_out, varid_rlatu, "long_name", "latitude")
351    
352     call nf95_def_var(ncid_out, "rlonv", nf90_float, dimid_rlonv, varid_rlonv)
353     call nf95_put_att(ncid_out, varid_rlonv, "units", "degrees_east")
354     call nf95_put_att(ncid_out, varid_rlonv, "long_name", "longitude")
355    
356     ! Primary variables:
357    
358     call nf95_def_var(ncid_out, "pls", nf90_float, &
359     (/dimid_rlonv, dimid_rlatu, dimid_sigs/), varid_pls)
360     call nf95_put_att(ncid_out, varid_pls, "units", "millibar")
361     call nf95_put_att(ncid_out, varid_pls, "long_name", &
362     "pressure at LMDZ mid-layers")
363    
364     do i = 1, n_o3_param
365     call nf95_def_var(ncid_out, name_out(i), nf90_float, &
366     (/dimid_rlonv, dimid_rlatu, dimid_sigs, dimid_month/),&
367     & varid_out(i))
368     ! The following commands may fail. That is OK. It should just
369     ! mean that the attribute is not defined in the input file.
370     ncerr = nf90_copy_att(ncid_in, varid_in(i), "long_name",&
371     & ncid_out, varid_out(i))
372     call handle_err_copy_att("long_name")
373     ncerr = nf90_copy_att(ncid_in, varid_in(i), "units", ncid_out,&
374     & varid_out(i))
375     call handle_err_copy_att("units")
376     ncerr = nf90_copy_att(ncid_in, varid_in(i), "standard_name", ncid_out,&
377     & varid_out(i))
378     call handle_err_copy_att("standard_name")
379     end do
380    
381     ! Global attributes:
382     call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global)
383     call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global)
384     call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
385    
386     call nf95_enddef(ncid_out)
387    
388     ! Coordinate variables:
389    
390     ncerr = nf90_put_var(ncid_out, varid_month, (/(i + 0.5, i = 1, 12)/))
391     call handle_err("nf90_put_var", ncerr, ncid_out)
392    
393     ncerr = nf90_put_var(ncid_out, varid_sigs, (/(i, i = 1, llm)/))
394     call handle_err("nf90_put_var", ncerr, ncid_out)
395    
396     ncerr = nf90_put_var(ncid_out, varid_rlatu, rlatu(jjm+1:1:-1) / pi * 180.)
397     ! (convert from rad to degrees and sort in increasing order)
398     call handle_err("nf90_put_var", ncerr, ncid_out)
399    
400 guez 4 ncerr = nf90_put_var(ncid_out, varid_rlonv, rlonv / pi * 180.)
401 guez 3 ! (convert from rad to degrees)
402     call handle_err("nf90_put_var", ncerr, ncid_out)
403    
404     ! Primary variable:
405    
406     ncerr = nf90_put_var(ncid_out, varid_pls, pls(:, jjm+1:1:-1, :) / 100.)
407     ! (convert from Pa to hPa)
408     call handle_err("nf90_put_var", ncerr, ncid_out)
409    
410     contains
411    
412     subroutine handle_err_copy_att(att_name)
413    
414     use netcdf95, only: nf90_noerr, nf90_strerror
415    
416     character(len=*), intent(in):: att_name
417    
418     !----------------------------------------
419    
420     if (ncerr /= nf90_noerr) then
421     print *, "prepare_out " // trim(name_out(i)) &
422     // " nf90_copy_att " // att_name // " -- " &
423     // trim(nf90_strerror(ncerr))
424     end if
425    
426     end subroutine handle_err_copy_att
427    
428     end subroutine prepare_out
429    
430     end module regr_coefoz_m

  ViewVC Help
Powered by ViewVC 1.1.21