source: LMDZ5/branches/AI-cosp/libf/phylmd/regr_lat_time_climoz_m.F90 @ 5360

Last change on this file since 5360 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

  • 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: 16.5 KB
Line 
1! $Id$
2module regr_lat_time_climoz_m
3
4  ! Author: Lionel GUEZ
5
6  implicit none
7
8  private
9  public regr_lat_time_climoz
10
11contains
12
13  subroutine regr_lat_time_climoz(read_climoz)
14
15    ! "regr_lat_time_climoz" stands for "regrid latitude time
16    ! climatology ozone".
17
18    ! This procedure reads a climatology of ozone from a NetCDF file,
19    ! regrids it in latitude and time, and writes the regridded field
20    ! to a new NetCDF file.
21
22    ! The input field depends on time, pressure level and latitude.
23
24    ! If the input field has missing values, they must be signaled by
25    ! the "missing_value" attribute.
26
27    ! We assume that the input field is a step function of latitude
28    ! and that the input latitude coordinate gives the centers of steps.
29    ! Regridding in latitude is made by averaging, with a cosine of
30    ! latitude factor.
31    ! The target LMDZ latitude grid is the "scalar" grid: "rlatu".
32    ! The values of "rlatu" are taken to be the centers of intervals.
33
34    ! We assume that in the input file:
35
36    ! -- Latitude is in degrees.
37
38    ! -- Latitude and pressure are strictly monotonic (as all NetCDF
39    ! coordinate variables should be).
40
41    ! -- The time coordinate is in ascending order (even though we do
42    ! not use its values).
43    ! The input file may contain either values for 12 months or values
44    ! for 14 months.
45    ! If there are 14 months then we assume that we have (in that order):
46    ! December, January, February, ..., November, December, January
47
48    ! -- Missing values are contiguous, at the bottom of
49    ! the vertical domain and at the latitudinal boundaries.
50
51    ! If values are all missing at a given latitude and date, then we
52    ! replace those missing values by values at the closest latitude,
53    ! equatorward, with valid values.
54    ! Then, at each latitude and each date, the missing values are replaced
55    ! by the lowest valid value above missing values.
56
57    ! Regridding in time is by linear interpolation.
58    ! Monthly values are processed to get daily values, on the basis
59    ! of a 360-day calendar.
60    ! If there are 14 months, we use the first December value to
61    ! interpolate values between January 1st and mid-January.
62    ! We use the last January value to interpolate values between
63    ! mid-December and end of December.
64    ! If there are only 12 months in the input file then we assume
65    ! periodicity for interpolation at the beginning and at the end of the
66    ! year.
67
68    use mod_grid_phy_lmdz, ONLY : nbp_lat
69    use regr1_step_av_m, only: regr1_step_av
70    use regr3_lint_m, only: regr3_lint
71    use netcdf95, only: handle_err, nf95_close, nf95_get_att, nf95_gw_var, &
72         nf95_inq_dimid, nf95_inq_varid, nf95_inquire_dimension, nf95_open, &
73         nf95_put_var
74    use netcdf, only: nf90_get_att, nf90_get_var, nf90_noerr, nf90_nowrite
75    use assert_m, only: assert
76    use regular_lonlat_mod, only : boundslat_reg, south
77    use nrtype, only: pi
78
79    integer, intent(in):: read_climoz ! read ozone climatology
80    ! Allowed values are 1 and 2
81    ! 1: read a single ozone climatology that will be used day and night
82    ! 2: read two ozone climatologies, the average day and night
83    ! climatology and the daylight climatology
84
85    ! Variables local to the procedure:
86
87    integer n_plev ! number of pressure levels in the input data
88    integer n_lat ! number of latitudes in the input data
89    integer n_month ! number of months in the input data
90
91    real, pointer:: latitude(:)
92    ! (of input data, converted to rad, sorted in strictly ascending order)
93
94    real, allocatable:: lat_in_edg(:)
95    ! (edges of latitude intervals for input data, in rad, in strictly
96    ! ascending order)
97
98    real, pointer:: plev(:)
99    ! pressure levels of input data, sorted in strictly ascending
100    ! order, converted to hPa
101
102    logical desc_lat ! latitude in descending order in the input file
103    logical desc_plev ! pressure levels in descending order in the input file
104
105    real, allocatable:: o3_in(:, :, :, :)
106    ! (n_lat, n_plev, n_month, read_climoz)
107    ! ozone climatologies from the input file
108    ! "o3_in(j, k, :, :)" is at latitude "latitude(j)" and pressure
109    ! level "plev(k)".
110    ! Third dimension is month index, first value may be December or January.
111    ! "o3_in(:, :, :, 1)" is for the day- night average, "o3_in(:, :, :, 2)"
112    ! is for daylight.
113
114    real missing_value
115
116    real, allocatable:: o3_regr_lat(:, :, :, :)
117    ! (jjm + 1, n_plev, 0:13, read_climoz)
118    ! mean of "o3_in" over a latitude interval of LMDZ
119    ! First dimension is latitude interval.
120    ! The latitude interval for "o3_regr_lat(j,:, :, :)" contains "rlatu(j)".
121    ! If "j" is between 2 and "jjm" then the interval is:
122    ! [rlatv(j), rlatv(j-1)]
123    ! If "j" is 1 or "jjm + 1" then the interval is:
124    ! [rlatv(1), pi / 2]
125    ! or:
126    ! [- pi / 2, rlatv(jjm)]
127    ! respectively.
128    ! "o3_regr_lat(:, k, :, :)" is for pressure level "plev(k)".
129    ! Third dimension is month number, 1 for January.
130    ! "o3_regr_lat(:, :, :, 1)" is average day and night,
131    ! "o3_regr_lat(:, :, :, 2)" is for daylight.
132
133    real, allocatable:: o3_out(:, :, :, :)
134    ! (jjm + 1, n_plev, 360, read_climoz)
135    ! regridded ozone climatology
136    ! "o3_out(j, k, l, :)" is at latitude "rlatu(j)", pressure
137    ! level "plev(k)" and date "January 1st 0h" + "tmidday(l)", in a
138    ! 360-day calendar.
139    ! "o3_out(:, :, :, 1)" is average day and night,
140    ! "o3_out(:, :, :, 2)" is for daylight.
141
142    integer j, k, l,m
143
144    ! For NetCDF:
145    integer ncid_in, ncid_out ! IDs for input and output files
146    integer varid_plev, varid_time, varid, ncerr, dimid
147    character(len=80) press_unit ! pressure unit
148
149    integer varid_in(read_climoz), varid_out(read_climoz)
150    ! index 1 is for average ozone day and night, index 2 is for
151    ! daylight ozone.
152
153    real, parameter:: tmidmonth(0:13) = (/(-15. + 30. * l, l = 0, 13)/)
154    ! (time to middle of month, in days since January 1st 0h, in a
155    ! 360-day calendar)
156    ! (We add values -15 and 375 so that, for example, day 3 of the year is
157    ! interpolated between the December and the January value.)
158
159    real, parameter:: tmidday(360) = (/(l + 0.5, l = 0, 359)/)
160    ! (time to middle of day, in days since January 1st 0h, in a
161    ! 360-day calendar)
162
163    !---------------------------------
164
165    print *, "Call sequence information: regr_lat_time_climoz"
166    call assert(read_climoz == 1 .or. read_climoz == 2, "regr_lat_time_climoz")
167
168    call nf95_open("climoz.nc", nf90_nowrite, ncid_in)
169
170    ! Get coordinates from the input file:
171
172    call nf95_inq_varid(ncid_in, "latitude", varid)
173    call nf95_gw_var(ncid_in, varid, latitude)
174    ! Convert from degrees to rad, because we will take the sine of latitude:
175    latitude = latitude / 180. * pi
176    n_lat = size(latitude)
177    ! We need to supply the latitudes to "regr1_step_av" in
178    ! ascending order, so invert order if necessary:
179    desc_lat = latitude(1) > latitude(n_lat)
180    if (desc_lat) latitude = latitude(n_lat:1:-1)
181
182    ! Compute edges of latitude intervals:
183    allocate(lat_in_edg(n_lat + 1))
184    lat_in_edg(1) = - pi / 2
185    forall (j = 2:n_lat) lat_in_edg(j) = (latitude(j - 1) + latitude(j)) / 2
186    lat_in_edg(n_lat + 1) = pi / 2
187    deallocate(latitude) ! pointer
188
189    call nf95_inq_varid(ncid_in, "plev", varid)
190    call nf95_gw_var(ncid_in, varid, plev)
191    n_plev = size(plev)
192    ! We only need the pressure coordinate to copy it to the output file.
193    ! The program "gcm" will assume that pressure levels are in
194    ! ascending order in the regridded climatology so invert order if
195    ! necessary:
196    desc_plev = plev(1) > plev(n_plev)
197    if (desc_plev) plev = plev(n_plev:1:-1)
198    call nf95_get_att(ncid_in, varid, "units", press_unit)
199    if (press_unit == "Pa") then
200       ! Convert to hPa:
201       plev = plev / 100.
202    elseif (press_unit /= "hPa") then
203       print *, "regr_lat_time_climoz: the only recognized units are Pa " &
204            // "and hPa."
205       stop 1
206    end if
207
208    ! Create the output file and get the variable IDs:
209    call prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
210         varid_time)
211
212    ! Write remaining coordinate variables:
213    call nf95_put_var(ncid_out, varid_plev, plev)
214    call nf95_put_var(ncid_out, varid_time, tmidday)
215
216    deallocate(plev) ! pointer
217
218    ! Get the  number of months:
219    call nf95_inq_dimid(ncid_in, "time", dimid)
220    call nf95_inquire_dimension(ncid_in, dimid, nclen=n_month)
221
222    allocate(o3_in(n_lat, n_plev, n_month, read_climoz))
223
224    call nf95_inq_varid(ncid_in, "tro3", varid_in(1))
225    ncerr = nf90_get_var(ncid_in, varid_in(1), o3_in(:, :, :, 1))
226    call handle_err("regr_lat_time_climoz nf90_get_var tro3", ncerr, ncid_in)
227
228    if (read_climoz == 2) then
229       call nf95_inq_varid(ncid_in, "tro3_daylight", varid_in(2))
230       ncerr = nf90_get_var(ncid_in, varid_in(2), o3_in(:, :, :, 2))
231       call handle_err("regr_lat_time_climoz nf90_get_var tro3_daylight", &
232            ncerr, ncid_in, varid_in(2))
233    end if
234
235    if (desc_lat) o3_in = o3_in(n_lat:1:-1, :, :, :)
236    if (desc_plev) o3_in = o3_in(:, n_plev:1:-1, :, :)
237
238    do m = 1, read_climoz
239       ncerr = nf90_get_att(ncid_in, varid_in(m), "missing_value", &
240            missing_value)
241       if (ncerr == nf90_noerr) then
242          do l = 1, n_month
243             ! Take care of latitudes where values are all missing:
244
245             ! Next to the south pole:
246             j = 1
247             do while (o3_in(j, 1, l, m) == missing_value)
248                j = j + 1
249             end do
250             if (j > 1) o3_in(:j-1, :, l, m) = &
251                  spread(o3_in(j, :, l, m), dim=1, ncopies=j-1)
252             
253             ! Next to the north pole:
254             j = n_lat
255             do while (o3_in(j, 1, l, m) == missing_value)
256                j = j - 1
257             end do
258             if (j < n_lat) o3_in(j+1:, :, l, m) = &
259                  spread(o3_in(j, :, l, m), dim=1, ncopies=n_lat-j)
260
261             ! Take care of missing values at high pressure:
262             do j = 1, n_lat
263                ! Find missing values, starting from top of atmosphere
264                ! and going down.
265                ! We have already taken care of latitudes full of
266                ! missing values so the highest level has a valid value.
267                k = 2
268                do while  (o3_in(j, k, l, m) /= missing_value .and. k < n_plev)
269                   k = k + 1
270                end do
271                ! Replace missing values with the valid value at the
272                ! lowest level above missing values:
273                if (o3_in(j, k, l, m) == missing_value) &
274                     o3_in(j, k:n_plev, l, m) = o3_in(j, k-1, l, m)
275             end do
276          end do
277       else
278          print *, "regr_lat_time_climoz: field ", m, &
279               ", no missing value attribute"
280       end if
281    end do
282
283    call nf95_close(ncid_in)
284
285    allocate(o3_regr_lat(nbp_lat, n_plev, 0:13, read_climoz))
286    allocate(o3_out(nbp_lat, n_plev, 360, read_climoz))
287
288    ! Regrid in latitude:
289    ! We average with respect to sine of latitude, which is
290    ! equivalent to weighting by cosine of latitude:
291    if (n_month == 12) then
292       print *, &
293            "Found 12 months in ozone climatologies, assuming periodicity..."
294       o3_regr_lat(nbp_lat:1:-1, :, 1:12, :) = regr1_step_av(o3_in, &
295            xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
296       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
297       ! in descending order)
298
299       ! Duplicate January and December values, in preparation of time
300       ! interpolation:
301       o3_regr_lat(:, :, 0, :) = o3_regr_lat(:, :, 12, :)
302       o3_regr_lat(:, :, 13, :) = o3_regr_lat(:, :, 1, :)
303    else
304       print *, "Using 14 months in ozone climatologies..."
305       o3_regr_lat(nbp_lat:1:-1, :, :, :) = regr1_step_av(o3_in, &
306            xs=sin(lat_in_edg), xt=sin((/- pi / 2, boundslat_reg(nbp_lat-1:1:-1,south), pi / 2/)))
307       ! (invert order of indices in "o3_regr_lat" because "rlatu" is
308       ! in descending order)
309    end if
310
311    ! Regrid in time by linear interpolation:
312    o3_out = regr3_lint(o3_regr_lat, tmidmonth, tmidday)
313
314    ! Write to file:
315    do m = 1, read_climoz
316       call nf95_put_var(ncid_out, varid_out(m), o3_out(nbp_lat:1:-1, :, :, m))
317       ! (The order of "rlatu" is inverted in the output file)
318    end do
319
320    call nf95_close(ncid_out)
321
322  end subroutine regr_lat_time_climoz
323
324  !********************************************
325
326  subroutine prepare_out(ncid_in, n_plev, ncid_out, varid_out, varid_plev, &
327       varid_time)
328
329    ! This subroutine creates the NetCDF output file, defines
330    ! dimensions and variables, and writes one of the coordinate variables.
331
332    use mod_grid_phy_lmdz, ONLY : nbp_lat
333    use netcdf95, only: nf95_create, nf95_def_dim, nf95_def_var, &
334         nf95_put_att, nf95_enddef, nf95_copy_att, nf95_put_var
335    use netcdf, only: nf90_clobber, nf90_float, nf90_global
336    use nrtype, only: pi
337    use regular_lonlat_mod, only : lat_reg
338
339    integer, intent(in):: ncid_in, n_plev
340    integer, intent(out):: ncid_out, varid_plev, varid_time
341
342    integer, intent(out):: varid_out(:) ! dim(1 or 2)
343    ! "varid_out(1)" is for average ozone day and night,
344    ! "varid_out(2)" is for daylight ozone.
345
346    ! Variables local to the procedure:
347
348    integer ncerr
349    integer dimid_rlatu, dimid_plev, dimid_time
350    integer varid_rlatu
351
352    !---------------------------
353
354    print *, "Call sequence information: prepare_out"
355
356    call nf95_create("climoz_LMDZ.nc", nf90_clobber, ncid_out)
357
358    ! Dimensions:
359    call nf95_def_dim(ncid_out, "time", 360, dimid_time)
360    call nf95_def_dim(ncid_out, "plev", n_plev, dimid_plev)
361    call nf95_def_dim(ncid_out, "rlatu", nbp_lat, dimid_rlatu)
362
363    ! Define coordinate variables:
364
365    call nf95_def_var(ncid_out, "time", nf90_float, dimid_time, varid_time)
366    call nf95_put_att(ncid_out, varid_time, "units", "days since 2000-1-1")
367    call nf95_put_att(ncid_out, varid_time, "calendar", "360_day")
368    call nf95_put_att(ncid_out, varid_time, "standard_name", "time")
369
370    call nf95_def_var(ncid_out, "plev", nf90_float, dimid_plev, varid_plev)
371    call nf95_put_att(ncid_out, varid_plev, "units", "millibar")
372    call nf95_put_att(ncid_out, varid_plev, "standard_name", "air_pressure")
373    call nf95_put_att(ncid_out, varid_plev, "long_name", "air pressure")
374
375    call nf95_def_var(ncid_out, "rlatu", nf90_float, dimid_rlatu, varid_rlatu)
376    call nf95_put_att(ncid_out, varid_rlatu, "units", "degrees_north")
377    call nf95_put_att(ncid_out, varid_rlatu, "standard_name", "latitude")
378
379    ! Define the primary variables:
380
381    call nf95_def_var(ncid_out, "tro3", nf90_float, &
382         (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(1))
383    call nf95_put_att(ncid_out, varid_out(1), "long_name", &
384         "ozone mole fraction")
385    call nf95_put_att(ncid_out, varid_out(1), "standard_name", &
386         "mole_fraction_of_ozone_in_air")
387
388    if (size(varid_out) == 2) then
389       call nf95_def_var(ncid_out, "tro3_daylight", nf90_float, &
390            (/dimid_rlatu, dimid_plev, dimid_time/), varid_out(2))
391       call nf95_put_att(ncid_out, varid_out(2), "long_name", &
392            "ozone mole fraction in daylight")
393    end if
394
395    ! Global attributes:
396
397    ! The following commands, copying attributes, may fail.
398    ! That is OK.
399    ! It should just mean that the attribute is not defined in the input file.
400
401    call nf95_copy_att(ncid_in, nf90_global, "Conventions", ncid_out, &
402         nf90_global, ncerr)
403    call handle_err_copy_att("Conventions")
404
405    call nf95_copy_att(ncid_in, nf90_global, "title", ncid_out, nf90_global, &
406         ncerr)
407    call handle_err_copy_att("title")
408
409    call nf95_copy_att(ncid_in, nf90_global, "institution", ncid_out, &
410         nf90_global, ncerr)
411    call handle_err_copy_att("institution")
412
413    call nf95_copy_att(ncid_in, nf90_global, "source", ncid_out, nf90_global, &
414         ncerr)
415    call handle_err_copy_att("source")
416
417    call nf95_put_att(ncid_out, nf90_global, "comment", "Regridded for LMDZ")
418
419    call nf95_enddef(ncid_out)
420
421    ! Write one of the coordinate variables:
422    call nf95_put_var(ncid_out, varid_rlatu, lat_reg(nbp_lat:1:-1) / pi * 180.)
423    ! (convert from rad to degrees and sort in ascending order)
424
425  contains
426
427    subroutine handle_err_copy_att(att_name)
428
429      use netcdf, only: nf90_noerr, nf90_strerror
430
431      character(len=*), intent(in):: att_name
432
433      !----------------------------------------
434
435      if (ncerr /= nf90_noerr) then
436         print *, "regr_lat_time_climoz_m prepare_out nf95_copy_att " &
437              // att_name // " -- " // trim(nf90_strerror(ncerr))
438      end if
439
440    end subroutine handle_err_copy_att
441
442  end subroutine prepare_out
443
444end module regr_lat_time_climoz_m
Note: See TracBrowser for help on using the repository browser.