source: LMDZ6/branches/LMDZ_ECRad/libf/obsolete/regr1_conserv_m.F90 @ 5098

Last change on this file since 5098 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.
File size: 10.6 KB
Line 
1module regr1_conserv_m
2
3  ! Author: Lionel GUEZ
4
5  use assert_eq_m, only: assert_eq
6  use assert_m, only: assert
7  use interpolation, only: locate
8
9  implicit none
10
11  interface regr1_conserv
12
13     ! This generic procedure regrids a piecewise linear function (not
14     ! necessarily continuous) by averaging it. This is a conservative
15     ! regridding. The regridding operation is done on the first
16     ! dimension of the input array. Input are positions of cell
17     ! edges. The target grid should be included in the source grid:
18     ! no extrapolation is allowed.
19
20     ! The only difference between the procedures is the rank of the
21     ! first argument.
22
23     ! real, intent(in), rank >= 1:: vs ! (ns, ...)
24     ! averages of cells of the source grid
25     ! vs(is, ...) for [xs(is), xs(is + 1)]
26
27     ! real, intent(in):: xs(:) ! (ns + 1)
28     ! edges of cells of the source grid, in strictly ascending order
29
30     ! real, intent(in):: xt(:) ! (nt + 1)
31     ! edges of cells of the target grid, in strictly ascending order
32
33     ! real, intent(in), optional, rank >= 1:: slope ! (ns, ...)
34     ! same rank as vs
35     ! slopes inside cells of the source grid
36     ! slope(is, ...) for [xs(is), xs(is + 1)]
37     ! If not present, slopes are taken equal to 0. The regridding is
38     ! then first order.
39
40     ! real, intent(out), rank >= 1:: vt(nt, ...)
41     ! has the same rank as vs and slope
42     ! averages of cells of the target grid
43     ! vt(it, ...) for  [xt(it), xt(it + 1)]
44
45     ! ns and nt must be >= 1.
46
47     ! See notes for explanations on the algorithm and justification
48     ! of algorithmic choices.
49
50     module procedure regr11_conserv, regr12_conserv, regr13_conserv, &
51          regr14_conserv
52  end interface regr1_conserv
53
54  private
55  public regr1_conserv
56
57contains
58
59  subroutine regr11_conserv(vs, xs, xt, vt, slope)
60
61    ! vs and slope have rank 1.
62
63    real, intent(in):: vs(:)
64    real, intent(in):: xs(:)
65    real, intent(in):: xt(:)
66    real, intent(out):: vt(:)
67    real, intent(in), optional:: slope(:)
68
69    ! Local:
70    integer is, it, ns, nt
71    logical slope_present
72
73    !---------------------------------------------
74
75    ns = assert_eq(size(vs), size(xs) - 1, "regr11_conserv ns")
76    nt = assert_eq(size(xt) - 1, size(vt), "regr11_conserv nt")
77
78    ! Quick check on sort order:
79    call assert(xs(1) < xs(2), "regr11_conserv xs bad order")
80    call assert(xt(1) < xt(2), "regr11_conserv xt bad order")
81
82    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
83         "regr11_conserv extrapolation")
84    slope_present = present(slope)
85
86    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
87    do it = 1, nt
88       ! 1 <= is <= ns
89       ! xs(is) <= xt(it) < xs(is + 1)
90       if (xt(it + 1) <= xs(is + 1)) then
91          vt(it) = mean_lin(xt(it), xt(it + 1))
92       else
93          vt(it) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
94          is = is + 1
95          do while (xs(is + 1) < xt(it + 1))
96             ! 1 <= is <= ns - 1
97             vt(it) = vt(it) + (xs(is + 1) - xs(is)) * vs(is)
98             is = is + 1
99          end do
100          ! 1 <= is <= ns
101          vt(it) = (vt(it) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) &
102               - xs(is))) / (xt(it + 1) - xt(it))
103       end if
104
105       if (xs(is + 1) == xt(it + 1)) is = is + 1
106       ! 1 <= is <= ns .or. it == nt
107    end do
108
109  contains
110
111    real function mean_lin(a, b)
112
113      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
114
115      real, intent(in):: a, b
116
117      !---------------------------------------------
118
119      if (slope_present) then
120         mean_lin = slope(is) / 2. * (a + b - xs(is) - xs(is + 1)) + vs(is)
121      else
122         mean_lin = vs(is)
123      end if
124
125    end function mean_lin
126
127  end subroutine regr11_conserv
128
129  !********************************************
130
131  subroutine regr12_conserv(vs, xs, xt, vt, slope)
132
133    ! vs and slope have rank 2.
134
135    real, intent(in):: vs(:, :)
136    real, intent(in):: xs(:)
137    real, intent(in):: xt(:)
138    real, intent(out):: vt(:, :)
139    real, intent(in), optional:: slope(:, :)
140
141    ! Local:
142    integer is, it, ns, nt, n2
143    logical slope_present
144
145    !---------------------------------------------
146
147    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_conserv ns")
148    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr12_conserv nt")
149    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr12_conserv n2")
150
151    ! Quick check on sort order:
152    call assert(xs(1) < xs(2), "regr12_conserv xs bad order")
153    call assert(xt(1) < xt(2), "regr12_conserv xt bad order")
154
155    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
156         "regr12_conserv extrapolation")
157    slope_present = present(slope)
158
159    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
160    do it = 1, nt
161       ! 1 <= is <= ns
162       ! xs(is) <= xt(it) < xs(is + 1)
163       if (xt(it + 1) <= xs(is + 1)) then
164          vt(it, :) = mean_lin(xt(it), xt(it + 1))
165       else
166          vt(it, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
167          is = is + 1
168          do while (xs(is + 1) < xt(it + 1))
169             ! 1 <= is <= ns - 1
170             vt(it, :) = vt(it, :) + (xs(is + 1) - xs(is)) * vs(is, :)
171             is = is + 1
172          end do
173          ! 1 <= is <= ns
174          vt(it, :) = (vt(it, :) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) &
175               - xs(is))) / (xt(it + 1) - xt(it))
176       end if
177
178       if (xs(is + 1) == xt(it + 1)) is = is + 1
179       ! 1 <= is <= ns .or. it == nt
180    end do
181
182  contains
183
184    function mean_lin(a, b)
185
186      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
187
188      real, intent(in):: a, b
189      real mean_lin(n2)
190
191      !---------------------------------------------
192
193      if (slope_present) then
194         mean_lin = slope(is, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
195              + vs(is, :)
196      else
197         mean_lin = vs(is, :)
198      end if
199
200    end function mean_lin
201
202  end subroutine regr12_conserv
203
204  !********************************************
205
206  subroutine regr13_conserv(vs, xs, xt, vt, slope)
207
208    ! vs and slope have rank 3.
209
210    real, intent(in):: vs(:, :, :)
211    real, intent(in):: xs(:)
212    real, intent(in):: xt(:)
213    real, intent(out):: vt(:, :, :)
214    real, intent(in), optional:: slope(:, :, :)
215
216    ! Local:
217    integer is, it, ns, nt, n2, n3
218    logical slope_present
219
220    !---------------------------------------------
221
222    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_conserv ns")
223    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr13_conserv nt")
224    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr13_conserv n2")
225    n3 = assert_eq(size(vs, 3), size(vt, 3), "regr13_conserv n3")
226
227    ! Quick check on sort order:
228    call assert(xs(1) < xs(2), "regr13_conserv xs bad order")
229    call assert(xt(1) < xt(2), "regr13_conserv xt bad order")
230
231    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
232         "regr13_conserv extrapolation")
233    slope_present = present(slope)
234
235    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
236    do it = 1, nt
237       ! 1 <= is <= ns
238       ! xs(is) <= xt(it) < xs(is + 1)
239       if (xt(it + 1) <= xs(is + 1)) then
240          vt(it, :, :) = mean_lin(xt(it), xt(it + 1))
241       else
242          vt(it, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
243          is = is + 1
244          do while (xs(is + 1) < xt(it + 1))
245             ! 1 <= is <= ns - 1
246             vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - xs(is)) * vs(is, :, :)
247             is = is + 1
248          end do
249          ! 1 <= is <= ns
250          vt(it, :, :) = (vt(it, :, :) + mean_lin(xs(is), xt(it + 1)) &
251               * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it))
252       end if
253
254       if (xs(is + 1) == xt(it + 1)) is = is + 1
255       ! 1 <= is <= ns .or. it == nt
256    end do
257
258  contains
259
260    function mean_lin(a, b)
261
262      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
263
264      real, intent(in):: a, b
265      real mean_lin(n2, n3)
266
267      !---------------------------------------------
268
269      if (slope_present) then
270         mean_lin = slope(is, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
271              + vs(is, :, :)
272      else
273         mean_lin = vs(is, :, :)
274      end if
275
276    end function mean_lin
277
278  end subroutine regr13_conserv
279
280  !********************************************
281
282  subroutine regr14_conserv(vs, xs, xt, vt, slope)
283
284    ! vs and slope have rank 4.
285
286    real, intent(in):: vs(:, :, :, :)
287    real, intent(in):: xs(:)
288    real, intent(in):: xt(:)
289    real, intent(out):: vt(:, :, :, :)
290    real, intent(in), optional:: slope(:, :, :, :)
291
292    ! Local:
293    integer is, it, ns, nt, n2, n3, n4
294    logical slope_present
295
296    !---------------------------------------------
297
298    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_conserv ns")
299    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr14_conserv nt")
300    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr14_conserv n2")
301    n3 = assert_eq(size(vs, 3), size(vt, 3), "regr14_conserv n3")
302    n4 = assert_eq(size(vs, 4), size(vt, 4), "regr14_conserv n4")
303
304    ! Quick check on sort order:
305    call assert(xs(1) < xs(2), "regr14_conserv xs bad order")
306    call assert(xt(1) < xt(2), "regr14_conserv xt bad order")
307
308    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
309         "regr14_conserv extrapolation")
310    slope_present = present(slope)
311
312    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
313    do it = 1, nt
314       ! 1 <= is <= ns
315       ! xs(is) <= xt(it) < xs(is + 1)
316       if (xt(it + 1) <= xs(is + 1)) then
317          vt(it, :, :, :) = mean_lin(xt(it), xt(it + 1))
318       else
319          vt(it, :, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
320          is = is + 1
321          do while (xs(is + 1) < xt(it + 1))
322             ! 1 <= is <= ns - 1
323             vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - xs(is)) &
324                  * vs(is, :, :, :)
325             is = is + 1
326          end do
327          ! 1 <= is <= ns
328          vt(it, :, :, :) = (vt(it, :, :, :) + mean_lin(xs(is), xt(it + 1)) &
329               * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it))
330       end if
331
332       if (xs(is + 1) == xt(it + 1)) is = is + 1
333       ! 1 <= is <= ns .or. it == nt
334    end do
335
336  contains
337
338    function mean_lin(a, b)
339
340      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
341
342      real, intent(in):: a, b
343      real mean_lin(n2, n3, n4)
344
345      !---------------------------------------------
346
347      if (slope_present) then
348         mean_lin = slope(is, :, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
349              + vs(is, :, :, :)
350      else
351         mean_lin = vs(is, :, :, :)
352      end if
353
354    end function mean_lin
355
356  end subroutine regr14_conserv
357
358end module regr1_conserv_m
Note: See TracBrowser for help on using the repository browser.