source: LMDZ5/branches/IPSLCM6.0.11.rc1/libf/misc/regr_conserv_m.F90 @ 5292

Last change on this file since 5292 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: 15.4 KB
Line 
1MODULE regr_conserv_m
2
3  USE assert_eq_m,   ONLY: assert_eq
4  USE assert_m,      ONLY: assert
5  USE interpolation, ONLY: locate
6
7  IMPLICIT NONE
8
9! Purpose: Each procedure regrids a piecewise linear function (not necessarily
10!          continuous) by averaging it. This is a conservative regridding.
11!          The regridding operation is done along dimension "ix" of the input
12!          field "vs". Input are positions of cell edges.
13! Remarks:
14!   * The target grid should be included in the source grid:
15!     no extrapolation is allowed.
16!   * Field on target grid "vt" has same rank, slopes and averages as "vs".
17!   * If optional argument "slope" is not given, 0 is assumed for slopes.
18!     Then the regridding is first order instead of second.
19!   * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt)
20
21! Argument                         Type         Description
22!-------------------------------------------------------------------------------
23!  INTEGER, INTENT(IN)  :: ix      Scalar       dimension regridded <=rank(vs)
24!  REAL,    INTENT(IN)  :: vs(*)   Rank>=1      averages in source grid cells
25!  REAL,    INTENT(IN)  :: xs(:)   Vector(ns+1) edges of source grid, asc. order
26!  REAL,    INTENT(IN)  :: xt(:)   Vector(nt+1) edges of target grid, asc. order
27!  REAL,    INTENT(OUT) :: vt(*)   Rank>=1      regridded field
28! [REAL,    INTENT(IN)  :: slope]  Rank>=1      slopes in source grid cells
29
30  INTERFACE regr_conserv
31    ! The procedures differ only from the rank of the input/output fields.
32    MODULE PROCEDURE regr1_conserv, regr2_conserv, regr3_conserv, &
33                     regr4_conserv, regr5_conserv
34  END INTERFACE
35
36  PRIVATE
37  PUBLIC :: regr_conserv
38
39CONTAINS
40
41!-------------------------------------------------------------------------------
42!
43SUBROUTINE regr1_conserv(ix, vs, xs, xt, vt, slope)
44!
45!-------------------------------------------------------------------------------
46! Arguments:
47  INTEGER,        INTENT(IN)  :: ix
48  REAL,           INTENT(IN)  :: vs(:)
49  REAL,           INTENT(IN)  :: xs(:), xt(:)
50  REAL,           INTENT(OUT) :: vt(:)
51  REAL, OPTIONAL, INTENT(IN)  :: slope(:)
52!-------------------------------------------------------------------------------
53! Local variables:
54  INTEGER :: is, it, ns, nt
55  REAL    :: co, idt
56  LOGICAL :: lslope
57!-------------------------------------------------------------------------------
58  lslope=PRESENT(slope)
59  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
60  is=locate(xs,xt(1))                    !--- 1<= is <= ns (no extrapolation)
61  vt(:)=0.
62  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
63    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
64      CALL mean_lin(xt(it),xt(it+1))
65    ELSE
66      CALL mean_lin(xt(it),xs(is+1)); is=is+1
67      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
68      CALL mean_lin(xs(is),xs(is+1)); is=is+1
69      END DO
70      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
71    END IF
72    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
73  END DO
74
75CONTAINS
76
77!-------------------------------------------------------------------------------
78SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
79!-------------------------------------------------------------------------------
80! Arguments:
81  REAL, INTENT(IN)    :: a, b
82!-------------------------------------------------------------------------------
83  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
84    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
85    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
86    vt(it) = vt(it)+idt*(b-a)*(vs(is)+co*slope(is))
87  ELSE
88    vt(it) = vt(it)+idt*(b-a)* vs(is)
89  END IF
90
91END SUBROUTINE mean_lin
92!-------------------------------------------------------------------------------
93
94END SUBROUTINE regr1_conserv
95!
96!-------------------------------------------------------------------------------
97
98
99!-------------------------------------------------------------------------------
100!
101SUBROUTINE regr2_conserv(ix, vs, xs, xt, vt, slope)
102!
103!-------------------------------------------------------------------------------
104! Arguments:
105  INTEGER,        INTENT(IN)  :: ix
106  REAL,           INTENT(IN)  :: vs(:,:)
107  REAL,           INTENT(IN)  :: xs(:), xt(:)
108  REAL,           INTENT(OUT) :: vt(:,:)
109  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:)
110!-------------------------------------------------------------------------------
111! Local variables:
112  INTEGER :: is, it, ns, nt
113  REAL    :: co, idt
114  LOGICAL :: lslope
115!-------------------------------------------------------------------------------
116  lslope=PRESENT(slope)
117  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
118  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
119  vt(:,:)=0.
120  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
121    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
122      CALL mean_lin(xt(it),xt(it+1))
123    ELSE
124      CALL mean_lin(xt(it),xs(is+1)); is=is+1
125      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
126      CALL mean_lin(xs(is),xs(is+1)); is=is+1
127      END DO
128      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
129    END IF
130    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
131  END DO
132
133CONTAINS
134
135!-------------------------------------------------------------------------------
136SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
137!-------------------------------------------------------------------------------
138! Arguments:
139  REAL, INTENT(IN)    :: a, b
140!-------------------------------------------------------------------------------
141  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
142    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
143    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
144    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)*(vs(is,:)+co*slope(is,:))
145    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)*(vs(:,is)+co*slope(:,is))
146  ELSE
147    IF(ix==1) vt(it,:) = vt(it,:)+idt*(b-a)* vs(is,:)
148    IF(ix==2) vt(:,it) = vt(:,it)+idt*(b-a)* vs(:,is)
149  END IF
150
151END SUBROUTINE mean_lin
152!-------------------------------------------------------------------------------
153
154END SUBROUTINE regr2_conserv
155!
156!-------------------------------------------------------------------------------
157
158
159!-------------------------------------------------------------------------------
160!
161SUBROUTINE regr3_conserv(ix, vs, xs, xt, vt, slope)
162!
163!-------------------------------------------------------------------------------
164! Arguments:
165  INTEGER,        INTENT(IN)  :: ix
166  REAL,           INTENT(IN)  :: vs(:,:,:)
167  REAL,           INTENT(IN)  :: xs(:), xt(:)
168  REAL,           INTENT(OUT) :: vt(:,:,:)
169  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:)
170!-------------------------------------------------------------------------------
171! Local variables:
172  INTEGER :: is, it, ns, nt
173  REAL    :: co, idt
174  LOGICAL :: lslope
175!-------------------------------------------------------------------------------
176  lslope=PRESENT(slope)
177  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
178  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
179  vt(:,:,:)=0.
180  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
181    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
182      CALL mean_lin(xt(it),xt(it+1))
183    ELSE
184      CALL mean_lin(xt(it),xs(is+1)); is=is+1
185      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
186      CALL mean_lin(xs(is),xs(is+1)); is=is+1
187      END DO
188      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
189    END IF
190    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
191  END DO
192
193CONTAINS
194
195!-------------------------------------------------------------------------------
196SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
197!-------------------------------------------------------------------------------
198! Arguments:
199  REAL, INTENT(IN)    :: a, b
200!-------------------------------------------------------------------------------
201  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
202    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
203    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
204    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)*(vs(is,:,:)+co*slope(is,:,:))
205    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)*(vs(:,is,:)+co*slope(:,is,:))
206    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)*(vs(:,:,is)+co*slope(:,:,is))
207  ELSE
208    IF(ix==1) vt(it,:,:) = vt(it,:,:)+idt*(b-a)* vs(is,:,:)
209    IF(ix==2) vt(:,it,:) = vt(:,it,:)+idt*(b-a)* vs(:,is,:)
210    IF(ix==3) vt(:,:,it) = vt(:,:,it)+idt*(b-a)* vs(:,:,is)
211  END IF
212
213END SUBROUTINE mean_lin
214!-------------------------------------------------------------------------------
215
216END SUBROUTINE regr3_conserv
217!
218!-------------------------------------------------------------------------------
219
220
221!-------------------------------------------------------------------------------
222!
223SUBROUTINE regr4_conserv(ix, vs, xs, xt, vt, slope)
224!
225!-------------------------------------------------------------------------------
226! Arguments:
227  INTEGER,        INTENT(IN)  :: ix
228  REAL,           INTENT(IN)  :: vs(:,:,:,:)
229  REAL,           INTENT(IN)  :: xs(:), xt(:)
230  REAL,           INTENT(OUT) :: vt(:,:,:,:)
231  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:)
232!-------------------------------------------------------------------------------
233! Local variables:
234  INTEGER :: is, it, ns, nt
235  REAL    :: co, idt
236  LOGICAL :: lslope
237!-------------------------------------------------------------------------------
238  lslope=PRESENT(slope)
239  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
240  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
241  vt(:,:,:,:)=0.
242  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
243    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
244      CALL mean_lin(xt(it),xt(it+1))
245    ELSE
246      CALL mean_lin(xt(it),xs(is+1)); is=is+1
247      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
248      CALL mean_lin(xs(is),xs(is+1)); is=is+1
249      END DO
250      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
251    END IF
252    IF(xs(is+1)==xt(it+1)) is=is+1       !--- 1<=is<=ns .OR. it==nt
253  END DO
254
255CONTAINS
256
257!-------------------------------------------------------------------------------
258SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
259!-------------------------------------------------------------------------------
260! Arguments:
261  REAL, INTENT(IN) :: a, b
262!-------------------------------------------------------------------------------
263  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
264    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
265    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
266    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)*(vs(is,:,:,:)+co*slope(is,:,:,:))
267    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)*(vs(:,is,:,:)+co*slope(:,is,:,:))
268    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)*(vs(:,:,is,:)+co*slope(:,:,is,:))
269    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)*(vs(:,:,:,is)+co*slope(:,:,:,is))
270  ELSE
271    IF(ix==1) vt(it,:,:,:) = vt(it,:,:,:)+idt*(b-a)* vs(is,:,:,:)
272    IF(ix==2) vt(:,it,:,:) = vt(:,it,:,:)+idt*(b-a)* vs(:,is,:,:)
273    IF(ix==3) vt(:,:,it,:) = vt(:,:,it,:)+idt*(b-a)* vs(:,:,is,:)
274    IF(ix==4) vt(:,:,:,it) = vt(:,:,:,it)+idt*(b-a)* vs(:,:,:,is)
275  END IF
276
277END SUBROUTINE mean_lin
278!-------------------------------------------------------------------------------
279
280END SUBROUTINE regr4_conserv
281!
282!-------------------------------------------------------------------------------
283
284
285!-------------------------------------------------------------------------------
286!
287SUBROUTINE regr5_conserv(ix, vs, xs, xt, vt, slope)
288!
289!-------------------------------------------------------------------------------
290! Arguments:
291  INTEGER,        INTENT(IN)  :: ix
292  REAL,           INTENT(IN)  :: vs(:,:,:,:,:)
293  REAL,           INTENT(IN)  :: xs(:), xt(:)
294  REAL,           INTENT(OUT) :: vt(:,:,:,:,:)
295  REAL, OPTIONAL, INTENT(IN)  :: slope(:,:,:,:,:)
296!-------------------------------------------------------------------------------
297! Local variables:
298  INTEGER :: is, it, ns, nt
299  REAL    :: co, idt
300  LOGICAL :: lslope
301!-------------------------------------------------------------------------------
302  lslope=PRESENT(slope)
303  CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt)
304  is=locate(xs,xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
305  vt(:,:,:,:,:)=0.
306  DO it=1, nt; idt=1./(xt(it+1)-xt(it))
307    IF(xt(it+1)<=xs(is+1)) THEN          !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1)
308      CALL mean_lin(xt(it),xt(it+1))
309    ELSE
310      CALL mean_lin(xt(it),xs(is+1)); is=is+1
311      DO WHILE(xs(is+1)<xt(it+1))        !--- 1<=is<=ns-1
312      CALL mean_lin(xs(is),xs(is+1)); is=is+1
313      END DO
314      CALL mean_lin(xs(is),xt(it+1))     !--- 1<=is<=ns
315    END IF
316    IF(xs(is+1)==xt(it+1)) is=is+1     !--- 1<=is<=ns .OR. it==nt
317  END DO
318
319CONTAINS
320
321!-------------------------------------------------------------------------------
322SUBROUTINE mean_lin(a,b)  ! mean [a,b] of the linear function in [xs(is),xs(is+1)]
323!-------------------------------------------------------------------------------
324! Arguments:
325  REAL, INTENT(IN) :: a, b
326!-------------------------------------------------------------------------------
327  IF(lslope.AND.(a/=xs(is).OR.b/=xs(is+1))) THEN; co=0.
328    IF(a==xt(it  )) co=co+xt(it  )-xs(is  )
329    IF(b==xt(it+1)) co=co+xt(it+1)-xs(is+1)
330    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)*(vs(is,:,:,:,:)+co*slope(is,:,:,:,:))
331    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)*(vs(:,is,:,:,:)+co*slope(:,is,:,:,:))
332    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)*(vs(:,:,is,:,:)+co*slope(:,:,is,:,:))
333    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)*(vs(:,:,:,is,:)+co*slope(:,:,:,is,:))
334    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)*(vs(:,:,:,:,is)+co*slope(:,:,:,:,is))
335  ELSE
336    IF(ix==1) vt(it,:,:,:,:) = vt(it,:,:,:,:)+idt*(b-a)* vs(is,:,:,:,:)
337    IF(ix==2) vt(:,it,:,:,:) = vt(:,it,:,:,:)+idt*(b-a)* vs(:,is,:,:,:)
338    IF(ix==3) vt(:,:,it,:,:) = vt(:,:,it,:,:)+idt*(b-a)* vs(:,:,is,:,:)
339    IF(ix==4) vt(:,:,:,it,:) = vt(:,:,:,it,:)+idt*(b-a)* vs(:,:,:,is,:)
340    IF(ix==5) vt(:,:,:,:,it) = vt(:,:,:,:,it)+idt*(b-a)* vs(:,:,:,:,is)
341  END IF
342
343END SUBROUTINE mean_lin
344!-------------------------------------------------------------------------------
345
346END SUBROUTINE regr5_conserv
347!
348!-------------------------------------------------------------------------------
349
350
351!-------------------------------------------------------------------------------
352!
353SUBROUTINE check_size(ix,svs,svt,xs,xt,ns,nt)
354!
355!-------------------------------------------------------------------------------
356! Arguments:
357  INTEGER, INTENT(IN)  :: ix, svs(:), svt(:)
358  REAL,    INTENT(IN)  :: xs(:), xt(:)
359  INTEGER, INTENT(OUT) :: ns,    nt
360!-------------------------------------------------------------------------------
361! Local variables:
362  INTEGER :: rk, is, ii, n
363  CHARACTER(LEN=80) :: sub, msg
364!-------------------------------------------------------------------------------
365  WRITE(sub,'(a,2i0,a)')"regr",SIZE(svs),ix,"_conserv"
366  rk=assert_eq(SIZE(svs),SIZE(svt),TRIM(sub)//': inconsistent ranks')
367  WRITE(msg,'(a,2(i0,a))')TRIM(sub)//": ix (",ix,") exceeds field rank (",rk,")"
368  CALL assert(ix>=1.AND.ix<=rk,msg)
369  ii=0
370  DO is=1,rk; IF(is==ix) CYCLE
371    WRITE(msg,'(a,i0)')TRIM(sub)//" n",is
372    n=assert_eq(svs(is),svt(is),msg)
373    ii=ii+1
374  END DO
375  ns=assert_eq(svs(ix),SIZE(xs)-1,TRIM(sub)//" ns")
376  nt=assert_eq(svt(ix),SIZE(xt)-1,TRIM(sub)//" nt")
377
378  !--- Quick check on sort order:
379  CALL assert(xs(1)<xs(2), TRIM(sub)//": xs bad order")
380  CALL assert(xt(1)<xt(2), TRIM(sub)//": xt bad order")
381  CALL assert(xs(1)<=xt(1).AND.xt(nt+1)<=xs(ns+1), TRIM(sub)//": extrapolation")
382
383END SUBROUTINE check_size
384!
385!-------------------------------------------------------------------------------
386
387
388END MODULE regr_conserv_m
Note: See TracBrowser for help on using the repository browser.