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 lmdz_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 | |
---|