source: LMDZ4/branches/LMDZ4-dev-20091210/libf/bibio/regr3_lint_m.F90 @ 5065

Last change on this file since 5065 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: 2.6 KB
Line 
1! $Id$
2module regr3_lint_m
3
4  ! Author: Lionel GUEZ
5
6  implicit none
7
8  interface regr3_lint
9     ! Each procedure regrids by linear interpolation.
10     ! The regridding operation is done on the third dimension of the
11     ! input array.
12     ! The difference betwwen the procedures is the rank of the first argument.
13     module procedure regr33_lint, regr34_lint
14  end interface
15
16  private
17  public regr3_lint
18
19contains
20
21  function regr33_lint(vs, xs, xt) result(vt)
22
23    ! "vs" has rank 3.
24
25    use assert_eq_m, only: assert_eq
26    use interpolation, only: hunt
27
28    real, intent(in):: vs(:, :, :)
29    ! (values of the function at source points "xs")
30
31    real, intent(in):: xs(:)
32    ! (abscissas of points in source grid, in strictly monotonic order)
33
34    real, intent(in):: xt(:)
35    ! (abscissas of points in target grid)
36
37    real vt(size(vs, 1), size(vs, 2), size(xt))
38    ! (values of the function on the target grid)
39
40    ! Variables local to the procedure:
41    integer is, it, ns
42    integer is_b ! "is" bound between 1 and "ns - 1"
43
44    !--------------------------------------
45
46    ns = assert_eq(size(vs, 3), size(xs), "regr33_lint ns")
47
48    is = -1 ! go immediately to bisection on first call to "hunt"
49
50    do it = 1, size(xt)
51       call hunt(xs, xt(it), is)
52       is_b = min(max(is, 1), ns - 1)
53       vt(:, :, it) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b) &
54            + (xt(it) - xs(is_b)) * vs(:, :, is_b+1)) / (xs(is_b+1) - xs(is_b))
55    end do
56
57  end function regr33_lint
58
59  !*********************************************************
60
61  function regr34_lint(vs, xs, xt) result(vt)
62
63    ! "vs" has rank 4.
64
65    use assert_eq_m, only: assert_eq
66    use interpolation, only: hunt
67
68    real, intent(in):: vs(:, :, :, :)
69    ! (values of the function at source points "xs")
70
71    real, intent(in):: xs(:)
72    ! (abscissas of points in source grid, in strictly monotonic order)
73
74    real, intent(in):: xt(:)
75    ! (abscissas of points in target grid)
76
77    real vt(size(vs, 1), size(vs, 2), size(xt), size(vs, 4))
78    ! (values of the function on the target grid)
79
80    ! Variables local to the procedure:
81    integer is, it, ns
82    integer is_b ! "is" bound between 1 and "ns - 1"
83
84    !--------------------------------------
85
86    ns = assert_eq(size(vs, 3), size(xs), "regr34_lint ns")
87
88    is = -1 ! go immediately to bisection on first call to "hunt"
89
90    do it = 1, size(xt)
91       call hunt(xs, xt(it), is)
92       is_b = min(max(is, 1), ns - 1)
93       vt(:, :, it, :) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b, :) &
94            + (xt(it) - xs(is_b)) * vs(:, :, is_b+1, :)) &
95            / (xs(is_b+1) - xs(is_b))
96    end do
97
98  end function regr34_lint
99
100end module regr3_lint_m
Note: See TracBrowser for help on using the repository browser.