source: LMDZ5/trunk/libf/misc/regr_lint_m.F90 @ 5305

Last change on this file since 5305 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: 8.1 KB
RevLine 
[2788]1MODULE regr_lint_m
2
3  USE assert_eq_m,   ONLY: assert_eq
4  USE assert_m,      ONLY: assert
5  USE interpolation, ONLY: hunt
6
7  IMPLICIT NONE
8
9! Purpose: Each procedure regrids by linear interpolation along dimension "ix"
10!          the input field "vs" to the output field "vt".
11! Remark:
12!   * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt)
13
14! Argument                         Type         Description
15!-------------------------------------------------------------------------------
16!  INTEGER, INTENT(IN)  :: ix      Scalar       dimension regridded <=rank(vs)
17!  REAL,    INTENT(IN)  :: vs(*)   Rank>=1      source grid field values
18!  REAL,    INTENT(IN)  :: xs(:)   Vector(ns)   centers of source grid, asc. order
19!  REAL,    INTENT(IN)  :: xt(:)   Vector(nt)   centers of target grid, asc. order
20!  REAL,    INTENT(OUT) :: vt(*)   Rank>=1      regridded field
21
22  INTERFACE regr_lint
23    ! The procedures differ only from the rank of the input/output fields.
24    MODULE PROCEDURE regr1_lint, regr2_lint, regr3_lint, &
25                     regr4_lint, regr5_lint
26  END INTERFACE
27
28  PRIVATE
29  PUBLIC :: regr_lint
30
31CONTAINS
32
33
34!-------------------------------------------------------------------------------
35!
36SUBROUTINE regr1_lint(ix, vs, xs, xt, vt)
37!
38!-------------------------------------------------------------------------------
39! Arguments:
40  INTEGER, INTENT(IN)  :: ix
41  REAL,    INTENT(IN)  :: vs(:)
42  REAL,    INTENT(IN)  :: xs(:)
43  REAL,    INTENT(IN)  :: xt(:)
44  REAL,    INTENT(OUT) :: vt(:)
45!-------------------------------------------------------------------------------
46! Local variables:
47  INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
48  REAL    :: r
49!-------------------------------------------------------------------------------
50  CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt)
51  is = -1 ! go immediately to bisection on first call to "hunt"
52  DO it=1,SIZE(xt)
53    CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1)
54    r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb))
55    vt(it)=r*vs(isb)+(1.-r)*vs(isb+1)
56  END DO
57
58END SUBROUTINE regr1_lint
59!
60!-------------------------------------------------------------------------------
61
62
63!-------------------------------------------------------------------------------
64!
65SUBROUTINE regr2_lint(ix, vs, xs, xt, vt)
66!
67!-------------------------------------------------------------------------------
68! Arguments:
69  INTEGER, INTENT(IN)  :: ix
70  REAL,    INTENT(IN)  :: vs(:,:)
71  REAL,    INTENT(IN)  :: xs(:)
72  REAL,    INTENT(IN)  :: xt(:)
73  REAL,    INTENT(OUT) :: vt(:,:)
74!-------------------------------------------------------------------------------
75! Local variables:
76  INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
77  REAL    :: r
78!-------------------------------------------------------------------------------
79  CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt)
80  is = -1 ! go immediately to bisection on first call to "hunt"
81  DO it=1,SIZE(xt)
82    CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1)
83    r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb))
84    IF(ix==1) vt(it,:)=r*vs(isb,:)+(1.-r)*vs(isb+1,:)
85    IF(ix==2) vt(:,it)=r*vs(:,isb)+(1.-r)*vs(:,isb+1)
86  END DO
87
88END SUBROUTINE regr2_lint
89!
90!-------------------------------------------------------------------------------
91
92
93!-------------------------------------------------------------------------------
94!
95SUBROUTINE regr3_lint(ix, vs, xs, xt, vt)
96!
97!-------------------------------------------------------------------------------
98! Arguments:
99  INTEGER, INTENT(IN)  :: ix
100  REAL,    INTENT(IN)  :: vs(:,:,:)
101  REAL,    INTENT(IN)  :: xs(:)
102  REAL,    INTENT(IN)  :: xt(:)
103  REAL,    INTENT(OUT) :: vt(:,:,:)
104!-------------------------------------------------------------------------------
105! Local variables:
106  INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
107  REAL    :: r
108!-------------------------------------------------------------------------------
109  CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt)
110  is = -1 ! go immediately to bisection on first call to "hunt"
111  DO it=1,SIZE(xt)
112    CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1)
113    r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb))
114    IF(ix==1) vt(it,:,:)=r*vs(isb,:,:)+(1.-r)*vs(isb+1,:,:)
115    IF(ix==2) vt(:,it,:)=r*vs(:,isb,:)+(1.-r)*vs(:,isb+1,:)
116    IF(ix==3) vt(:,:,it)=r*vs(:,:,isb)+(1.-r)*vs(:,:,isb+1)
117  END DO
118
119END SUBROUTINE regr3_lint
120!
121!-------------------------------------------------------------------------------
122
123
124!-------------------------------------------------------------------------------
125!
126SUBROUTINE regr4_lint(ix, vs, xs, xt, vt)
127!
128!-------------------------------------------------------------------------------
129! Arguments:
130  INTEGER, INTENT(IN)  :: ix
131  REAL,    INTENT(IN)  :: vs(:,:,:,:)
132  REAL,    INTENT(IN)  :: xs(:)
133  REAL,    INTENT(IN)  :: xt(:)
134  REAL,    INTENT(OUT) :: vt(:,:,:,:)
135!-------------------------------------------------------------------------------
136! Local variables:
137  INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
138  REAL    :: r
139!-------------------------------------------------------------------------------
140  CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt)
141  is = -1 ! go immediately to bisection on first call to "hunt"
142  DO it=1,SIZE(xt)
143    CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1)
144    r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb))
145    IF(ix==1) vt(it,:,:,:)=r*vs(isb,:,:,:)+(1.-r)*vs(isb+1,:,:,:)
146    IF(ix==2) vt(:,it,:,:)=r*vs(:,isb,:,:)+(1.-r)*vs(:,isb+1,:,:)
147    IF(ix==3) vt(:,:,it,:)=r*vs(:,:,isb,:)+(1.-r)*vs(:,:,isb+1,:)
148    IF(ix==4) vt(:,:,:,it)=r*vs(:,:,:,isb)+(1.-r)*vs(:,:,:,isb+1)
149  END DO
150
151END SUBROUTINE regr4_lint
152!
153!-------------------------------------------------------------------------------
154
155
156!-------------------------------------------------------------------------------
157!
158SUBROUTINE regr5_lint(ix, vs, xs, xt, vt)
159!
160!-------------------------------------------------------------------------------
161! Arguments:
162  INTEGER, INTENT(IN)  :: ix
163  REAL,    INTENT(IN)  :: vs(:,:,:,:,:)
164  REAL,    INTENT(IN)  :: xs(:)
165  REAL,    INTENT(IN)  :: xt(:)
166  REAL,    INTENT(OUT) :: vt(:,:,:,:,:)
167!-------------------------------------------------------------------------------
168! Local variables:
169  INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1"
170  REAL    :: r
171!-------------------------------------------------------------------------------
172  CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt)
173  is = -1 ! go immediately to bisection on first call to "hunt"
174  DO it=1,SIZE(xt)
175    CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1)
176    r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb))
177    IF(ix==1) vt(it,:,:,:,:)=r*vs(isb,:,:,:,:)+(1.-r)*vs(isb+1,:,:,:,:)
178    IF(ix==2) vt(:,it,:,:,:)=r*vs(:,isb,:,:,:)+(1.-r)*vs(:,isb+1,:,:,:)
179    IF(ix==3) vt(:,:,it,:,:)=r*vs(:,:,isb,:,:)+(1.-r)*vs(:,:,isb+1,:,:)
180    IF(ix==4) vt(:,:,:,it,:)=r*vs(:,:,:,isb,:)+(1.-r)*vs(:,:,:,isb+1,:)
181    IF(ix==5) vt(:,:,:,:,it)=r*vs(:,:,:,:,isb)+(1.-r)*vs(:,:,:,:,isb+1)
182  END DO
183
184END SUBROUTINE regr5_lint
185!
186!-------------------------------------------------------------------------------
187
188
189!-------------------------------------------------------------------------------
190!
191SUBROUTINE check_size(ix,svs,svt,nxs,nxt,ns,nt)
192!
193!-------------------------------------------------------------------------------
194! Arguments:
195  INTEGER, INTENT(IN)  :: ix, svs(:), svt(:), nxs, nxt
196  INTEGER, INTENT(OUT) :: ns, nt
197!-------------------------------------------------------------------------------
198! Local variables:
199  INTEGER :: rk, is
200  CHARACTER(LEN=80) :: sub, msg
201!-------------------------------------------------------------------------------
202  rk=SIZE(svs)
203  WRITE(sub,'(a,2i0,a)')"regr",rk,ix,"_lint"
204  CALL assert(ix>=1.AND.ix<=rk,TRIM(sub)//": ix exceeds fields rank")
205  DO is=1,rk; IF(is==ix) CYCLE
206    WRITE(msg,'(a,i1)')TRIM(sub)//" n",is
207    CALL assert(svs(is)==svt(is),msg)
208  END DO
209  ns=assert_eq(svs(ix),nxs,TRIM(sub)//" ns")
210  nt=assert_eq(svt(ix),nxt,TRIM(sub)//" nt")
211
212END SUBROUTINE check_size
213!
214!-------------------------------------------------------------------------------
215
216END MODULE regr_lint_m
217!
218!-------------------------------------------------------------------------------
219
Note: See TracBrowser for help on using the repository browser.