source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/conf_dat_m.F90 @ 5122

Last change on this file since 5122 was 5113, checked in by abarral, 4 months ago

Rename modules in misc from *_mod > lmdz_*
Put cbrt.f90, ch*.f90, pch*.f90 in new lmdz_libmath_pch.f90
File size: 8.9 KB
Line 
1MODULE conf_dat_m
2
3!*******************************************************************************
4
5  PRIVATE
6  PUBLIC :: conf_dat2d, conf_dat3d
7
8
9CONTAINS
10
11
12!-------------------------------------------------------------------------------
13
14SUBROUTINE conf_dat2d(title, xd, yd, xf, yf, champd, interbar)
15
16!-------------------------------------------------------------------------------
17! Author: P. Le Van
18!-------------------------------------------------------------------------------
19! Purpose: Configure the 2D data field "champd" so that:
20!     - Longitudes are in [ -pi     pi   ]
21!     - Latitudes  are in [  pi/2. -pi/2 ]
22!  * xd / yd are initial lon / lats.
23!  * xf / yf are output  lon / lats, possibly modified to satisfy configuration.
24!  * interbar is TRUE for barycentric interpolation.
25!-------------------------------------------------------------------------------
26  USE lmdz_assert_eq, ONLY: assert_eq
27  IMPLICIT NONE
28!-------------------------------------------------------------------------------
29! Arguments:
30  CHARACTER(LEN=*), INTENT(IN)    :: title
31  REAL,             INTENT(IN)    :: xd(:), yd(:)  ! dim (lons) (lats)
32  REAL,             INTENT(INOUT) :: xf(:), yf(:)  ! dim (lons) (lats)
33  REAL,             INTENT(INOUT) :: champd(:,:)   ! dim (lons,lats)
34  LOGICAL,          INTENT(IN)    :: interbar
35!-------------------------------------------------------------------------------
36! Local variables:
37  CHARACTER(LEN=256) :: modname="conf_dat2d"
38  INTEGER :: i, j, ip180, ind, lons, lats
39  LOGICAL :: radlon, invlon ,radlat, invlat
40  REAL    :: pi, pis2, depi, rlatmin, rlatmax, oldxd1
41  REAL, ALLOCATABLE :: xtemp(:) , ytemp(:), champf(:,:)
42!-------------------------------------------------------------------------------
43  lons=assert_eq(SIZE(xd),SIZE(xf),SIZE(champd,1),TRIM(modname)//" lons")
44  lats=assert_eq(SIZE(yd),SIZE(yf),SIZE(champd,2),TRIM(modname)//" lats")
45  ALLOCATE(xtemp(lons),ytemp(lats)); xtemp(:)=xd(:); ytemp(:)=yd(:)
46  ALLOCATE(champf(lons,lats))
47  pi = 2. * ASIN(1.)
48  pis2 = pi/2.
49  depi = 2. * pi
50  radlon=.FALSE.; invlon=.FALSE.
51  IF     (xtemp(1)>=-pi-0.5.AND.xtemp(lons)<=  pi+0.5) THEN
52    radlon = .TRUE.;  invlon = .FALSE.
53  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=depi+0.5) THEN
54    radlon = .TRUE.;  invlon = .TRUE.
55  ELSE IF(xtemp(1)>= -180.5.AND.xtemp(lons)<=   180.5) THEN
56    radlon = .FALSE.; invlon = .FALSE.
57  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=   360.5) THEN
58    radlon = .FALSE.; invlon = .TRUE.
59  ELSE; WRITE(6,*) 'Problems with data longitudes for file '//TRIM(title)
60  END IF
61  invlat = ytemp(1)<ytemp(lats)
62  rlatmin = MIN( ytemp(1), ytemp(lats) )
63  rlatmax = MAX( ytemp(1), ytemp(lats) )
64     
65  IF     (rlatmin>=-pis2-0.5.AND.rlatmax<=pis2+0.5) THEN; radlat = .TRUE.
66  ELSE IF(rlatmin>= -90.-0.5.AND.rlatmax<= 90.+0.5) THEN; radlat = .FALSE.
67  ELSE; WRITE(6,*) 'Problems with data latitudes for file '//TRIM(title)
68  END IF
69
70  IF(.NOT.radlon) xtemp(:)=xtemp(:)*pi/180.
71  IF(.NOT.radlat) ytemp(:)=ytemp(:)*pi/180.
72
73!--- FLIPPED LONGITUDES
74  IF(invlon) THEN
75    champf(:,:)=champd(:,:); xf(:)=xtemp(:)
76
77  !--- Longitudes rotated to get them in [ -pi pi] interval.
78    DO i=1,lons;      IF(xf(i)>pi) EXIT;                   END DO; ip180 = i
79    DO i=1,lons;      IF(xf(i)>pi) xf(i)=xf(i)-depi;       END DO
80    DO i= ip180,lons; ind=i-ip180+1; xtemp(ind)=xf(i);     END DO
81    DO i= ind+1,lons;                xtemp(i  )=xf(i-ind); END DO
82
83  !--- Longitudes rotated in champf
84    DO j=1,lats
85      DO i=ip180,lons; ind=i-ip180+1; champd(ind,j)=champf(i    ,j); END DO
86      DO i=ind+1,lons;                champd(i  ,j)=champf(i-ind,j); END DO
87    END DO
88  END IF
89
90!--- FLIPPED LATITUDES
91  IF(invlat) THEN
92    yf(:)=ytemp(:)
93    champf(:,:)=champd(:,:)
94    ytemp(lats:1:-1)=yf(:)
95    DO j=1,lats; champd(:,lats-j+1)=champf(:,j); END DO
96  END IF
97
98!--- FOR BARYCENTRIC INTERPOLATION
99  IF(interbar) THEN
100    oldxd1 = xtemp(1)
101    xtemp(1:lons-1)=0.5*(xtemp(1:lons-1)+xtemp(2:lons))
102    xtemp(  lons  )=0.5*(xtemp(  lons  )+oldxd1+depi)
103    ytemp(1:lats-1)=0.5*(ytemp(1:lats-1)+ytemp(2:lats))
104  END IF
105  DEALLOCATE(champf)
106  xf(:)=xtemp(:); DEALLOCATE(xtemp)
107  yf(:)=ytemp(:); DEALLOCATE(ytemp)
108
109END SUBROUTINE conf_dat2d
110
111!-------------------------------------------------------------------------------
112
113
114!-------------------------------------------------------------------------------
115
116SUBROUTINE conf_dat3d(title, xd, yd, zd, xf, yf, zf, champd, interbar)
117
118!-------------------------------------------------------------------------------
119! Author: P. Le Van
120!-------------------------------------------------------------------------------
121! Purpose: Configure the 3D data field "champd" so that:
122!     - Longitudes are in [ -pi     pi   ]
123!     - Latitudes  are in [  pi/2. -pi/2 ]
124!     - Vertical levels from ground to model top (in Pascals)
125!  * xd / yd are initial lon / lats.
126!  * xf / yf are output  lon / lats, possibly modified to satisfy configuration.
127!  * zf      are output pressures
128!  * interbar is TRUE for barycentric interpolation.
129!-------------------------------------------------------------------------------
130  USE lmdz_assert_eq, ONLY: assert_eq
131  IMPLICIT NONE
132!-------------------------------------------------------------------------------
133! Arguments:
134  CHARACTER(LEN=*), INTENT(IN)    :: title
135  REAL,             INTENT(IN)    :: xd(:), yd(:), zd(:) ! (lons) (lats) (levs)
136  REAL,             INTENT(INOUT) :: xf(:), yf(:), zf(:) ! (lons) (lats) (levs)
137  REAL,             INTENT(INOUT) :: champd(:,:,:)       ! (lons,lats,levs)
138  LOGICAL,          INTENT(IN)    :: interbar
139!-------------------------------------------------------------------------------
140! Local variables:
141  CHARACTER(LEN=256) :: modname="conf_dat3d"
142  INTEGER :: i, j, l, ip180, ind, lons, lats, levs
143  LOGICAL :: radlon, invlon ,radlat, invlat, invlev
144  REAL    :: pi, pis2, depi, presmax, rlatmin, rlatmax, oldxd1
145  REAL, ALLOCATABLE :: xtemp(:) , ytemp(:), ztemp(:), champf(:,:,:)
146!-------------------------------------------------------------------------------
147  lons=assert_eq(SIZE(xd),SIZE(xf),SIZE(champd,1),TRIM(modname)//" lons")
148  lats=assert_eq(SIZE(yd),SIZE(yf),SIZE(champd,2),TRIM(modname)//" lats")
149  levs=assert_eq(SIZE(zd),SIZE(zf),SIZE(champd,3),TRIM(modname)//" levs")
150  ALLOCATE(xtemp(lons),ytemp(lats),ztemp(levs),champf(lons,lats,levs))
151  xtemp(:)=xd(:); ytemp(:)=yd(:); ztemp(:)=zd(:)
152  pi = 2. * ASIN(1.)
153  pis2 = pi/2.
154  depi = 2. * pi
155  radlon=.FALSE.; invlon=.FALSE.
156  IF     (xtemp(1)>=-pi-0.5.AND.xtemp(lons)<=  pi+0.5) THEN
157    radlon = .TRUE.;  invlon = .FALSE.
158  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=depi+0.5) THEN
159    radlon = .TRUE.;  invlon = .TRUE.
160  ELSE IF(xtemp(1)>= -180.5.AND.xtemp(lons)<=   180.5) THEN
161    radlon = .FALSE.; invlon = .FALSE.
162  ELSE IF(xtemp(1)>=   -0.5.AND.xtemp(lons)<=   360.5) THEN
163    radlon = .FALSE.; invlon = .TRUE.
164  ELSE; WRITE(6,*) 'Problems with data longitudes for file '//TRIM(title)
165  END IF
166  invlat = ytemp(1)<ytemp(lats)
167  rlatmin = MIN( ytemp(1), ytemp(lats) )
168  rlatmax = MAX( ytemp(1), ytemp(lats) )
169     
170  IF     (rlatmin>=-pis2-0.5.AND.rlatmax<=pis2+0.5) THEN; radlat = .TRUE.
171  ELSE IF(rlatmin>= -90.-0.5.AND.rlatmax<= 90.+0.5) THEN; radlat = .FALSE.
172  ELSE; WRITE(6,*) 'Problems with data latitudes for file '//TRIM(title)
173  END IF
174
175  IF(.NOT.radlon) xtemp(:)=xtemp(:)*pi/180.
176  IF(.NOT.radlat) ytemp(:)=ytemp(:)*pi/180.
177
178!--- FLIPPED LONGITUDES
179  IF(invlon) THEN
180    champf(:,:,:)=champd(:,:,:); xf(:)=xtemp(:)
181
182  !--- Longitudes rotated to get them in [ -pi pi] interval.
183    DO i=1,lons;      IF(xf(i)>pi) EXIT;                   END DO; ip180 = i
184    DO i=1,lons;      IF(xf(i)>pi) xf(i)=xf(i)-depi;       END DO
185    DO i= ip180,lons; ind=i-ip180+1; xtemp(ind)=xf(i);     END DO
186    DO i= ind+1,lons;                xtemp(i  )=xf(i-ind); END DO
187
188  !--- Longitudes rotated in champf
189    DO l=1,levs
190    DO j=1,lats
191      DO i=ip180,lons; ind=i-ip180+1; champd(ind,j,l)=champf(i    ,j,l); END DO
192      DO i=ind+1,lons;                champd(i  ,j,l)=champf(i-ind,j,l); END DO
193    END DO
194    END DO
195  END IF
196
197!--- FLIPPED LATITUDES
198  IF(invlat) THEN
199    yf(:)=ytemp(:)
200    champf(:,:,:)=champd(:,:,:)
201    ytemp(lats:1:-1)=yf(:)
202    DO l=1,levs
203      DO j=1,lats; champd(:,lats-j+1,l)=champf(:,j,l); END DO
204    END DO
205  END IF
206
207!--- FOR BARYCENTRIC INTERPOLATION
208  IF(interbar) THEN
209    oldxd1 = xtemp(1)
210    xtemp(1:lons-1)=0.5*(xtemp(1:lons-1)+xtemp(2:lons))
211    xtemp(  lons  )=0.5*(xtemp(  lons  )+oldxd1+depi)
212    ytemp(1:lats-1)=0.5*(ytemp(1:lats-1)+ytemp(2:lats))
213  END IF
214
215!--- FLIPPED LEVELS
216  invlev=ztemp(1)<ztemp(levs)
217  IF(MAX(ztemp(1),ztemp(levs))<1200.) ztemp(:)=ztemp(:)*100.
218  IF(invlev) THEN
219    zf(:)=ztemp(:)
220    champf(:,:,:)=champd(:,:,:)
221    ztemp(levs:1:-1)=zf(:)
222    DO l=1,levs; champd(:,:,levs+1-l)=champf(:,:,l); END DO
223  END IF
224
225  DEALLOCATE(champf)
226  xf(:)=xtemp(:); DEALLOCATE(xtemp)
227  yf(:)=ytemp(:); DEALLOCATE(ytemp)
228  zf(:)=ztemp(:); DEALLOCATE(ztemp)
229
230END SUBROUTINE conf_dat3d
231
232!-------------------------------------------------------------------------------
233
234END MODULE conf_dat_m
235
Note: See TracBrowser for help on using the repository browser.