1 |
|
|
!******************************************************************************* |
2 |
|
|
! |
3 |
|
|
MODULE 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 |
|
|
|
12 |
|
|
CONTAINS |
13 |
|
|
|
14 |
|
|
!------------------------------------------------------------------------------- |
15 |
|
|
! |
16 |
|
|
SUBROUTINE 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,m2, n1,n2, mx,my, nn, i,j |
32 |
|
|
LOGICAL :: li, lo, first=.TRUE. |
33 |
|
|
REAL :: inc, cpa, spa, crlo(SIZE(x_i)) |
34 |
|
|
REAL, SAVE :: pi, hpi |
35 |
|
|
INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot |
36 |
|
|
LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask |
37 |
|
|
REAL, DIMENSION(SIZE(x_i),SIZE(y_i)) :: dist |
38 |
|
|
REAL, DIMENSION(SIZE(x_o)) :: a, b |
39 |
|
|
REAL, DIMENSION(SIZE(y_o)) :: c, d |
40 |
|
|
REAL, PARAMETER :: thresh=1.E-5 |
41 |
|
|
!------------------------------------------------------------------------------- |
42 |
|
|
IF(first) THEN; pi=4.0*ATAN(1.0); hpi=pi/2.0; first=.FALSE.; END IF |
43 |
|
|
mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o) |
44 |
|
|
m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1) |
45 |
|
|
n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2) |
46 |
|
|
li=PRESENT(d_i ); IF(li) THEN; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF |
47 |
|
|
lo=PRESENT(d_o2); IF(lo) THEN; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF |
48 |
|
|
mi=assert_eq(mi,m1,TRIM(modname)//" mi") |
49 |
|
|
ni=assert_eq(ni,n1,TRIM(modname)//" ni") |
50 |
|
|
mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),TRIM(modname)//" mo") |
51 |
|
|
no=assert_eq(no,n2,my,SIZE(d_o1,2),TRIM(modname)//" no") |
52 |
|
|
mask(:,:)=.TRUE.; IF(PRESENT(msk)) mask(:,:)=msk(:,:) |
53 |
|
|
|
54 |
|
|
!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID |
55 |
|
|
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 |
56 |
|
|
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 |
57 |
|
|
a(1 )=x_o(1 )-(x_o(2 )-x_o(1 ))/2.0; a(2:mo )= b(1:mo-1) |
58 |
|
|
c(1 )=y_o(1 )-(y_o(2 )-y_o(1 ))/2.0; c(2:no )= d(1:no-1) |
59 |
|
|
|
60 |
|
|
!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID |
61 |
|
|
d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0 |
62 |
|
|
IF(lo) d_o2(:,:)=0. |
63 |
|
|
DO ji = 1, ni |
64 |
|
|
DO ii = 1, mi |
65 |
|
|
IF(li) inc=d_i(ii,ji) |
66 |
|
|
DO jo = 1, no |
67 |
|
|
IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND. & |
68 |
|
|
(y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) CYCLE |
69 |
|
|
DO io = 1, mo |
70 |
|
|
IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. & |
71 |
|
|
(x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) CYCLE |
72 |
|
|
num_tot(io,jo)=num_tot(io,jo)+1 |
73 |
|
|
IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc |
74 |
|
|
IF(.NOT.lo) CYCLE |
75 |
|
|
IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc |
76 |
|
|
END DO |
77 |
|
|
END DO |
78 |
|
|
END DO |
79 |
|
|
END DO |
80 |
|
|
|
81 |
|
|
!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL |
82 |
|
|
found(:,:)=num_tot(:,:)/=0 |
83 |
|
|
WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:)) |
84 |
|
|
IF(PRESENT(d_o2)) THEN |
85 |
|
|
WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:)) |
86 |
|
|
RETURN |
87 |
|
|
END IF |
88 |
|
|
nn=COUNT(found); IF(nn==0) RETURN |
89 |
|
|
|
90 |
|
|
!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2) |
91 |
|
|
DO io = 1, mo |
92 |
|
|
DO jo = 1, no |
93 |
|
|
IF(found(io,jo)) CYCLE |
94 |
|
|
! IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo |
95 |
|
|
crlo(:)=COS(x_o(io)-x_i(:)) !--- COS of points 1 and 2 angle |
96 |
|
|
cpa=COS(y_o(jo)); spa=SIN(y_o(jo)) |
97 |
|
|
DO j=1,ni; dist(:,j)=ACOS(spa*SIN(y_i(j))+cpa*COS(y_i(j))*crlo(:)); END DO |
98 |
|
|
nr=MINLOC(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr |
99 |
|
|
inc=1.0; IF(li) inc=d_i(nr(1),nr(2)) |
100 |
|
|
IF(mask(io,jo)) d_o1(io,jo)=inc |
101 |
|
|
END DO |
102 |
|
|
END DO |
103 |
|
|
|
104 |
|
|
END SUBROUTINE fine2coarse |
105 |
|
|
! |
106 |
|
|
!------------------------------------------------------------------------------- |
107 |
|
|
|
108 |
|
|
|
109 |
|
|
!------------------------------------------------------------------------------- |
110 |
|
|
! |
111 |
|
|
SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie) |
112 |
|
|
! |
113 |
|
|
!------------------------------------------------------------------------------- |
114 |
|
|
! Author: Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead) |
115 |
|
|
!------------------------------------------------------------------------------- |
116 |
|
|
! Purpose: Naive method to regrid on a coarser grid. |
117 |
|
|
! Value at a new point is an average of the old points lcoated in a zone |
118 |
|
|
! surrounding the considered new point. |
119 |
|
|
! No ponderation is used (see grille_p) |
120 |
|
|
! |
121 |
|
|
! (c) |
122 |
|
|
! ----d----- |
123 |
|
|
! | . . . .| |
124 |
|
|
! | | |
125 |
|
|
! (b)a . * . .b(a) |
126 |
|
|
! | | |
127 |
|
|
! | . . . .| |
128 |
|
|
! ----c----- |
129 |
|
|
! (d) |
130 |
|
|
! |
131 |
|
|
!------------------------------------------------------------------------------- |
132 |
|
|
IMPLICIT NONE |
133 |
|
|
!------------------------------------------------------------------------------- |
134 |
|
|
! Arguments: |
135 |
|
|
REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
136 |
|
|
REAL, INTENT(IN) :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT FIELD |
137 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
138 |
|
|
REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
139 |
|
|
!------------------------------------------------------------------------------- |
140 |
|
|
CALL fine2coarse(xdata,ydata,x,y,sortie,entree) |
141 |
|
|
|
142 |
|
|
END SUBROUTINE grille_m |
143 |
|
|
! |
144 |
|
|
!------------------------------------------------------------------------------- |
145 |
|
|
|
146 |
|
|
|
147 |
|
|
!------------------------------------------------------------------------------- |
148 |
|
|
! |
149 |
|
|
SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask) |
150 |
|
|
! |
151 |
|
|
!------------------------------------------------------------------------------- |
152 |
|
|
! Author: Z.X. Li (april 1st 1994) |
153 |
|
|
!------------------------------------------------------------------------------- |
154 |
|
|
! Purpose: Remap rugosity length ; constant value (0.001) on oceans. |
155 |
|
|
! Naive method (see grille_m) |
156 |
|
|
!------------------------------------------------------------------------------- |
157 |
|
|
IMPLICIT NONE |
158 |
|
|
!------------------------------------------------------------------------------- |
159 |
|
|
! Arguments: |
160 |
|
|
REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
161 |
|
|
REAL, INTENT(IN) :: entree(SIZE(xdata),SIZE(ydata)) !--- INPUT FIELD |
162 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
163 |
|
|
REAL, INTENT(OUT) :: sortie(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
164 |
|
|
REAL, INTENT(IN) :: mask (SIZE(x),SIZE(y)) !--- MASK |
165 |
|
|
!------------------------------------------------------------------------------- |
166 |
|
|
CALL fine2coarse(xdata,ydata,x,y,sortie,LOG(entree)) |
167 |
|
|
WHERE(NINT(mask)==1) |
168 |
|
|
sortie(:,:)=EXP(sortie(:,:)) |
169 |
|
|
ELSE WHERE |
170 |
|
|
sortie(:,:)=0.001 |
171 |
|
|
END WHERE |
172 |
|
|
|
173 |
|
|
END SUBROUTINE rugosite |
174 |
|
|
! |
175 |
|
|
!------------------------------------------------------------------------------- |
176 |
|
|
|
177 |
|
|
|
178 |
|
|
!------------------------------------------------------------------------------- |
179 |
|
|
! |
180 |
|
|
SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice) |
181 |
|
|
! |
182 |
|
|
!------------------------------------------------------------------------------- |
183 |
|
|
! Author: Z.X. Li (april 1st 1994) |
184 |
|
|
! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract. |
185 |
|
|
! field (between 0 and 1). |
186 |
|
|
! Naive method (see grille_m) |
187 |
|
|
!------------------------------------------------------------------------------- |
188 |
|
|
IMPLICIT NONE |
189 |
|
|
!------------------------------------------------------------------------------- |
190 |
|
|
! Arguments: |
191 |
|
|
REAL, INTENT(IN) :: xdata(:), ydata(:) !--- INPUT FIELD X AND Y COORD. |
192 |
|
|
REAL, INTENT(IN) :: glace01(:,:) !--- INPUT FIELD |
193 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT FIELD X AND Y COORD. |
194 |
|
|
REAL, INTENT(OUT) :: frac_ice(SIZE(x),SIZE(y)) !--- OUTPUT FIELD |
195 |
|
|
!------------------------------------------------------------------------------- |
196 |
|
|
CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=NINT(glace01)==1) |
197 |
|
|
|
198 |
|
|
END SUBROUTINE sea_ice |
199 |
|
|
! |
200 |
|
|
!------------------------------------------------------------------------------- |
201 |
|
|
|
202 |
|
|
|
203 |
|
|
!------------------------------------------------------------------------------- |
204 |
|
|
! |
205 |
|
|
SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs) |
206 |
|
|
! |
207 |
|
|
!------------------------------------------------------------------------------- |
208 |
|
|
! Author: Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995. |
209 |
|
|
! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell |
210 |
|
|
!------------------------------------------------------------------------------- |
211 |
|
|
IMPLICIT NONE |
212 |
|
|
!------------------------------------------------------------------------------- |
213 |
|
|
! Arguments: |
214 |
|
|
REAL, INTENT(IN) :: xrel(:), yrel(:) !--- INPUT FIELD X AND Y COORD. |
215 |
|
|
REAL, INTENT(IN) :: relief(:,:) !--- INPUT FIELD |
216 |
|
|
REAL, INTENT(IN) :: xmod(:), ymod(:) !--- OUTPUT FIELD X AND Y COORD. |
217 |
|
|
REAL, INTENT(OUT) :: rugs(SIZE(xmod),SIZE(ymod)) !--- OUTPUT FIELD |
218 |
|
|
!------------------------------------------------------------------------------- |
219 |
|
|
! Local variable: |
220 |
|
|
INTEGER :: k, nn |
221 |
|
|
INTEGER, PARAMETER:: itmp=360, jtmp=180 |
222 |
|
|
REAL :: out(SIZE(xmod),SIZE(xmod)), amin, amax |
223 |
|
|
REAL :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp) |
224 |
|
|
!------------------------------------------------------------------------------- |
225 |
|
|
|
226 |
|
|
!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION |
227 |
|
|
xtmp(:)=4.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)] |
228 |
|
|
ytmp(:)=2.0*ATAN(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)] |
229 |
|
|
CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp) |
230 |
|
|
cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2 |
231 |
|
|
nn=COUNT(cham2tmp<=-7.5) |
232 |
|
|
IF(nn/=0) THEN |
233 |
|
|
PRINT*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn |
234 |
|
|
CALL ABORT_GCM("", "", 1) |
235 |
|
|
END IF |
236 |
|
|
nn=COUNT(cham2tmp<0) |
237 |
|
|
IF(nn/=0) PRINT*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn |
238 |
|
|
WHERE(cham2tmp<0.0) cham2tmp=0.0 |
239 |
|
|
cham2tmp(:,:)=SQRT(cham2tmp(:,:)) |
240 |
|
|
amin=MINVAL(cham2tmp); amax=MAXVAL(cham2tmp) |
241 |
|
|
PRINT*, 'Ecart-type 1x1:', amin, amax |
242 |
|
|
|
243 |
|
|
!--- COMPUTE RUGOSITY AT REQUIRED SCALE |
244 |
|
|
WHERE(cham2tmp<0.001) cham2tmp=0.001 |
245 |
|
|
CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(LOG(cham2tmp))) |
246 |
|
|
out=EXP(out) |
247 |
|
|
amin=MINVAL(out); amax=MAXVAL(out) |
248 |
|
|
PRINT*, 'Ecart-type du modele:', amin, amax |
249 |
|
|
out=out/amax*20.0 |
250 |
|
|
amin=MINVAL(out); amax=MAXVAL(out) |
251 |
|
|
PRINT*, 'Longueur de rugosite du modele:', amin, amax |
252 |
|
|
rugs=REAL(out) |
253 |
|
|
|
254 |
|
|
END SUBROUTINE rugsoro |
255 |
|
|
! |
256 |
|
|
!------------------------------------------------------------------------------- |
257 |
|
|
|
258 |
|
|
END MODULE grid_atob_m |
259 |
|
|
! |
260 |
|
|
!******************************************************************************* |
261 |
|
|
|