source: LMDZ5/branches/AI-cosp/libf/dynphy_lonlat/grid_atob_m.f90 @ 3485

Last change on this file since 3485 was 2361, checked in by dcugnet, 9 years ago

In etat0dyn: removed few useless lines: "masque" is always known because etat0dyn is called after etat0phys.
In grid_atob: shape error in routine fine2coarse now fixed: "msk" argument and local variable mask must have the dimensions of the output grid. Working unit of dist_sphe is no longer degree, but radian.
File size: 12.1 KB
Line 
1!*******************************************************************************
2!
3MODULE grid_atob_m
4!
5!*******************************************************************************
6
7  USE assert_eq_m, ONLY: assert_eq
8
9  PRIVATE
10  PUBLIC :: grille_m, rugosite, sea_ice, rugsoro
11
12CONTAINS
13
14!-------------------------------------------------------------------------------
15!
16SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2)
17!
18!-------------------------------------------------------------------------------
19  IMPLICIT NONE
20!-------------------------------------------------------------------------------
21! Arguments:
22  REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
23  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
24  REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
25  REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
26  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
27  REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
28!-------------------------------------------------------------------------------
29! Local variables:
30  CHARACTER(LEN=256) :: modname="fine2coarse"
31  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
32  LOGICAL :: li, lo
33  REAL    :: inc
34  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
35  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
36  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
37  REAL,    DIMENSION(SIZE(x_o))           :: a, b
38  REAL,    DIMENSION(SIZE(y_o))           :: c, d
39  REAL,    PARAMETER :: thresh=1.E-5
40!-------------------------------------------------------------------------------
41  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
42  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
43  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
44  li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
45  lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
46  mi=assert_eq(mi,m1,TRIM(modname)//" mi")
47  ni=assert_eq(ni,n1,TRIM(modname)//" ni")
48  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo")
49  no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no")
50  mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)
51
52!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
53  b(mo)=x_o(mo)+(x_o(mo)-x_o(mo-1))/2.0; b(1:mo-1)=(x_o(1:mo-1)+x_o(2:mo))/2.0
54  d(no)=y_o(no)+(y_o(no)-y_o(no-1))/2.0; d(1:no-1)=(y_o(1:no-1)+y_o(2:no))/2.0
55  a(1 )=x_o(1 )-(x_o(2 )-x_o(1   ))/2.0; a(2:mo  )=   b(1:mo-1)
56  c(1 )=y_o(1 )-(y_o(2 )-y_o(1   ))/2.0; c(2:no  )=   d(1:no-1)
57
58!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID
59  d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0
60  IF(lo) d_o2(:,:)=0.
61  DO ji = 1, ni
62    DO ii = 1, mi
63      IF(li) inc=d_i(ii,ji)
64      DO jo = 1, no
65        IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND.   &
66           (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE
67        DO io = 1, mo
68          IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. &
69             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE
70          num_tot(io,jo)=num_tot(io,jo)+1
71          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
72          IF(.NOT.lo) CYCLE
73          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
74        END DO
75      END DO
76    END DO
77  END DO
78
79!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL
80  found(:,:)=num_tot(:,:)/=0
81  WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:))
82  IF(PRESENT(d_o2)) THEN
83    WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:))
84    RETURN
85  END IF
86  nn=COUNT(found); IF(nn==0) RETURN
87
88!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2)
89  DO io = 1, mo
90    DO jo = 1, no
91      IF(found(io,jo)) CYCLE
92!      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
93      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
94      nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
95      inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
96      IF(mask(io,jo)) d_o1(io,jo)=inc
97    END DO
98  END DO
99
100END SUBROUTINE fine2coarse
101!
102!-------------------------------------------------------------------------------
103
104
105!-------------------------------------------------------------------------------
106!
107SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie)
108!
109!-------------------------------------------------------------------------------
110! Author:  Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead)
111!-------------------------------------------------------------------------------
112! Purpose: Naive method to regrid on a coarser grid.
113!   Value at a new point is an average of the old points lcoated in a zone
114!   surrounding the considered new point.
115!   No ponderation is used (see grille_p)
116!
117!           (c)
118!        ----d-----
119!        | . . . .|
120!        |        |
121!     (b)a . * . .b(a)
122!        |        |
123!        | . . . .|
124!        ----c-----
125!           (d)
126!
127!-------------------------------------------------------------------------------
128  IMPLICIT NONE
129!-------------------------------------------------------------------------------
130! Arguments:
131  REAL, INTENT(IN)  :: xdata(:), ydata(:)       !--- INPUT  FIELD X AND Y COORD.
132  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
133  REAL, INTENT(IN)  :: x(:), y(:)               !--- OUTPUT FIELD X AND Y COORD.
134  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y))  !--- OUTPUT FIELD
135!-------------------------------------------------------------------------------
136  CALL fine2coarse(xdata,ydata,x,y,sortie,entree)
137
138END SUBROUTINE grille_m
139!
140!-------------------------------------------------------------------------------
141
142
143!-------------------------------------------------------------------------------
144!
145SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask)
146!
147!-------------------------------------------------------------------------------
148! Author:  Z.X. Li (april 1st 1994)
149!-------------------------------------------------------------------------------
150! Purpose: Remap rugosity length ; constant value (0.001) on oceans.
151! Naive method  (see grille_m)
152!-------------------------------------------------------------------------------
153  IMPLICIT NONE
154!-------------------------------------------------------------------------------
155! Arguments:
156  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
157  REAL, INTENT(IN)  :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT  FIELD
158  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
159  REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD
160  REAL, INTENT(IN)  :: mask  (SIZE(x),SIZE(y)) !--- MASK
161!-------------------------------------------------------------------------------
162  CALL fine2coarse(xdata,ydata,x,y,sortie,LOG(entree))
163  WHERE(NINT(mask)==1)
164    sortie(:,:)=EXP(sortie(:,:))
165  ELSE WHERE
166    sortie(:,:)=0.001
167  END WHERE
168
169END SUBROUTINE rugosite
170!
171!-------------------------------------------------------------------------------
172
173
174!-------------------------------------------------------------------------------
175!
176SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice)
177!
178!-------------------------------------------------------------------------------
179! Author:  Z.X. Li (april 1st 1994)
180! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract.
181! field (between 0 and 1).
182! Naive method  (see grille_m)
183!-------------------------------------------------------------------------------
184  IMPLICIT NONE
185!-------------------------------------------------------------------------------
186! Arguments:
187  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
188  REAL, INTENT(IN)  :: glace01(:,:)            !--- INPUT  FIELD
189  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
190  REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD
191!-------------------------------------------------------------------------------
192  CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=NINT(glace01)==1)
193
194END SUBROUTINE sea_ice
195!
196!-------------------------------------------------------------------------------
197
198
199!-------------------------------------------------------------------------------
200!
201SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs)
202!
203!-------------------------------------------------------------------------------
204! Author:  Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995.
205! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell
206!-------------------------------------------------------------------------------
207  IMPLICIT NONE
208!-------------------------------------------------------------------------------
209! Arguments:
210  REAL, INTENT(IN)  :: xrel(:), yrel(:)        !--- INPUT  FIELD X AND Y COORD.
211  REAL, INTENT(IN)  :: relief(:,:)             !--- INPUT  FIELD
212  REAL, INTENT(IN)  :: xmod(:), ymod(:)        !--- OUTPUT FIELD X AND Y COORD.
213  REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD
214!-------------------------------------------------------------------------------
215! Local variable:
216  INTEGER           :: k, nn
217  INTEGER, PARAMETER:: itmp=360, jtmp=180
218  REAL  :: out(SIZE(xmod),SIZE(xmod)), amin, amax
219  REAL  :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp)
220!-------------------------------------------------------------------------------
221
222!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION
223  xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)]
224  ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)]
225  CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp)
226  cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2
227  nn=COUNT(cham2tmp<=-7.5)
228  IF(nn/=0) THEN
229    PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn
230    CALL ABORT_GCM("", "", 1)
231  END IF
232  nn=COUNT(cham2tmp<0)
233  IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn
234  WHERE(cham2tmp<0.0) cham2tmp=0.0
235  cham2tmp(:,:)=SQRT(cham2tmp(:,:))
236  amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp)
237  PRINT*, 'Ecart-type 1x1:', amin, amax
238
239!--- COMPUTE RUGOSITY AT REQUIRED SCALE
240  WHERE(cham2tmp<0.001) cham2tmp=0.001
241  CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp)))
242  out=EXP(out)
243  amin=MINVAL(out); amax=MAXVAL(out)
244  PRINT*, 'Ecart-type du modele:', amin, amax
245  out=out/amax*20.0
246  amin=MINVAL(out); amax=MAXVAL(out)
247  PRINT*, 'Longueur de rugosite du modele:', amin, amax
248  rugs=REAL(out)
249
250END SUBROUTINE rugsoro
251!
252!-------------------------------------------------------------------------------
253
254
255!-------------------------------------------------------------------------------
256!
257SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
258!
259!-------------------------------------------------------------------------------
260! Author:  Laurent Li (december 30th 1996).
261! Purpose: Compute min. distance (along big circle) between 2 points in radians.
262!-------------------------------------------------------------------------------
263  IMPLICIT NONE
264!-------------------------------------------------------------------------------
265! Arguments:
266  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
267  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
268  REAL, INTENT(OUT):: distance(SIZE(rlon),SIZE(rlat)) !--- Distance    (radians)
269!-------------------------------------------------------------------------------
270! Local variables:
271  LOGICAL, SAVE :: first=.TRUE.
272  REAL,    SAVE :: pi, hpi
273  REAL    :: pa, pb, cpa, spa, crlo(SIZE(rlon))
274  INTEGER :: i, j
275!-------------------------------------------------------------------------------
276  IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF
277  crlo(:)=COS(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
278  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
279  cpa=COS(pa); spa=SIN(pa)
280  DO j=1,SIZE(rlat)
281    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
282    distance(:,j)=ACOS(cpa*COS(pb)+spa*SIN(pb)*crlo(:))
283  END DO
284
285END SUBROUTINE dist_sphe
286!
287!-------------------------------------------------------------------------------
288
289END MODULE grid_atob_m
290!
291!*******************************************************************************
292
Note: See TracBrowser for help on using the repository browser.