source: LMDZ5/trunk/libf/misc/slopes_m.F90 @ 5086

Last change on this file since 5086 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: 9.9 KB
RevLine 
[2788]1MODULE slopes_m
[2440]2
3  ! Author: Lionel GUEZ
[2788]4  ! Extension / factorisation: David CUGNET
[2440]5
[2788]6  IMPLICIT NONE
[2440]7
[2788]8  ! Those generic function computes second order slopes with Van
9  ! Leer slope-limiting, given cell averages. Reference: Dukowicz,
10  ! 1987, SIAM Journal on Scientific and Statistical Computing, 8,
11  ! 305.
[2440]12
[2788]13  ! The only difference between the specific functions is the rank
14  ! of the first argument and the equal rank of the result.
[2440]15
[2788]16  ! slope(ix,...) acts on ix th dimension.
[2440]17
[2788]18  ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1
19  ! real, intent(in):: x(:) ! (n + 1) cell edges
20  ! real slopes, same shape as f ! (n, ...)
21  INTERFACE slopes
22     MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5
23  END INTERFACE
[2440]24
[2788]25  PRIVATE
26  PUBLIC :: slopes
[2440]27
[2788]28CONTAINS
[2440]29
[2788]30!-------------------------------------------------------------------------------
31!
32PURE FUNCTION slopes1(ix, f, x)
33!
34!-------------------------------------------------------------------------------
35! Arguments:
36  INTEGER, INTENT(IN) :: ix
37  REAL,    INTENT(IN) :: f(:)
38  REAL,    INTENT(IN) :: x(:)
39  REAL :: slopes1(SIZE(f,1))
40!-------------------------------------------------------------------------------
41! Local:
42  INTEGER :: n, i, j, sta(2), sto(2)
43  REAL :: xc(SIZE(f,1))                             ! (n) cell centers
44  REAL :: h(2:SIZE(f,1)-1), delta_xc(2:SIZE(f,1)-1) ! (2:n-1)
45  REAL :: fm, ff, fp, dx
46!-------------------------------------------------------------------------------
47  n=SIZE(f,ix)
48  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
49  FORALL(i=2:n-1)
50    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
51  END FORALL
52  slopes1(:)=0.
53  DO i=2,n-1
54    ff=f(i); fm=f(i-1); fp=f(i+1)
55    IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
56      slopes1(i)=0.; CYCLE           !--- Local extremum
57      !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
58      slopes1(i)=(fp-fm)/delta_xc(i)
59      !--- Slope limitation
60      slopes1(i) = SIGN(MIN(ABS(slopes1(i)), &
61        ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes1(i) )
62     END IF
63  END DO
[2440]64
[2788]65END FUNCTION slopes1
66!
67!-------------------------------------------------------------------------------
[2440]68
69
[2788]70!-------------------------------------------------------------------------------
71!
72PURE FUNCTION slopes2(ix, f, x)
73!
74!-------------------------------------------------------------------------------
75! Arguments:
76  INTEGER, INTENT(IN) :: ix
77  REAL,    INTENT(IN) :: f(:, :)
78  REAL,    INTENT(IN) :: x(:)
79  REAL :: slopes2(SIZE(f,1),SIZE(f,2))
80!-------------------------------------------------------------------------------
81! Local:
82  INTEGER :: n, i, j, sta(2), sto(2)
83  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
84  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
85  REAL :: fm, ff, fp, dx
86!-------------------------------------------------------------------------------
87  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
88  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
89  FORALL(i=2:n-1)
90    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
91  END FORALL
92  slopes2(:,:)=0.
93  sta=[1,1]; sta(ix)=2
94  sto=SHAPE(f); sto(ix)=n-1
95  DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
96    DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
97      ff=f(i,j)
98      SELECT CASE(ix)
99        CASE(1); fm=f(i-1,j); fp=f(i+1,j)
100        CASE(2); fm=f(i,j-1); fp=f(i,j+1)
101      END SELECT
102      IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
103        slopes2(i,j)=0.; CYCLE           !--- Local extremum
104        !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
105        slopes2(i,j)=(fp-fm)/dx
106        !--- Slope limitation
107        slopes2(i,j) = SIGN(MIN(ABS(slopes2(i,j)), &
108          ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes2(i,j) )
109       END IF
110    END DO
111  END DO
112  DEALLOCATE(xc,h,delta_xc)
[2440]113
[2788]114END FUNCTION slopes2
115!
116!-------------------------------------------------------------------------------
[2440]117
118
[2788]119!-------------------------------------------------------------------------------
120!
121PURE FUNCTION slopes3(ix, f, x)
122!
123!-------------------------------------------------------------------------------
124! Arguments:
125  INTEGER, INTENT(IN) :: ix
126  REAL,    INTENT(IN) :: f(:, :, :)
127  REAL,    INTENT(IN) :: x(:)
128  REAL :: slopes3(SIZE(f,1),SIZE(f,2),SIZE(f,3))
129!-------------------------------------------------------------------------------
130! Local:
131  INTEGER :: n, i, j, k, sta(3), sto(3)
132  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
133  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
134  REAL :: fm, ff, fp, dx
135!-------------------------------------------------------------------------------
136  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
137  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
138  FORALL(i=2:n-1)
139    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
140  END FORALL
141  slopes3(:,:,:)=0.
142  sta=[1,1,1]; sta(ix)=2
143  sto=SHAPE(f); sto(ix)=n-1
144  DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
145    DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
146      DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
147        ff=f(i,j,k)
148        SELECT CASE(ix)
149          CASE(1); fm=f(i-1,j,k); fp=f(i+1,j,k)
150          CASE(2); fm=f(i,j-1,k); fp=f(i,j+1,k)
151          CASE(3); fm=f(i,j,k-1); fp=f(i,j,k+1)
152        END SELECT
153        IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
154          slopes3(i,j,k)=0.; CYCLE           !--- Local extremum
155          !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
156          slopes3(i,j,k)=(fp-fm)/dx
157          !--- Slope limitation
158          slopes3(i,j,k) = SIGN(MIN(ABS(slopes3(i,j,k)), &
159            ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes3(i,j,k) )
160         END IF
161      END DO
162    END DO
163  END DO
164  DEALLOCATE(xc,h,delta_xc)
[2440]165
[2788]166END FUNCTION slopes3
167!
168!-------------------------------------------------------------------------------
[2440]169
170
[2788]171!-------------------------------------------------------------------------------
172!
173PURE FUNCTION slopes4(ix, f, x)
174!
175!-------------------------------------------------------------------------------
176! Arguments:
177  INTEGER, INTENT(IN) :: ix
178  REAL,    INTENT(IN) :: f(:, :, :, :)
179  REAL,    INTENT(IN) :: x(:)
180  REAL :: slopes4(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4))
181!-------------------------------------------------------------------------------
182! Local:
183  INTEGER :: n, i, j, k, l, m, sta(4), sto(4)
184  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
185  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
186  REAL :: fm, ff, fp, dx
187!-------------------------------------------------------------------------------
188  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
189  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
190  FORALL(i=2:n-1)
191    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
192  END FORALL
193  slopes4(:,:,:,:)=0.
194  sta=[1,1,1,1]; sta(ix)=2
195  sto=SHAPE(f); sto(ix)=n-1
196  DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l)
197    DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
198      DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
199        DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
200          ff=f(i,j,k,l)
201          SELECT CASE(ix)
202            CASE(1); fm=f(i-1,j,k,l); fp=f(i+1,j,k,l)
203            CASE(2); fm=f(i,j-1,k,l); fp=f(i,j+1,k,l)
204            CASE(3); fm=f(i,j,k-1,l); fp=f(i,j,k+1,l)
205            CASE(4); fm=f(i,j,k,l-1); fp=f(i,j,k,l+1)
206          END SELECT
207          IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
208            slopes4(i,j,k,l)=0.; CYCLE           !--- Local extremum
209            !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
210            slopes4(i,j,k,l)=(fp-fm)/dx
211            !--- Slope limitation
212            slopes4(i,j,k,l) = SIGN(MIN(ABS(slopes4(i,j,k,l)), &
213              ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes4(i,j,k,l) )
214           END IF
215        END DO
216      END DO
217    END DO
218  END DO
219  DEALLOCATE(xc,h,delta_xc)
[2440]220
[2788]221END FUNCTION slopes4
222!
223!-------------------------------------------------------------------------------
[2440]224
225
[2788]226!-------------------------------------------------------------------------------
227!
228PURE FUNCTION slopes5(ix, f, x)
229!
230!-------------------------------------------------------------------------------
231! Arguments:
232  INTEGER, INTENT(IN) :: ix
233  REAL,    INTENT(IN) :: f(:, :, :, :, :)
234  REAL,    INTENT(IN) :: x(:)
235  REAL :: slopes5(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4),SIZE(f,5))
236!-------------------------------------------------------------------------------
237! Local:
238  INTEGER :: n, i, j, k, l, m, sta(5), sto(5)
239  REAL, ALLOCATABLE :: xc(:)                        ! (n) cell centers
240  REAL, ALLOCATABLE :: h(:), delta_xc(:)            ! (2:n-1)
241  REAL :: fm, ff, fp, dx
242!-------------------------------------------------------------------------------
243  n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-1))
244  FORALL(i=1:n) xc(i)=(x(i)+x(i+1))/2.
245  FORALL(i=2:n-1)
246    h(i)=ABS(x(i+1)-xc(i)) ; delta_xc(i)=xc(i+1)-xc(i-1)
247  END FORALL
248  slopes5(:,:,:,:,:)=0.
249  sta=[1,1,1,1,1]; sta(ix)=2
250  sto=SHAPE(f);    sto(ix)=n-1
251  DO m=sta(5),sto(5); IF(ix==5) dx=delta_xc(m)
252    DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l)
253      DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k)
254        DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j)
255          DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i)
256            ff=f(i,j,k,l,m)
257            SELECT CASE(ix)
258              CASE(1); fm=f(i-1,j,k,l,m); fp=f(i+1,j,k,l,m)
259              CASE(2); fm=f(i,j-1,k,l,m); fp=f(i,j+1,k,l,m)
260              CASE(3); fm=f(i,j,k-1,l,m); fp=f(i,j,k+1,l,m)
261              CASE(4); fm=f(i,j,k,l-1,m); fp=f(i,j,k,l+1,m)
262              CASE(5); fm=f(i,j,k,l,m-1); fp=f(i,j,k,l,m+1)
263            END SELECT
264            IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN
265              slopes5(i,j,k,l,m)=0.; CYCLE           !--- Local extremum
266              !--- 2nd order slope ; (fm, ff, fp) strictly monotonous
267              slopes5(i,j,k,l,m)=(fp-fm)/dx
268              !--- Slope limitation
269              slopes5(i,j,k,l,m) = SIGN(MIN(ABS(slopes5(i,j,k,l,m)), &
270                ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes5(i,j,k,l,m) )
271            END IF
272          END DO
273        END DO
274      END DO
275    END DO
276  END DO
277  DEALLOCATE(xc,h,delta_xc)
[2440]278
[2788]279END FUNCTION slopes5
280!
281!-------------------------------------------------------------------------------
[2440]282
[2788]283END MODULE slopes_m
Note: See TracBrowser for help on using the repository browser.