source: LMDZ6/branches/contrails/libf/phylmd/regr_lat_time_coefoz_m.f90 @ 5426

Last change on this file since 5426 was 5268, checked in by abarral, 2 months ago

.f90 <-> .F90 depending on cpp key use

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