source: LMDZ6/trunk/libf/phylmd/regr_lat_time_coefoz_m.F90 @ 4484

Last change on this file since 4484 was 2788, checked in by dcugnet, 8 years ago

Changes in ce0l about the way ozone forcing files are generated:

1) 3D raw input files "climoz.nc" are now handled.
2) Default behaviour is now to let the gcm interpolate in time online.

This helps to avoid huge forcing files (in particular for 3D fields).
In this case, the output files "climoz_LMDZ.nc" all have 14 records:

  • records 2-13 are obtained with records 1-12 of "climoz.nc".
  • records 1 and 14 are obtained respectively with:
    • record 12 of "climoz_m.nc" if available, of "climoz.nc" otherwise.
    • record 1 of "climoz_p.nc" if available, of "climoz.nc" otherwise.

3) If ok_daily_climoz key is TRUE, the time interpolation (one record

a day) is forced, using the 14 records described below.
This now depends on the calendar (it was on a 360 days basis only).

Changes in the gcm about the way zone forcing files are read/interpolated:

1) 3D horizontally interpolated "climoz_LMDZ.nc" files are now handled.
2) Daily files (already interpolated in time) are still handled, but their

number of records must match the expected number of days, that depends
on the calendar (records step is no longer 1/360 year).

3) 14 records monthly files are now handled (and prefered). This reduces

the I/O to a minimum and the aditional computational cost is low (simple
online linear time interpolation).

4) If adjust_tropopause key is TRUE, the input fields are stretched using

following method:

  • LMDZ dynamical tropopause is detected: Ptrop_lmdz = MAX ( P(Potential Vorticity==2PVU), P(theta==380K) )
  • file chemical tropopause is detected: Ptrop_file = P( tro3 == o3t ), where:

o3t = 91. + 28. * SIN(PI*(month-2)/6) (ppbV)

This formula comes from Thouret & al., ACP 6, 1033-1051, 2006.
The second term of the expression is multiplied by TANH(lat_deg/20.)
to account for latitude dependency.

  • File profile is streched in a +/- 5kms zone around the mean tropopause to ensure resulting tropopause matches the one of LMDZ. See procedure regr_pr_time_av for more details.
  • 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: 12.0 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, handle_err, &
46         nf95_put_var, nf95_gw_var
47    use netcdf, only: nf90_nowrite, nf90_get_var
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, pointer:: 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, pointer:: 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) ! pointer
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) ! pointer
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       ncerr = nf90_get_var(ncid_in, varid_in(i_v), o3_par_in)
204       call handle_err("nf90_get_var", ncerr, ncid_in)
205
206       if (desc_lat) o3_par_in = o3_par_in(n_lat:1:-1, :, :)
207
208       ! Regrid in latitude:
209       ! We average with respect to sine of latitude, which is
210       ! equivalent to weighting by cosine of latitude:
211       call regr_conserv(1, o3_par_in, xs = sin(lat_in_edg), &
212            xt = (/-1., sin((/boundslat_reg(nbp_lat-1:1:-1,south)/)), 1./), &
213            vt = v_regr_lat(nbp_lat:1:-1, :, 1:12))
214       ! (invert order of indices in "v_regr_lat" because "rlatu" is
215       ! in descending order)
216
217       ! Duplicate January and December values, in preparation of time
218       ! interpolation:
219       v_regr_lat(:, :, 0) = v_regr_lat(:, :, 12)
220       v_regr_lat(:, :, 13) = v_regr_lat(:, :, 1)
221
222       ! Regrid in time by linear interpolation:
223       call regr_lint(3, v_regr_lat, tmidmonth, tmidday, o3_par_out)
224
225       ! Write to file:
226       call nf95_put_var(ncid_out, varid_out(i_v), &
227            o3_par_out(nbp_lat:1:-1, :, :))
228       ! (The order of "rlatu" is inverted in the output file)
229    end do
230
231    call nf95_close(ncid_out)
232    call nf95_close(ncid_in)
233
234  end subroutine regr_lat_time_coefoz
235
236  !********************************************
237
238  subroutine prepare_out(ncid_in, varid_in, n_plev, name_out, ncid_out, &
239       varid_out, varid_plev, varid_time)
240
241    ! This subroutine creates the NetCDF output file, defines
242    ! dimensions and variables, and writes one of the coordinate variables.
243
244    use mod_grid_phy_lmdz, ONLY : nbp_lat
245    use assert_eq_m, only: assert_eq
246
247    use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
248         nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
249    use netcdf, only: nf90_clobber, nf90_float, nf90_copy_att, nf90_global
250    use nrtype, only: pi
251    use regular_lonlat_mod, only : lat_reg
252
253    integer, intent(in):: ncid_in, varid_in(:), n_plev
254    character(len=*), intent(in):: name_out(:) ! of NetCDF variables
255    integer, intent(out):: ncid_out, varid_out(:), varid_plev, varid_time
256
257    ! Variables local to the procedure:
258
259    integer ncerr
260    integer dimid_rlatu, dimid_plev, dimid_time
261    integer varid_rlatu
262    integer i, n_o3_param
263
264    !---------------------------
265
266    print *, "Call sequence information: prepare_out"
267    n_o3_param = assert_eq(size(varid_in), size(varid_out), &
268         size(name_out), "prepare_out")
269
270    call nf95_create("coefoz_LMDZ.nc", nf90_clobber, ncid_out)
271
272    ! Dimensions:
273    call nf95_def_dim(ncid_out, "time", 360, dimid_time)
274    call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
275    call nf95_def_dim(ncid_out, "rlatu", nbp_lat, dimid_rlatu)
276
277    ! Define coordinate variables:
278
279    call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
280    call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
281    call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
282    call nf95_put_att(ncid_out, varid_time, "standard_name", "time")
283
284    call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
285    call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
286    call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
287    call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")
288
289    call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
290    call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
291    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
292
293    ! Define primary variables:
294
295    do i = 1, n_o3_param
296       call nf95_def_var(ncid_out, name_out(i), nf90_float, &
297            (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(i))
298
299       ! The following commands may fail. That is OK. It should just
300       ! mean that the attribute is not defined in the input file.
301
302       ncerr = nf90_copy_att(ncid_in, varid_in(i), "long_name",&
303            & ncid_out, varid_out(i))
304       call handle_err_copy_att("long_name")
305
306       ncerr = nf90_copy_att(ncid_in, varid_in(i), "units", ncid_out,&
307            & varid_out(i))
308       call handle_err_copy_att("units")
309
310       ncerr = nf90_copy_att(ncid_in, varid_in(i), "standard_name", ncid_out,&
311            & varid_out(i))
312       call handle_err_copy_att("standard_name")
313    end do
314
315    ! Global attributes:
316    call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, &
317         nf90_global)
318    call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global)
319    call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global)
320    call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
321
322    call nf95_enddef(ncid_out)
323
324    ! Write one of the coordinate variables:
325    call nf95_put_var(ncid_out, varid_rlatu, lat_reg(nbp_lat:1:-1) / pi * 180.)
326    ! (convert from rad to degrees and sort in ascending order)
327
328  contains
329
330    subroutine handle_err_copy_att(att_name)
331
332      use netcdf, only: nf90_noerr, nf90_strerror
333
334      character(len=*), intent(in):: att_name
335
336      !----------------------------------------
337
338      if (ncerr /= nf90_noerr) then
339         print *, "prepare_out " // trim(name_out(i)) &
340              // " nf90_copy_att " // att_name // " -- " &
341              // trim(nf90_strerror(ncerr))
342      end if
343
344    end subroutine handle_err_copy_att
345
346  end subroutine prepare_out
347
348end module regr_lat_time_coefoz_m
Note: See TracBrowser for help on using the repository browser.