1 |
|
|
! |
2 |
|
|
! $Id$ |
3 |
|
|
! |
4 |
|
|
MODULE grid_noro_m |
5 |
|
|
! |
6 |
|
|
!******************************************************************************* |
7 |
|
|
|
8 |
|
|
USE print_control_mod, ONLY: lunout |
9 |
|
|
USE assert_eq_m, ONLY: assert_eq |
10 |
|
|
PRIVATE |
11 |
|
|
PUBLIC :: grid_noro, grid_noro0, read_noro |
12 |
|
|
|
13 |
|
|
|
14 |
|
|
CONTAINS |
15 |
|
|
|
16 |
|
|
|
17 |
|
|
!------------------------------------------------------------------------------- |
18 |
|
|
! |
19 |
|
|
SUBROUTINE grid_noro(xd,yd,zd,x,y,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask) |
20 |
|
|
! |
21 |
|
|
!------------------------------------------------------------------------------- |
22 |
|
|
! Author: F. Lott (see also Z.X. Li, A. Harzallah et L. Fairhead) |
23 |
|
|
!------------------------------------------------------------------------------- |
24 |
|
|
! Purpose: Compute the Parameters of the SSO scheme as described in LOTT &MILLER |
25 |
|
|
! (1997) and LOTT(1999). |
26 |
|
|
!------------------------------------------------------------------------------- |
27 |
|
|
! Comments: |
28 |
|
|
! * Target points are on a rectangular grid: |
29 |
|
|
! iim+1 latitudes including North and South Poles; |
30 |
|
|
! jjm+1 longitudes, with periodicity jjm+1=1. |
31 |
|
|
! * At the poles, the fields value is repeated jjm+1 time. |
32 |
|
|
! * The parameters a,b,c,d represent the limits of the target gridpoint region. |
33 |
|
|
! The means over this region are calculated from USN data, ponderated by a |
34 |
|
|
! weight proportional to the surface occupated by the data inside the model |
35 |
|
|
! gridpoint area. In most circumstances, this weight is the ratio between the |
36 |
|
|
! surfaces of the USN gridpoint area and the model gridpoint area. |
37 |
|
|
! |
38 |
|
|
! (c) |
39 |
|
|
! ----d----- |
40 |
|
|
! | . . . .| |
41 |
|
|
! | | |
42 |
|
|
! (b)a . * . .b(a) |
43 |
|
|
! | | |
44 |
|
|
! | . . . .| |
45 |
|
|
! ----c----- |
46 |
|
|
! (d) |
47 |
|
|
! * Hard-coded US Navy dataset dimensions (imdp=2160 ; jmdp=1080) have been |
48 |
|
|
! removed (ALLOCATABLE used). |
49 |
|
|
! * iext (currently 10% of imdp) represents the margin to ensure output cells |
50 |
|
|
! on the edge are contained in input cells. |
51 |
|
|
!=============================================================================== |
52 |
|
|
IMPLICIT NONE |
53 |
|
|
!------------------------------------------------------------------------------- |
54 |
|
|
! Arguments: |
55 |
|
|
REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) |
56 |
|
|
REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp) |
57 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) |
58 |
|
|
REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) |
59 |
|
|
REAL, INTENT(OUT) :: zmea(:,:) !--- MEAN OROGRAPHY (imar+1,jmar) |
60 |
|
|
REAL, INTENT(OUT) :: zstd(:,:) !--- STANDARD DEVIATION (imar+1,jmar) |
61 |
|
|
REAL, INTENT(OUT) :: zsig(:,:) !--- SLOPE (imar+1,jmar) |
62 |
|
|
REAL, INTENT(OUT) :: zgam(:,:) !--- ANISOTROPY (imar+1,jmar) |
63 |
|
|
REAL, INTENT(OUT) :: zthe(:,:) !--- SMALL AXIS ORIENTATION (imar+1,jmar) |
64 |
|
|
REAL, INTENT(OUT) :: zpic(:,:) !--- MAXIMUM ALTITITUDE (imar+1,jmar) |
65 |
|
|
REAL, INTENT(OUT) :: zval(:,:) !--- MINIMUM ALTITITUDE (imar+1,jmar) |
66 |
|
|
REAL, INTENT(OUT) :: mask(:,:) !--- MASK (imar+1,jmar) |
67 |
|
|
!------------------------------------------------------------------------------- |
68 |
|
|
! Local variables: |
69 |
|
|
CHARACTER(LEN=256) :: modname="grid_noro" |
70 |
|
|
REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) |
71 |
|
|
REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2) |
72 |
|
|
! CORRELATIONS OF OROGRAPHY GRADIENT ! dim (imar+1,jmar) |
73 |
|
|
REAL, ALLOCATABLE :: ztz(:,:), zxtzx(:,:), zytzy(:,:), zxtzy(:,:), weight(:,:) |
74 |
|
|
! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imar+2*iext,jmdp+2) |
75 |
|
|
REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:) |
76 |
|
|
REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imar+1,jmar) |
77 |
|
|
REAL, ALLOCATABLE :: a(:), b(:) ! dim (imar+1) |
78 |
|
|
REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmar) |
79 |
|
|
LOGICAL :: masque_lu |
80 |
|
|
INTEGER :: i, ii, imdp, imar, iext |
81 |
|
|
INTEGER :: j, jj, jmdp, jmar, nn |
82 |
|
|
REAL :: xpi, zdeltax, zlenx, weighx, xincr, zweinor, xk, xl, xm |
83 |
|
|
REAL :: rad, zdeltay, zleny, weighy, masque, zweisud, xp, xq, xw |
84 |
|
|
|
85 |
|
|
|
86 |
|
|
|
87 |
|
|
!------------------------------------------------------------------------------- |
88 |
|
|
imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") |
89 |
|
|
jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") |
90 |
|
|
imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), & |
91 |
|
|
SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), & |
92 |
|
|
SIZE(mask,1)],TRIM(modname)//" imar")-1 |
93 |
|
|
jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), & |
94 |
|
|
SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), & |
95 |
|
|
SIZE(mask,2)],TRIM(modname)//" jmar") |
96 |
|
|
! IF(imar/=iim) CALL abort_physic(TRIM(modname),'imar/=iim' ,1) |
97 |
|
|
! IF(jmar/=jjm+1) CALL abort_physic(TRIM(modname),'jmar/=jjm+1',1) |
98 |
|
|
iext=imdp/10 !--- OK up to 36 degrees cell |
99 |
|
|
xpi = ACOS(-1.) |
100 |
|
|
rad = 6371229. |
101 |
|
|
zdeltay=2.*xpi/REAL(jmdp)*rad |
102 |
|
|
WRITE(lunout,*)"*** Orography parameters at sub-cell scale ***" |
103 |
|
|
|
104 |
|
|
!--- ARE WE USING A READ MASK ? |
105 |
|
|
masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 |
106 |
|
|
WRITE(lunout,*)'Masque lu: ',masque_lu |
107 |
|
|
|
108 |
|
|
!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: |
109 |
|
|
ALLOCATE(xusn(imdp+2*iext)) |
110 |
|
|
xusn(1 +iext:imdp +iext)=xd(:) |
111 |
|
|
xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi |
112 |
|
|
xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi |
113 |
|
|
|
114 |
|
|
ALLOCATE(yusn(jmdp+2)) |
115 |
|
|
yusn(1 )=yd(1) +(yd(1) -yd(2)) |
116 |
|
|
yusn(2:jmdp+1)=yd(:) |
117 |
|
|
yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) |
118 |
|
|
|
119 |
|
|
ALLOCATE(zusn(imdp+2*iext,jmdp+2)) |
120 |
|
|
zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) |
121 |
|
|
zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) |
122 |
|
|
zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) |
123 |
|
|
zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) |
124 |
|
|
zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) |
125 |
|
|
zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) |
126 |
|
|
zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) |
127 |
|
|
|
128 |
|
|
!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) |
129 |
|
|
ALLOCATE(a(imar+1),b(imar+1)) |
130 |
|
|
b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 |
131 |
|
|
b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 |
132 |
|
|
a(1)=x(1)-(x(2)-x(1))/2.0 |
133 |
|
|
a(2:imar+1)= b(1:imar) |
134 |
|
|
|
135 |
|
|
ALLOCATE(c(jmar),d(jmar)) |
136 |
|
|
d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 |
137 |
|
|
d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 |
138 |
|
|
c(1)=y(1)-(y(2)-y(1))/2.0 |
139 |
|
|
c(2:jmar)=d(1:jmar-1) |
140 |
|
|
|
141 |
|
|
!--- INITIALIZATIONS: |
142 |
|
|
ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0 |
143 |
|
|
ALLOCATE(zxtzx (imar+1,jmar)); zxtzx (:,:)= 0.0 |
144 |
|
|
ALLOCATE(zytzy (imar+1,jmar)); zytzy (:,:)= 0.0 |
145 |
|
|
ALLOCATE(zxtzy (imar+1,jmar)); zxtzy (:,:)= 0.0 |
146 |
|
|
ALLOCATE(ztz (imar+1,jmar)); ztz (:,:)= 0.0 |
147 |
|
|
zmea(:,:)= 0.0 |
148 |
|
|
zpic(:,:)=-1.E+10 |
149 |
|
|
zval(:,:)= 1.E+10 |
150 |
|
|
|
151 |
|
|
!--- COMPUTE SLOPES CORRELATIONS ON USN GRID |
152 |
|
|
! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imdp+2*iext,jmdp+2) |
153 |
|
|
ALLOCATE(zytzyusn(imdp+2*iext,jmdp+2)); zytzyusn(:,:)=0.0 |
154 |
|
|
ALLOCATE(zxtzxusn(imdp+2*iext,jmdp+2)); zxtzxusn(:,:)=0.0 |
155 |
|
|
ALLOCATE(zxtzyusn(imdp+2*iext,jmdp+2)); zxtzyusn(:,:)=0.0 |
156 |
|
|
DO j = 2, jmdp+1 |
157 |
|
|
zdeltax=zdeltay*cos(yusn(j)) |
158 |
|
|
DO i = 2, imdp+2*iext-1 |
159 |
|
|
zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2 |
160 |
|
|
zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2 |
161 |
|
|
zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1)) /zdeltay & |
162 |
|
|
& *(zusn(i+1,j)-zusn(i-1,j)) /zdeltax |
163 |
|
|
END DO |
164 |
|
|
END DO |
165 |
|
|
|
166 |
|
|
!--- SUMMATION OVER GRIDPOINT AREA |
167 |
|
|
zleny=xpi/REAL(jmdp)*rad |
168 |
|
|
xincr=xpi/REAL(jmdp)/2. |
169 |
|
|
ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. |
170 |
|
|
ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. |
171 |
|
|
DO ii = 1, imar+1 |
172 |
|
|
DO jj = 1, jmar |
173 |
|
|
DO j = 2,jmdp+1 |
174 |
|
|
zlenx=zleny*COS(yusn(j)) |
175 |
|
|
zdeltax=zdeltay*COS(yusn(j)) |
176 |
|
|
weighy=(xincr+AMIN1(c(jj)-yusn(j),yusn(j)-d(jj)))*rad |
177 |
|
|
weighy=AMAX1(0.,AMIN1(weighy,zleny)) |
178 |
|
|
|
179 |
|
|
IF(weighy==0.) CYCLE |
180 |
|
|
DO i = 2, imdp+2*iext-1 |
181 |
|
|
weighx=(xincr+AMIN1(xusn(i)-a(ii),b(ii)-xusn(i)))*rad*COS(yusn(j)) |
182 |
|
|
weighx=AMAX1(0.,AMIN1(weighx,zlenx)) |
183 |
|
|
|
184 |
|
|
IF(weighx==0.) CYCLE |
185 |
|
|
num_tot(ii,jj)=num_tot(ii,jj)+1.0 |
186 |
|
|
IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 |
187 |
|
|
weight(ii,jj)=weight(ii,jj)+weighx*weighy |
188 |
|
|
zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy |
189 |
|
|
zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy |
190 |
|
|
zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy |
191 |
|
|
ztz (ii,jj)= ztz(ii,jj)+zusn(i,j)*zusn(i,j)*weighx*weighy |
192 |
|
|
zmea (ii,jj)= zmea(ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN |
193 |
|
|
zpic (ii,jj)=AMAX1(zpic(ii,jj),zusn(i,j)) !--- PEAKS |
194 |
|
|
zval (ii,jj)=AMIN1(zval(ii,jj),zusn(i,j)) !--- VALLEYS |
195 |
|
|
END DO |
196 |
|
|
END DO |
197 |
|
|
END DO |
198 |
|
|
END DO |
199 |
|
|
|
200 |
|
|
!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME |
201 |
|
|
IF(.NOT.masque_lu) THEN |
202 |
|
|
WHERE(weight(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:) |
203 |
|
|
END IF |
204 |
|
|
nn=COUNT(weight(:,:)==0.0) |
205 |
|
|
IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn |
206 |
|
|
WHERE(weight(:,:)/=0.0) |
207 |
|
|
zmea (:,:)=zmea (:,:)/weight(:,:) |
208 |
|
|
zxtzx(:,:)=zxtzx(:,:)/weight(:,:) |
209 |
|
|
zytzy(:,:)=zytzy(:,:)/weight(:,:) |
210 |
|
|
zxtzy(:,:)=zxtzy(:,:)/weight(:,:) |
211 |
|
|
ztz (:,:)=ztz (:,:)/weight(:,:) |
212 |
|
|
zstd (:,:)=ztz (:,:)-zmea(:,:)**2 |
213 |
|
|
END WHERE |
214 |
|
|
WHERE(zstd(:,:)<0) zstd(:,:)=0. |
215 |
|
|
zstd (:,:)=SQRT(zstd(:,:)) |
216 |
|
|
|
217 |
|
|
!--- CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES: |
218 |
|
|
zxtzx(:, 1)=zxtzx(:, 2) |
219 |
|
|
zxtzx(:,jmar)=zxtzx(:,jmar-1) |
220 |
|
|
zxtzy(:, 1)=zxtzy(:, 2) |
221 |
|
|
zxtzy(:,jmar)=zxtzy(:,jmar-1) |
222 |
|
|
zytzy(:, 1)=zytzy(:, 2) |
223 |
|
|
zytzy(:,jmar)=zytzy(:,jmar-1) |
224 |
|
|
|
225 |
|
|
!=== FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME. |
226 |
|
|
!--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS. |
227 |
|
|
!------------------------------------------------------------------------------- |
228 |
|
|
zphi(:,:)=zmea(:,:) ! GK211005 (CG) UNSMOOTHED TOPO |
229 |
|
|
|
230 |
|
|
CALL MVA9(zmea); CALL MVA9(zstd); CALL MVA9(zpic); CALL MVA9(zval) |
231 |
|
|
CALL MVA9(zxtzx); CALL MVA9(zxtzy); CALL MVA9(zytzy) |
232 |
|
|
|
233 |
|
|
!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS) |
234 |
|
|
WHERE(weight(:,:)==0.0.OR.mask<0.1) |
235 |
|
|
zphi(:,:)=0.0; zmea(:,:)=0.0; zpic(:,:)=0.0; zval(:,:)=0.0; zstd(:,:)=0.0 |
236 |
|
|
END WHERE |
237 |
|
|
DO ii = 1, imar |
238 |
|
|
DO jj = 1, jmar |
239 |
|
|
IF(weight(ii,jj)==0.0) CYCLE |
240 |
|
|
!--- Coefficients K, L et M: |
241 |
|
|
xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2. |
242 |
|
|
xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2. |
243 |
|
|
xm=zxtzy(ii,jj) |
244 |
|
|
xp=xk-SQRT(xl**2+xm**2) |
245 |
|
|
xq=xk+SQRT(xl**2+xm**2) |
246 |
|
|
xw=1.e-8 |
247 |
|
|
IF(xp<=xw) xp=0. |
248 |
|
|
IF(xq<=xw) xq=xw |
249 |
|
|
IF(ABS(xm)<=xw) xm=xw*SIGN(1.,xm) |
250 |
|
|
!--- SLOPE, ANISOTROPY AND THETA ANGLE |
251 |
|
|
zsig(ii,jj)=SQRT(xq) |
252 |
|
|
zgam(ii,jj)=xp/xq |
253 |
|
|
zthe(ii,jj)=90.*ATAN2(xm,xl)/xpi |
254 |
|
|
END DO |
255 |
|
|
END DO |
256 |
|
|
WHERE(weight(:,:)==0.0.OR.mask<0.1) |
257 |
|
|
zsig(:,:)=0.0; zgam(:,:)=0.0; zthe(:,:)=0.0 |
258 |
|
|
END WHERE |
259 |
|
|
|
260 |
|
|
WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) |
261 |
|
|
WRITE(lunout,*)' ST. DEV.:' ,MAXVAL(zstd) |
262 |
|
|
WRITE(lunout,*)' PENTE:' ,MAXVAL(zsig) |
263 |
|
|
WRITE(lunout,*)' ANISOTROP:',MAXVAL(zgam) |
264 |
|
|
WRITE(lunout,*)' ANGLE:' ,MINVAL(zthe),MAXVAL(zthe) |
265 |
|
|
WRITE(lunout,*)' pic:' ,MAXVAL(zpic) |
266 |
|
|
WRITE(lunout,*)' val:' ,MAXVAL(zval) |
267 |
|
|
|
268 |
|
|
!--- Values at redundant longitude |
269 |
|
|
zmea(imar+1,:)=zmea(1,:) |
270 |
|
|
zphi(imar+1,:)=zphi(1,:) |
271 |
|
|
zpic(imar+1,:)=zpic(1,:) |
272 |
|
|
zval(imar+1,:)=zval(1,:) |
273 |
|
|
zstd(imar+1,:)=zstd(1,:) |
274 |
|
|
zsig(imar+1,:)=zsig(1,:) |
275 |
|
|
zgam(imar+1,:)=zgam(1,:) |
276 |
|
|
zthe(imar+1,:)=zthe(1,:) |
277 |
|
|
|
278 |
|
|
!--- Values at north pole |
279 |
|
|
zweinor =SUM(weight(1:imar,1)) |
280 |
|
|
zmea(:,1)=SUM(weight(1:imar,1)*zmea(1:imar,1))/zweinor |
281 |
|
|
zphi(:,1)=SUM(weight(1:imar,1)*zphi(1:imar,1))/zweinor |
282 |
|
|
zpic(:,1)=SUM(weight(1:imar,1)*zpic(1:imar,1))/zweinor |
283 |
|
|
zval(:,1)=SUM(weight(1:imar,1)*zval(1:imar,1))/zweinor |
284 |
|
|
zstd(:,1)=SUM(weight(1:imar,1)*zstd(1:imar,1))/zweinor |
285 |
|
|
zsig(:,1)=SUM(weight(1:imar,1)*zsig(1:imar,1))/zweinor |
286 |
|
|
zgam(:,1)=1.; zthe(:,1)=0. |
287 |
|
|
|
288 |
|
|
!--- Values at south pole |
289 |
|
|
zweisud =SUM(weight(1:imar,jmar),DIM=1) |
290 |
|
|
zmea(:,jmar)=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar))/zweisud |
291 |
|
|
zphi(:,jmar)=SUM(weight(1:imar,jmar)*zphi(1:imar,jmar))/zweisud |
292 |
|
|
zpic(:,jmar)=SUM(weight(1:imar,jmar)*zpic(1:imar,jmar))/zweisud |
293 |
|
|
zval(:,jmar)=SUM(weight(1:imar,jmar)*zval(1:imar,jmar))/zweisud |
294 |
|
|
zstd(:,jmar)=SUM(weight(1:imar,jmar)*zstd(1:imar,jmar))/zweisud |
295 |
|
|
zsig(:,jmar)=SUM(weight(1:imar,jmar)*zsig(1:imar,jmar))/zweisud |
296 |
|
|
zgam(:,jmar)=1.; zthe(:,jmar)=0. |
297 |
|
|
|
298 |
|
|
END SUBROUTINE grid_noro |
299 |
|
|
! |
300 |
|
|
!------------------------------------------------------------------------------- |
301 |
|
|
|
302 |
|
|
|
303 |
|
|
!------------------------------------------------------------------------------- |
304 |
|
|
! |
305 |
|
|
SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask) |
306 |
|
|
! |
307 |
|
|
!=============================================================================== |
308 |
|
|
! Purpose: Extracted from grid_noro to provide geopotential height for dynamics |
309 |
|
|
! without any call to physics subroutines. |
310 |
|
|
!=============================================================================== |
311 |
|
|
IMPLICIT NONE |
312 |
|
|
!------------------------------------------------------------------------------- |
313 |
|
|
! Arguments: |
314 |
|
|
REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp) |
315 |
|
|
REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp, jmdp) |
316 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) |
317 |
|
|
REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) |
318 |
|
|
REAL, INTENT(OUT) :: mask(:,:) !--- MASK (imar+1,jmar) |
319 |
|
|
!------------------------------------------------------------------------------- |
320 |
|
|
! Local variables: |
321 |
|
|
CHARACTER(LEN=256) :: modname="grid_noro0" |
322 |
|
|
REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2) |
323 |
|
|
REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext, jmdp+2) |
324 |
|
|
REAL, ALLOCATABLE :: weight(:,:) ! dim (imar+1,jmar) |
325 |
|
|
REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imar+1,jmar) |
326 |
|
|
REAL, ALLOCATABLE :: a(:), b(:) ! dim (imar+1) |
327 |
|
|
REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmar) |
328 |
|
|
|
329 |
|
|
LOGICAL :: masque_lu |
330 |
|
|
INTEGER :: i, ii, imdp, imar, iext |
331 |
|
|
INTEGER :: j, jj, jmdp, jmar, nn |
332 |
|
|
REAL :: xpi, zlenx, zleny, weighx, weighy, xincr, masque, rad |
333 |
|
|
|
334 |
|
|
!------------------------------------------------------------------------------- |
335 |
|
|
imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp") |
336 |
|
|
jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp") |
337 |
|
|
imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1 |
338 |
|
|
jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar") |
339 |
|
|
iext=imdp/10 |
340 |
|
|
xpi = ACOS(-1.) |
341 |
|
|
rad = 6371229. |
342 |
|
|
|
343 |
|
|
!--- ARE WE USING A READ MASK ? |
344 |
|
|
masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 |
345 |
|
|
WRITE(lunout,*)'Masque lu: ',masque_lu |
346 |
|
|
|
347 |
|
|
!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES: |
348 |
|
|
ALLOCATE(xusn(imdp+2*iext)) |
349 |
|
|
xusn(1 +iext:imdp +iext)=xd(:) |
350 |
|
|
xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi |
351 |
|
|
xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi |
352 |
|
|
|
353 |
|
|
ALLOCATE(yusn(jmdp+2)) |
354 |
|
|
yusn(1 )=yd(1) +(yd(1) -yd(2)) |
355 |
|
|
yusn(2:jmdp+1)=yd(:) |
356 |
|
|
yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1)) |
357 |
|
|
|
358 |
|
|
ALLOCATE(zusn(imdp+2*iext,jmdp+2)) |
359 |
|
|
zusn(1 +iext:imdp +iext,2:jmdp+1)=zd (: , :) |
360 |
|
|
zusn(1 : iext,2:jmdp+1)=zd (imdp-iext+1:imdp , :) |
361 |
|
|
zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd (1:iext , :) |
362 |
|
|
zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2) |
363 |
|
|
zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2) |
364 |
|
|
zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1) |
365 |
|
|
zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1) |
366 |
|
|
|
367 |
|
|
!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID) |
368 |
|
|
ALLOCATE(a(imar+1),b(imar+1)) |
369 |
|
|
b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0 |
370 |
|
|
b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0 |
371 |
|
|
a(1)=x(1)-(x(2)-x(1))/2.0 |
372 |
|
|
a(2:imar+1)= b(1:imar) |
373 |
|
|
|
374 |
|
|
ALLOCATE(c(jmar),d(jmar)) |
375 |
|
|
d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0 |
376 |
|
|
d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0 |
377 |
|
|
c(1)=y(1)-(y(2)-y(1))/2.0 |
378 |
|
|
c(2:jmar)=d(1:jmar-1) |
379 |
|
|
|
380 |
|
|
!--- INITIALIZATIONS: |
381 |
|
|
ALLOCATE(weight(imar+1,jmar)); weight(:,:)=0.0; zphi(:,:)=0.0 |
382 |
|
|
|
383 |
|
|
!--- SUMMATION OVER GRIDPOINT AREA |
384 |
|
|
zleny=xpi/REAL(jmdp)*rad |
385 |
|
|
xincr=xpi/REAL(jmdp)/2. |
386 |
|
|
ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0. |
387 |
|
|
ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0. |
388 |
|
|
DO ii = 1, imar+1 |
389 |
|
|
DO jj = 1, jmar |
390 |
|
|
DO j = 2,jmdp+1 |
391 |
|
|
zlenx=zleny*COS(yusn(j)) |
392 |
|
|
weighy=(xincr+AMIN1(c(jj)-yusn(j),yusn(j)-d(jj)))*rad |
393 |
|
|
weighy=AMAX1(0.,AMIN1(weighy,zleny)) |
394 |
|
|
IF(weighy/=0) CYCLE |
395 |
|
|
DO i = 2, imdp+2*iext-1 |
396 |
|
|
weighx=(xincr+AMIN1(xusn(i)-a(ii),b(ii)-xusn(i)))*rad*COS(yusn(j)) |
397 |
|
|
weighx=AMAX1(0.,AMIN1(weighx,zlenx)) |
398 |
|
|
IF(weighx/=0) CYCLE |
399 |
|
|
num_tot(ii,jj)=num_tot(ii,jj)+1.0 |
400 |
|
|
IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0 |
401 |
|
|
weight(ii,jj)=weight(ii,jj)+weighx*weighy |
402 |
|
|
zphi (ii,jj)=zphi (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN |
403 |
|
|
END DO |
404 |
|
|
END DO |
405 |
|
|
END DO |
406 |
|
|
END DO |
407 |
|
|
|
408 |
|
|
!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME |
409 |
|
|
IF(.NOT.masque_lu) THEN |
410 |
|
|
WHERE(weight(:,:)/=0.0) mask=num_lan(:,:)/num_tot(:,:) |
411 |
|
|
END IF |
412 |
|
|
nn=COUNT(weight(:,:)==0.0) |
413 |
|
|
IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn |
414 |
|
|
WHERE(weight/=0.0) zphi(:,:)=zphi(:,:)/weight(:,:) |
415 |
|
|
|
416 |
|
|
!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS) |
417 |
|
|
WHERE(weight(:,:)==0.0.OR.mask<0.1) zphi(:,:)=0.0 |
418 |
|
|
WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zphi) |
419 |
|
|
|
420 |
|
|
!--- Values at redundant longitude and at poles |
421 |
|
|
zphi(imar+1,:)=zphi(1,:) |
422 |
|
|
zphi(:, 1)=SUM(weight(1:imar, 1)*zphi(1:imar, 1))/SUM(weight(1:imar, 1)) |
423 |
|
|
zphi(:,jmar)=SUM(weight(1:imar,jmar)*zphi(1:imar,jmar))/SUM(weight(1:imar,jmar)) |
424 |
|
|
|
425 |
|
|
END SUBROUTINE grid_noro0 |
426 |
|
|
! |
427 |
|
|
!------------------------------------------------------------------------------- |
428 |
|
|
|
429 |
|
|
|
430 |
|
|
!------------------------------------------------------------------------------- |
431 |
|
|
! |
432 |
|
|
SUBROUTINE read_noro(x,y,fname,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask) |
433 |
|
|
! |
434 |
|
|
!------------------------------------------------------------------------------- |
435 |
|
|
! Purpose: Read parameters usually determined with grid_noro from a file. |
436 |
|
|
!=============================================================================== |
437 |
|
|
USE netcdf, ONLY: NF90_OPEN, NF90_INQ_DIMID, NF90_INQUIRE_DIMENSION, & |
438 |
|
|
NF90_NOERR, NF90_CLOSE, NF90_INQ_VARID, NF90_GET_VAR, NF90_STRERROR, & |
439 |
|
|
NF90_NOWRITE |
440 |
|
|
IMPLICIT NONE |
441 |
|
|
!------------------------------------------------------------------------------- |
442 |
|
|
! Arguments: |
443 |
|
|
REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar) |
444 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: fname ! PARAMETERS FILE NAME |
445 |
|
|
REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar) |
446 |
|
|
REAL, INTENT(OUT) :: zmea(:,:) !--- MEAN OROGRAPHY (imar+1,jmar) |
447 |
|
|
REAL, INTENT(OUT) :: zstd(:,:) !--- STANDARD DEVIATION (imar+1,jmar) |
448 |
|
|
REAL, INTENT(OUT) :: zsig(:,:) !--- SLOPE (imar+1,jmar) |
449 |
|
|
REAL, INTENT(OUT) :: zgam(:,:) !--- ANISOTROPY (imar+1,jmar) |
450 |
|
|
REAL, INTENT(OUT) :: zthe(:,:) !--- SMALL AXIS ORIENTATION (imar+1,jmar) |
451 |
|
|
REAL, INTENT(OUT) :: zpic(:,:) !--- MAXIMUM ALTITUDE (imar+1,jmar) |
452 |
|
|
REAL, INTENT(OUT) :: zval(:,:) !--- MINIMUM ALTITUDE (imar+1,jmar) |
453 |
|
|
REAL, INTENT(OUT) :: mask(:,:) !--- MASK (imar+1,jmar) |
454 |
|
|
!------------------------------------------------------------------------------- |
455 |
|
|
! Local variables: |
456 |
|
|
CHARACTER(LEN=256) :: modname="read_noro" |
457 |
|
|
INTEGER :: imar, jmar, fid, did, vid |
458 |
|
|
LOGICAL :: masque_lu |
459 |
|
|
REAL :: xpi, d2r |
460 |
|
|
!------------------------------------------------------------------------------- |
461 |
|
|
imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), & |
462 |
|
|
SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), & |
463 |
|
|
SIZE(mask,1)],TRIM(modname)//" imar")-1 |
464 |
|
|
jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), & |
465 |
|
|
SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), & |
466 |
|
|
SIZE(mask,2)],TRIM(modname)//" jmar") |
467 |
|
|
xpi=ACOS(-1.0); d2r=xpi/180. |
468 |
|
|
WRITE(lunout,*)"*** Orography parameters at sub-cell scale from file ***" |
469 |
|
|
|
470 |
|
|
!--- ARE WE USING A READ MASK ? |
471 |
|
|
masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0 |
472 |
|
|
WRITE(lunout,*)'Masque lu: ',masque_lu |
473 |
|
|
CALL ncerr(NF90_OPEN(fname,NF90_NOWRITE,fid)) |
474 |
|
|
CALL check_dim('x','longitude',x(1:imar)) |
475 |
|
|
CALL check_dim('y','latitude' ,y(1:jmar)) |
476 |
|
|
IF(.NOT.masque_lu) CALL get_fld('mask',mask) |
477 |
|
|
CALL get_fld('Zphi',zphi) |
478 |
|
|
CALL get_fld('Zmea',zmea) |
479 |
|
|
CALL get_fld('mu' ,zstd) |
480 |
|
|
CALL get_fld('Zsig',zsig) |
481 |
|
|
CALL get_fld('Zgam',zgam) |
482 |
|
|
CALL get_fld('Zthe',zthe) |
483 |
|
|
zpic=zmea+2*zstd |
484 |
|
|
zval=MAX(0.,zmea-2.*zstd) |
485 |
|
|
CALL ncerr(NF90_CLOSE(fid)) |
486 |
|
|
WRITE(lunout,*)' MEAN ORO:' ,MAXVAL(zmea) |
487 |
|
|
WRITE(lunout,*)' ST. DEV.:' ,MAXVAL(zstd) |
488 |
|
|
WRITE(lunout,*)' PENTE:' ,MAXVAL(zsig) |
489 |
|
|
WRITE(lunout,*)' ANISOTROP:',MAXVAL(zgam) |
490 |
|
|
WRITE(lunout,*)' ANGLE:' ,MINVAL(zthe),MAXVAL(zthe) |
491 |
|
|
WRITE(lunout,*)' pic:' ,MAXVAL(zpic) |
492 |
|
|
WRITE(lunout,*)' val:' ,MAXVAL(zval) |
493 |
|
|
|
494 |
|
|
CONTAINS |
495 |
|
|
|
496 |
|
|
|
497 |
|
|
SUBROUTINE get_fld(var,fld) |
498 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: var |
499 |
|
|
REAL, INTENT(INOUT) :: fld(:,:) |
500 |
|
|
CALL ncerr(NF90_INQ_VARID(fid,var,vid),var) |
501 |
|
|
CALL ncerr(NF90_GET_VAR(fid,vid,fld(1:imar,:)),var) |
502 |
|
|
fld(imar+1,:)=fld(1,:) |
503 |
|
|
END SUBROUTINE get_fld |
504 |
|
|
|
505 |
|
|
SUBROUTINE check_dim(dimd,nam,dimv) |
506 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: dimd |
507 |
|
|
CHARACTER(LEN=*), INTENT(IN) :: nam |
508 |
|
|
REAL, INTENT(IN) :: dimv(:) |
509 |
|
|
REAL, ALLOCATABLE :: tmp(:) |
510 |
|
|
INTEGER :: n |
511 |
|
|
CALL ncerr(NF90_INQ_DIMID(fid,dimd,did)) |
512 |
|
|
CALL ncerr(NF90_INQUIRE_DIMENSION(fid,did,len=n)); ALLOCATE(tmp(n)) |
513 |
|
|
CALL ncerr(NF90_INQ_VARID(fid,dimd,did)) |
514 |
|
|
CALL ncerr(NF90_GET_VAR(fid,did,tmp)) |
515 |
|
|
IF(MAXVAL(tmp)>xpi) tmp=tmp*d2r |
516 |
|
|
IF(n/=SIZE(dimv).OR.ANY(ABS(tmp-dimv)>1E-6)) THEN |
517 |
|
|
WRITE(lunout,*)'Problem with file "'//TRIM(fname)//'".' |
518 |
|
|
CALL abort_physic(modname,'Grid differs from LMDZ for '//TRIM(nam)//'.',1) |
519 |
|
|
END IF |
520 |
|
|
END SUBROUTINE check_dim |
521 |
|
|
|
522 |
|
|
SUBROUTINE ncerr(ncres,var) |
523 |
|
|
IMPLICIT NONE |
524 |
|
|
INTEGER, INTENT(IN) :: ncres |
525 |
|
|
CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: var |
526 |
|
|
CHARACTER(LEN=256) :: mess |
527 |
|
|
IF(ncres/=NF90_NOERR) THEN |
528 |
|
|
mess='Problem with file "'//TRIM(fname)//'"' |
529 |
|
|
IF(PRESENT(var)) mess=TRIM(mess)//' and variable "'//TRIM(var)//'"' |
530 |
|
|
WRITE(lunout,*)TRIM(mess)//'.' |
531 |
|
|
CALL abort_physic(modname,NF90_STRERROR(ncres),1) |
532 |
|
|
END IF |
533 |
|
|
END SUBROUTINE ncerr |
534 |
|
|
|
535 |
|
|
END SUBROUTINE read_noro |
536 |
|
|
! |
537 |
|
|
!------------------------------------------------------------------------------- |
538 |
|
|
|
539 |
|
|
|
540 |
|
|
!------------------------------------------------------------------------------- |
541 |
|
|
! |
542 |
|
|
SUBROUTINE MVA9(x) |
543 |
|
|
! |
544 |
|
|
!------------------------------------------------------------------------------- |
545 |
|
|
IMPLICIT NONE |
546 |
|
|
! MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS |
547 |
|
|
!------------------------------------------------------------------------------- |
548 |
|
|
! Arguments: |
549 |
|
|
REAL, INTENT(INOUT) :: x(:,:) |
550 |
|
|
!------------------------------------------------------------------------------- |
551 |
|
|
! Local variables: |
552 |
|
|
REAL :: xf(SIZE(x,DIM=1),SIZE(x,DIM=2)), WEIGHTpb(-1:1,-1:1) |
553 |
|
|
INTEGER :: i, j, imar, jmar |
554 |
|
|
!------------------------------------------------------------------------------- |
555 |
|
|
WEIGHTpb=RESHAPE([((1./REAL((1+i**2)*(1+j**2)),i=-1,1),j=-1,1)],SHAPE=[3,3]) |
556 |
|
|
WEIGHTpb=WEIGHTpb/SUM(WEIGHTpb) |
557 |
|
|
imar=SIZE(X,DIM=1); jmar=SIZE(X,DIM=2) |
558 |
|
|
DO j=2,jmar-1 |
559 |
|
|
DO i=2,imar-1 |
560 |
|
|
xf(i,j)=SUM(x(i-1:i+1,j-1:j+1)*WEIGHTpb(:,:)) |
561 |
|
|
END DO |
562 |
|
|
END DO |
563 |
|
|
DO j=2,jmar-1 |
564 |
|
|
xf(1,j)=SUM(x(imar-1,j-1:j+1)*WEIGHTpb(-1,:)) |
565 |
|
|
xf(1,j)=xf(1,j)+SUM(x(1:2,j-1:j+1)*WEIGHTpb(0:1,-1:1)) |
566 |
|
|
xf(imar,j)=xf(1,j) |
567 |
|
|
END DO |
568 |
|
|
xf(:, 1)=xf(:, 2) |
569 |
|
|
xf(:,jmar)=xf(:,jmar-1) |
570 |
|
|
x(:,:)=xf(:,:) |
571 |
|
|
|
572 |
|
|
END SUBROUTINE MVA9 |
573 |
|
|
! |
574 |
|
|
!------------------------------------------------------------------------------- |
575 |
|
|
|
576 |
|
|
|
577 |
|
|
END MODULE grid_noro_m |
578 |
|
|
|
579 |
|
|
|