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

  ViewVC Help
Powered by ViewVC 1.1.21