source: LMDZ4/branches/LMDZ4-dev-20091210/libf/bibio/regr1_step_av_m.F90 @ 5074

Last change on this file since 5074 was 1263, checked in by lguez, 15 years ago

1) Reactivated ability to read ozone (that was deactivated because of
dependency on version of IOIPSL). Added ability to read a pressure
coordinate in Pa in "regr_lat_time_climoz".

2) Added the ability to read a second ozone climatology, corresponding to
daylight ozone:

-- "read_climoz" is now an integer variable, instead of a logical
variable.

-- Added argument "read_climoz" to "phys_state_var_init",
"phys_output_open" and "regr_lat_time_climoz".

-- Created new variable "ozone_daylight" for "hist*.nc" output files.

-- Added a third dimension to variable "wo" in module
"phys_state_var_mod" and variable "POZON" in "radlwsw": index 1 for
average day-night ozone, index 2 for daylight ozone.

-- Added a fourth dimension to variables "o3_in", "o3_regr_lat" and
"o3_out" in "regr_lat_time_climoz": index 1 for average day-night
ozone, index 2 for daylight ozone.

-- In "physiq", moved call to "conf_phys" before call to
"phys_state_var_init". Thus, "conf_phys" is now inside the block "if
(first)" instead of "IF (debut)". There were definitions of "bl95_b0"
and "bl95_b1" that were useless because the variables were overwritten
by "conf_phys". Removed those definitions.

-- In "radlwsw", we pass the average day-night ozone to "LW_LMDAR4"
and the daylight ozone, if we have it, to "SW_LMDAR4" or
"SW_AEROAR4". If we do not have a specific field for daylight ozone
then "SW_LMDAR4" or "SW_AEROAR4" just get the average day-night ozone.

-- "regr_lat_time_climoz" now manages latitudes where the input ozone
field is missing at all levels (polar night).

-- Encapsulated "radlwsw" in a module.

3) Modifications to make sequential and parallel versions of
"create_etat0_limit" almost identical:

-- In "dyn3dpar/create_etat0_limit.F". No need to call
"phys_state_var_init", removed "use phys_state_var_mod" statement. No
need for "clesphys.h", removed "include" statement.

-- In "dyn3dpar/etat0_netcdf.F". Added argument "tau_ratqs" in call to
"conf_phys" (this bug was already corrected in "dyn3d"). Moved call to
"inifilr" after call to "infotrac_init" (as in "dyn3d").

4) Other peripheral modifications:

-- Added procedures "nf95_get_att" and "nf95_def_var_scalar" in
NetCDF95 interface. Overloaded "nf95_put_var" with three more
procedures: "nf95_put_var_FourByteReal", "nf95_put_var_FourByteInt",
"nf95_put_var_1D_FourByteInt".

-- Overloaded "regr1_step_av" with one more procedure:
"regr14_step_av". Overloaded "regr3_lint" with one more procedure:
"regr34_lint".

-- Corrected call to "Init_Phys_lmdz" in "dyn3d/create_etat0_limit.F":
the last argument should be an array, not a scalar.

-- Encapsulated "conf_phys" in a module.

-- Splitted module "regr_pr" into "regr_pr_av_m" and "regr_pr_int_m".

5) Tests:

This revision was compared to revision 1259, with optimization options
"debug" and "dev", parallelization options "none", "mpi", "omp" and
"mpi_omp", 1 and 2 MPI processes, 1 and 2 OpenMP threads, with the
compiler "FORTRAN90/SX Version 2.0 for SX-8". Both programs
"create_etat0_limit" and "gcm" were tested. In all cases,
parallelization does not change the results. With "read_climoz = 0" in
the ".def" files, the results of revision 1259 and of this revision
are the same.
File size: 8.1 KB
Line 
1! $Id$
2module regr1_step_av_m
3
4  ! Author: Lionel GUEZ
5
6  implicit none
7
8  interface regr1_step_av
9
10     ! Each procedure regrids a step function by averaging it.
11     ! The regridding operation is done on the first dimension of the
12     ! input array.
13     ! Source grid contains edges of steps.
14     ! Target grid contains positions of cell edges.
15     ! The target grid should be included in the source grid: no
16     ! extrapolation is allowed.
17     ! The difference between the procedures is the rank of the first argument.
18
19     module procedure regr11_step_av, regr12_step_av, regr13_step_av, &
20          regr14_step_av
21  end interface
22
23  private
24  public regr1_step_av
25
26contains
27
28  function regr11_step_av(vs, xs, xt) result(vt)
29
30    ! "vs" has rank 1.
31
32    use assert_eq_m, only: assert_eq
33    use assert_m, only: assert
34    use interpolation, only: locate
35
36    real, intent(in):: vs(:) ! values of steps on the source grid
37    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
38
39    real, intent(in):: xs(:)
40    ! (edges of of steps on the source grid, in strictly increasing order)
41
42    real, intent(in):: xt(:)
43    ! (edges of cells of the target grid, in strictly increasing order)
44
45    real vt(size(xt) - 1) ! average values on the target grid
46    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
47
48    ! Variables local to the procedure:
49    integer is, it, ns, nt
50    real left_edge
51
52    !---------------------------------------------
53
54    ns = assert_eq(size(vs), size(xs) - 1, "regr11_step_av ns")
55    nt = size(xt) - 1
56    ! Quick check on sort order:
57    call assert(xs(1) < xs(2), "regr11_step_av xs bad order")
58    call assert(xt(1) < xt(2), "regr11_step_av xt bad order")
59
60    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
61         "regr11_step_av extrapolation")
62
63    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
64    do it = 1, nt
65       ! 1 <= is <= ns
66       ! xs(is) <= xt(it) < xs(is + 1)
67       ! Compute "vt(it)":
68       left_edge = xt(it)
69       vt(it) = 0.
70       do while (xs(is + 1) < xt(it + 1))
71          ! 1 <= is <= ns - 1
72          vt(it) = vt(it) + (xs(is + 1) - left_edge) * vs(is)
73          is = is + 1
74          left_edge = xs(is)
75       end do
76       ! 1 <= is <= ns
77       vt(it) = (vt(it) + (xt(it + 1) - left_edge) * vs(is)) &
78            / (xt(it + 1) - xt(it))
79       if (xs(is + 1) == xt(it + 1)) is = is + 1
80       ! 1 <= is <= ns .or. it == nt
81    end do
82
83  end function regr11_step_av
84
85  !********************************************
86
87  function regr12_step_av(vs, xs, xt) result(vt)
88
89    ! "vs" has rank 2.
90
91    use assert_eq_m, only: assert_eq
92    use assert_m, only: assert
93    use interpolation, only: locate
94
95    real, intent(in):: vs(:, :) ! values of steps on the source grid
96    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
97
98    real, intent(in):: xs(:)
99    ! (edges of steps on the source grid, in strictly increasing order)
100
101    real, intent(in):: xt(:)
102    ! (edges of cells of the target grid, in strictly increasing order)
103
104    real vt(size(xt) - 1, size(vs, 2)) ! average values on the target grid
105    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
106
107    ! Variables local to the procedure:
108    integer is, it, ns, nt
109    real left_edge
110
111    !---------------------------------------------
112
113    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_step_av ns")
114    nt = size(xt) - 1
115
116    ! Quick check on sort order:
117    call assert(xs(1) < xs(2), "regr12_step_av xs bad order")
118    call assert(xt(1) < xt(2), "regr12_step_av xt bad order")
119
120    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
121         "regr12_step_av extrapolation")
122
123    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
124    do it = 1, nt
125       ! 1 <= is <= ns
126       ! xs(is) <= xt(it) < xs(is + 1)
127       ! Compute "vt(it, :)":
128       left_edge = xt(it)
129       vt(it, :) = 0.
130       do while (xs(is + 1) < xt(it + 1))
131          ! 1 <= is <= ns - 1
132          vt(it, :) = vt(it, :) + (xs(is + 1) - left_edge) * vs(is, :)
133          is = is + 1
134          left_edge = xs(is)
135       end do
136       ! 1 <= is <= ns
137       vt(it, :) = (vt(it, :) + (xt(it + 1) - left_edge) * vs(is, :)) &
138            / (xt(it + 1) - xt(it))
139       if (xs(is + 1) == xt(it + 1)) is = is + 1
140       ! 1 <= is <= ns .or. it == nt
141    end do
142
143  end function regr12_step_av
144
145  !********************************************
146
147  function regr13_step_av(vs, xs, xt) result(vt)
148
149    ! "vs" has rank 3.
150
151    use assert_eq_m, only: assert_eq
152    use assert_m, only: assert
153    use interpolation, only: locate
154
155    real, intent(in):: vs(:, :, :) ! values of steps on the source grid
156    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
157
158    real, intent(in):: xs(:)
159    ! (edges of steps on the source grid, in strictly increasing order)
160
161    real, intent(in):: xt(:)
162    ! (edges of cells of the target grid, in strictly increasing order)
163
164    real vt(size(xt) - 1, size(vs, 2), size(vs, 3))
165    ! (average values on the target grid)
166    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
167
168    ! Variables local to the procedure:
169    integer is, it, ns, nt
170    real left_edge
171
172    !---------------------------------------------
173
174    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_step_av ns")
175    nt = size(xt) - 1
176
177    ! Quick check on sort order:
178    call assert(xs(1) < xs(2), "regr13_step_av xs bad order")
179    call assert(xt(1) < xt(2), "regr13_step_av xt bad order")
180
181    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
182         "regr13_step_av extrapolation")
183
184    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
185    do it = 1, nt
186       ! 1 <= is <= ns
187       ! xs(is) <= xt(it) < xs(is + 1)
188       ! Compute "vt(it, :, :)":
189       left_edge = xt(it)
190       vt(it, :, :) = 0.
191       do while (xs(is + 1) < xt(it + 1))
192          ! 1 <= is <= ns - 1
193          vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - left_edge) * vs(is, :, :)
194          is = is + 1
195          left_edge = xs(is)
196       end do
197       ! 1 <= is <= ns
198       vt(it, :, :) = (vt(it, :, :) &
199            + (xt(it + 1) - left_edge) * vs(is, :, :)) / (xt(it + 1) - xt(it))
200       if (xs(is + 1) == xt(it + 1)) is = is + 1
201       ! 1 <= is <= ns .or. it == nt
202    end do
203
204  end function regr13_step_av
205
206  !********************************************
207
208  function regr14_step_av(vs, xs, xt) result(vt)
209
210    ! "vs" has rank 4.
211
212    use assert_eq_m, only: assert_eq
213    use assert_m, only: assert
214    use interpolation, only: locate
215
216    real, intent(in):: vs(:, :, :, :) ! values of steps on the source grid
217    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)
218
219    real, intent(in):: xs(:)
220    ! (edges of steps on the source grid, in strictly increasing order)
221
222    real, intent(in):: xt(:)
223    ! (edges of cells of the target grid, in strictly increasing order)
224
225    real vt(size(xt) - 1, size(vs, 2), size(vs, 3), size(vs, 4))
226    ! (average values on the target grid)
227    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
228
229    ! Variables local to the procedure:
230    integer is, it, ns, nt
231    real left_edge
232
233    !---------------------------------------------
234
235    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_step_av ns")
236    nt = size(xt) - 1
237
238    ! Quick check on sort order:
239    call assert(xs(1) < xs(2), "regr14_step_av xs bad order")
240    call assert(xt(1) < xt(2), "regr14_step_av xt bad order")
241
242    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
243         "regr14_step_av extrapolation")
244
245    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
246    do it = 1, nt
247       ! 1 <= is <= ns
248       ! xs(is) <= xt(it) < xs(is + 1)
249       ! Compute "vt(it, :, :, :)":
250       left_edge = xt(it)
251       vt(it, :, :, :) = 0.
252       do while (xs(is + 1) < xt(it + 1))
253          ! 1 <= is <= ns - 1
254          vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - left_edge) &
255               * vs(is, :, :, :)
256          is = is + 1
257          left_edge = xs(is)
258       end do
259       ! 1 <= is <= ns
260       vt(it, :, :, :) = (vt(it, :, :, :) + (xt(it + 1) - left_edge) &
261            * vs(is, :, :, :)) / (xt(it + 1) - xt(it))
262       if (xs(is + 1) == xt(it + 1)) is = is + 1
263       ! 1 <= is <= ns .or. it == nt
264    end do
265
266  end function regr14_step_av
267
268end module regr1_step_av_m
Note: See TracBrowser for help on using the repository browser.