LMDZ
grid_noro_m.F90
Go to the documentation of this file.
1 MODULE grid_noro_m
2 !
3 !*******************************************************************************
4 
5  PRIVATE
6  PUBLIC :: grid_noro
7 
8 
9 CONTAINS
10 
11 
12 !-------------------------------------------------------------------------------
13 !
14 SUBROUTINE grid_noro(xd,yd,zd,x,y,zphi,zmea,zstd,zsig,zgam,zthe,zpic,zval,mask)
15 !
16 !-------------------------------------------------------------------------------
17 ! Author: F. Lott (see also Z.X. Li, A. Harzallah et L. Fairhead)
18 !-------------------------------------------------------------------------------
19 ! Purpose: Compute the Parameters of the SSO scheme as described in LOTT &MILLER
20 ! (1997) and LOTT(1999).
21 !-------------------------------------------------------------------------------
22 ! Comments:
23 ! * Target points are on a rectangular grid:
24 ! iim+1 latitudes including North and South Poles;
25 ! jjm+1 longitudes, with periodicity jjm+1=1.
26 ! * At the poles, the fields value is repeated jjm+1 time.
27 ! * The parameters a,b,c,d represent the limits of the target gridpoint region.
28 ! The means over this region are calculated from USN data, ponderated by a
29 ! weight proportional to the surface occupated by the data inside the model
30 ! gridpoint area. In most circumstances, this weight is the ratio between the
31 ! surfaces of the USN gridpoint area and the model gridpoint area.
32 !
33 ! (c)
34 ! ----d-----
35 ! | . . . .|
36 ! | |
37 ! (b)a . * . .b(a)
38 ! | |
39 ! | . . . .|
40 ! ----c-----
41 ! (d)
42 ! * Hard-coded US Navy dataset dimensions (imdp=2160 ; jmdp=1080) have been
43 ! removed (ALLOCATABLE used).
44 ! * iext (currently 10% of imdp) represents the margin to ensure output cells
45 ! on the edge are contained in input cells.
46 !===============================================================================
47  USE assert_eq_m, ONLY: assert_eq
48  USE print_control_mod, ONLY: lunout
49  IMPLICIT NONE
50  REAL, PARAMETER :: epsfra = 1.e-5
51 !-------------------------------------------------------------------------------
52 ! Arguments:
53  REAL, INTENT(IN) :: xd(:), yd(:) !--- INPUT COORDINATES (imdp) (jmdp)
54  REAL, INTENT(IN) :: zd(:,:) !--- INPUT FIELD (imdp,jmdp)
55  REAL, INTENT(IN) :: x(:), y(:) !--- OUTPUT COORDINATES (imar+1) (jmar)
56  REAL, INTENT(OUT) :: zphi(:,:) !--- GEOPOTENTIAL (imar+1,jmar)
57  REAL, INTENT(OUT) :: zmea(:,:) !--- MEAN OROGRAPHY (imar+1,jmar)
58  REAL, INTENT(OUT) :: zstd(:,:) !--- STANDARD DEVIATION (imar+1,jmar)
59  REAL, INTENT(OUT) :: zsig(:,:) !--- SLOPE (imar+1,jmar)
60  REAL, INTENT(OUT) :: zgam(:,:) !--- ANISOTROPY (imar+1,jmar)
61  REAL, INTENT(OUT) :: zthe(:,:) !--- SMALL AXIS ORIENTATION (imar+1,jmar)
62  REAL, INTENT(OUT) :: zpic(:,:) !--- MAXIMUM ALTITITUDE (imar+1,jmar)
63  REAL, INTENT(OUT) :: zval(:,:) !--- MINIMUM ALTITITUDE (imar+1,jmar)
64  REAL, INTENT(OUT) :: mask(:,:) !--- MASK (imar+1,jmar)
65 !-------------------------------------------------------------------------------
66 ! Local variables:
67  CHARACTER(LEN=256) :: modname="grid_noro"
68  REAL, ALLOCATABLE :: xusn(:), yusn(:) ! dim (imdp+2*iext) (jmdp+2)
69  REAL, ALLOCATABLE :: zusn(:,:) ! dim (imdp+2*iext,jmdp+2)
70 ! CORRELATIONS OF OROGRAPHY GRADIENT ! dim (imar+1,jmar)
71  REAL, ALLOCATABLE :: ztz(:,:), zxtzx(:,:), zytzy(:,:), zxtzy(:,:), weight(:,:)
72 ! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imar+2*iext,jmdp+2)
73  REAL, ALLOCATABLE :: zxtzxusn(:,:), zytzyusn(:,:), zxtzyusn(:,:)
74  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea0(:,:) ! dim (imar+1,jmar)
75  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
76  REAL, ALLOCATABLE :: a(:), b(:) ! dim (imax)
77  REAL, ALLOCATABLE :: c(:), d(:) ! dim (jmax)
78  LOGICAL :: masque_lu
79  INTEGER :: i, ii, imdp, imar, iext
80  INTEGER :: j, jj, jmdp, jmar, nn
81  REAL :: xpi, zdeltax, zlenx, weighx, xincr, zmeanor0
82  REAL :: rad, zdeltay, zleny, weighy, masque, zmeasud0
83  REAL :: zbordnor, zmeanor, zstdnor, zsignor, zweinor, zpicnor, zvalnor
84  REAL :: zbordsud, zmeasud, zstdsud, zsigsud, zweisud, zpicsud, zvalsud
85  REAL :: zbordest, zbordoue, xk, xl, xm, xp, xq, xw
86 !-------------------------------------------------------------------------------
87  imdp=assert_eq(SIZE(xd),SIZE(zd,1),trim(modname)//" imdp")
88  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),trim(modname)//" jmdp")
89  imar=assert_eq([SIZE(x),SIZE(zphi,1),SIZE(zmea,1),SIZE(zstd,1),SIZE(zsig,1), &
90  SIZE(zgam,1),SIZE(zthe,1),SIZE(zpic,1),SIZE(zval,1), &
91  SIZE(mask,1)],trim(modname)//" imar")-1
92  jmar=assert_eq([SIZE(y),SIZE(zphi,2),SIZE(zmea,2),SIZE(zstd,2),SIZE(zsig,2), &
93  SIZE(zgam,2),SIZE(zthe,2),SIZE(zpic,2),SIZE(zval,2), &
94  SIZE(mask,2)],trim(modname)//" jmar")
95 ! IF(imar/=iim) CALL abort_physic(TRIM(modname),'imar/=iim' ,1)
96 ! IF(jmar/=jjm+1) CALL abort_physic(TRIM(modname),'jmar/=jjm+1',1)
97  iext=imdp/10 !--- OK up to 36 degrees cell
98  xpi = acos(-1.)
99  rad = 6371229.
100  zdeltay=2.*xpi/REAL(jmdp)*rad
101  WRITE(lunout,*)"*** Orography parameters at sub-cell scale ***"
102 
103 !--- ARE WE USING A READ MASK ?
104  masque_lu=any(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
105  WRITE(lunout,*)'Masque lu: ',masque_lu
106 
107 !--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
108  ALLOCATE(xusn(imdp+2*iext))
109  xusn(1 +iext:imdp +iext)=xd(:)
110  xusn(1 : iext)=xd(1+imdp-iext:imdp)-2.*xpi
111  xusn(1+imdp+iext:imdp+2*iext)=xd(1 :iext)+2.*xpi
112 
113  ALLOCATE(yusn(jmdp+2))
114  yusn(1 )=yd(1) +(yd(1) -yd(2))
115  yusn(2:jmdp+1)=yd(:)
116  yusn( jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
117 
118  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
119  zusn(1 +iext:imdp +iext,2:jmdp+1)=zd(: , :)
120  zusn(1 : iext,2:jmdp+1)=zd(imdp-iext+1:imdp , :)
121  zusn(1+imdp +iext:imdp+2*iext,2:jmdp+1)=zd(1:iext , :)
122  zusn(1 :imdp/2+iext, 1)=zusn(1+imdp/2:imdp +iext, 2)
123  zusn(1+imdp/2+iext:imdp+2*iext, 1)=zusn(1 :imdp/2+iext, 2)
124  zusn(1 :imdp/2+iext, jmdp+2)=zusn(1+imdp/2:imdp +iext,jmdp+1)
125  zusn(1+imdp/2+iext:imdp+2*iext, jmdp+2)=zusn(1 :imdp/2+iext,jmdp+1)
126 
127 !--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
128  ALLOCATE(a(imar+1),b(imar+1))
129  b(1:imar)=(x(1:imar )+ x(2:imar+1))/2.0
130  b(imar+1)= x( imar+1)+(x( imar+1)-x(imar))/2.0
131  a(1)=x(1)-(x(2)-x(1))/2.0
132  a(2:imar+1)= b(1:imar)
133 
134  ALLOCATE(c(jmar),d(jmar))
135  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
136  d( jmar )= y( jmar )+(y( jmar)-y(jmar-1))/2.0
137  c(1)=y(1)-(y(2)-y(1))/2.0
138  c(2:jmar)=d(1:jmar-1)
139 
140 !--- INITIALIZATIONS:
141  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
142  ALLOCATE(zxtzx(imar+1,jmar)); zxtzx(:,:)= 0.0
143  ALLOCATE(zytzy(imar+1,jmar)); zytzy(:,:)= 0.0
144  ALLOCATE(zxtzy(imar+1,jmar)); zxtzy(:,:)= 0.0
145  ALLOCATE(ztz(imar+1,jmar)); ztz(:,:)= 0.0
146  zmea(:,:)= 0.0
147  zpic(:,:)=-1.e+10
148  zval(:,:)= 1.e+10
149 
150 !--- COMPUTE SLOPES CORRELATIONS ON USN GRID
151 ! CORRELATIONS OF USN OROGRAPHY GRADIENTS ! dim (imdp+2*iext,jmdp+2)
152  ALLOCATE(zytzyusn(imdp+2*iext,jmdp+2)); zytzyusn(:,:)=0.0
153  ALLOCATE(zxtzxusn(imdp+2*iext,jmdp+2)); zxtzxusn(:,:)=0.0
154  ALLOCATE(zxtzyusn(imdp+2*iext,jmdp+2)); zxtzyusn(:,:)=0.0
155  DO j = 2, jmdp+1
156  zdeltax=zdeltay*cos(yusn(j))
157  DO i = 2, imdp+2*iext-1
158  zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
159  zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
160  zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1)) /zdeltay &
161  & *(zusn(i+1,j)-zusn(i-1,j)) /zdeltax
162  END DO
163  END DO
164 
165 !--- SUMMATION OVER GRIDPOINT AREA
166  zleny=xpi/REAL(jmdp)*rad
167  xincr=xpi/REAL(jmdp)/2.
168  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
169  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
170  DO ii = 1, imar+1
171  DO jj = 1, jmar
172  DO j = 2,jmdp+1
173  zlenx =zleny *cos(yusn(j))
174  zdeltax=zdeltay*cos(yusn(j))
175  zbordnor=(xincr+c(jj)-yusn(j))*rad
176  zbordsud=(xincr-d(jj)+yusn(j))*rad
177  weighy=amax1(0.,amin1(zbordnor,zbordsud,zleny))
178  IF(weighy==0.) cycle
179  DO i = 2, imdp+2*iext-1
180  zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
181  zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
182  weighx=amax1(0.,amin1(zbordest,zbordoue,zlenx))
183  IF(weighx==0.) cycle
184  num_tot(ii,jj)=num_tot(ii,jj)+1.0
185  IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
186  weight(ii,jj)=weight(ii,jj)+weighx*weighy
187  zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
188  zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
189  zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
190  ztz(ii,jj)= ztz(ii,jj)+zusn(i,j)*zusn(i,j)*weighx*weighy
191  zmea(ii,jj)= zmea(ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
192  zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j)) !--- PEAKS
193  zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j)) !--- VALLEYS
194  END DO
195  END DO
196  END DO
197  END DO
198 
199 !--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
200  IF(.NOT.masque_lu) THEN
201  WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
202  END IF
203  nn=count(weight(:,1:jmar-1)==0.0)
204  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
205  WHERE(weight(:,:)/=0.0)
206  zmea(:,:)=zmea(:,:)/weight(:,:)
207  zxtzx(:,:)=zxtzx(:,:)/weight(:,:)
208  zytzy(:,:)=zytzy(:,:)/weight(:,:)
209  zxtzy(:,:)=zxtzy(:,:)/weight(:,:)
210  ztz(:,:)=ztz(:,:)/weight(:,:)
211  zstd(:,:)=ztz(:,:)-zmea(:,:)**2
212  END WHERE
213  WHERE(zstd(:,:)<0) zstd(:,:)=0.
214  zstd(:,:)=sqrt(zstd(:,:))
215 
216 !--- CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
217  zxtzx(:, 1)=zxtzx(:, 2)
218  zxtzx(:,jmar)=zxtzx(:,jmar-1)
219  zxtzy(:, 1)=zxtzy(:, 2)
220  zxtzy(:,jmar)=zxtzy(:,jmar-1)
221  zytzy(:, 1)=zytzy(:, 2)
222  zytzy(:,jmar)=zytzy(:,jmar-1)
223 
224 !=== FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
225 !--- FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
226 !-------------------------------------------------------------------------------
227  ALLOCATE(zmea0(imar+1,jmar))
228  zmea0(:,:)=zmea(:,:) ! GK211005 (CG) UNSMOOTHED TOPO
229  CALL mva9(zmea); CALL mva9(zstd); CALL mva9(zpic); CALL mva9(zval)
230  CALL mva9(zxtzx); CALL mva9(zxtzy); CALL mva9(zytzy)
231 
232 !--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD. (SURFACE PARAMS MEANINGLESS)
233  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
234  WHERE(mask>=0.1) mask_tmp = 1.
235  WHERE(weight(:,:)/=0.0)
236 ! zphi (:,:)= mask_tmp(:,:)*zmea (:,:) ! GK211005 (CG) not necessarly smoothed
237  zphi(:,:)= mask_tmp(:,:)*zmea0(:,:)
238  zmea0(:,:)= mask_tmp(:,:)*zmea0(:,:)
239  zmea(:,:)= mask_tmp(:,:)*zmea(:,:)
240  zpic(:,:)= mask_tmp(:,:)*zpic(:,:)
241  zval(:,:)= mask_tmp(:,:)*zval(:,:)
242  zstd(:,:)= mask_tmp(:,:)*zstd(:,:)
243  END WHERE
244  DO ii = 1, imar
245  DO jj = 1, jmar
246  IF (weight(ii,jj)/=0.0) THEN
247  !--- Coefficients K, L et M:
248  xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
249  xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
250  xm=zxtzy(ii,jj)
251  xp=xk-sqrt(xl**2+xm**2)
252  xq=xk+sqrt(xl**2+xm**2)
253  xw=1.e-8
254  IF(xp<=xw) xp=0.
255  IF(xq<=xw) xq=xw
256  IF(abs(xm)<=xw) xm=xw*sign(1.,xm)
257  !--- SLOPE
258  zsig(ii,jj)=sqrt(xq)*mask_tmp(ii,jj)
259  !---ISOTROPY
260  zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
261  !--- THETA ANGLE
262  zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask_tmp(ii,jj)
263  END IF
264  END DO
265  END DO
266  WRITE(lunout,*)' MEAN ORO:' ,maxval(zmea)
267  WRITE(lunout,*)' ST. DEV.:' ,maxval(zstd)
268  WRITE(lunout,*)' PENTE:' ,maxval(zsig)
269  WRITE(lunout,*)' ANISOTROP:',maxval(zgam)
270  WRITE(lunout,*)' ANGLE:' ,minval(zthe),maxval(zthe)
271  WRITE(lunout,*)' pic:' ,maxval(zpic)
272  WRITE(lunout,*)' val:' ,maxval(zval)
273 
274 !--- Values at poles
275  zmea0(imar+1,:)=zmea0(1,:)
276  zmea(imar+1,:)=zmea(1,:)
277  zphi(imar+1,:)=zphi(1,:)
278  zpic(imar+1,:)=zpic(1,:)
279  zval(imar+1,:)=zval(1,:)
280  zstd(imar+1,:)=zstd(1,:)
281  zsig(imar+1,:)=zsig(1,:)
282  zgam(imar+1,:)=zgam(1,:)
283  zthe(imar+1,:)=zthe(1,:)
284 
285  zweinor =sum(weight(1:imar, 1),dim=1)
286  zweisud =sum(weight(1:imar,jmar),dim=1)
287  zmeanor0=sum(weight(1:imar, 1)*zmea0(1:imar, 1),dim=1)
288  zmeasud0=sum(weight(1:imar,jmar)*zmea0(1:imar,jmar),dim=1)
289  zmeanor =sum(weight(1:imar, 1)*zmea(1:imar, 1),dim=1)
290  zmeasud =sum(weight(1:imar,jmar)*zmea(1:imar,jmar),dim=1)
291  zstdnor =sum(weight(1:imar, 1)*zstd(1:imar, 1),dim=1)
292  zstdsud =sum(weight(1:imar,jmar)*zstd(1:imar,jmar),dim=1)
293  zsignor =sum(weight(1:imar, 1)*zsig(1:imar, 1),dim=1)
294  zsigsud =sum(weight(1:imar,jmar)*zsig(1:imar,jmar),dim=1)
295  zpicnor =sum(weight(1:imar, 1)*zpic(1:imar, 1),dim=1)
296  zpicsud =sum(weight(1:imar,jmar)*zpic(1:imar,jmar),dim=1)
297  zvalnor =sum(weight(1:imar, 1)*zval(1:imar, 1),dim=1)
298  zvalsud =sum(weight(1:imar,jmar)*zval(1:imar,jmar),dim=1)
299 
300  zmea(:,1)=zmeanor /zweinor; zmea(:,jmar)=zmeasud /zweisud
301 ! zphi(:,1)=zmeanor0/zweinor; zphi(:,jmar)=zmeasud0/zweisud TO COMMIT
302  zphi(:,1)=zmeanor /zweinor; zphi(:,jmar)=zmeasud /zweisud
303  zpic(:,1)=zpicnor /zweinor; zpic(:,jmar)=zpicsud /zweisud
304  zval(:,1)=zvalnor /zweinor; zval(:,jmar)=zvalsud /zweisud
305  zstd(:,1)=zstdnor /zweinor; zstd(:,jmar)=zstdsud /zweisud
306  zsig(:,1)=zsignor /zweinor; zsig(:,jmar)=zsigsud /zweisud
307  zgam(:,1)=1.; zgam(:,jmar)=1.
308  zthe(:,1)=0.; zthe(:,jmar)=0.
309 
310 END SUBROUTINE grid_noro
311 !
312 !-------------------------------------------------------------------------------
313 
314 
315 !-------------------------------------------------------------------------------
316 !
317 SUBROUTINE mva9(x)
318 !
319 !-------------------------------------------------------------------------------
320  IMPLICIT NONE
321 ! MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
322 !-------------------------------------------------------------------------------
323 ! Arguments:
324  REAL, INTENT(INOUT) :: x(:,:)
325 !-------------------------------------------------------------------------------
326 ! Local variables:
327  REAL :: xf(size(x,dim=1),size(x,dim=2)), weightpb(-1:1,-1:1)
328  INTEGER :: i, j, imar, jmar
329 !-------------------------------------------------------------------------------
330  weightpb=reshape([((1./REAL((1+i**2)*(1+j**2)),i=-1,1),j=-1,1)],shape=[3,3])
331  weightpb=weightpb/sum(weightpb)
332  imar=SIZE(x,dim=1); jmar=SIZE(x,dim=2)
333  DO j=2,jmar-1
334  DO i=2,imar-1
335  xf(i,j)=sum(x(i-1:i+1,j-1:j+1)*weightpb(:,:))
336  END DO
337  END DO
338  DO j=2,jmar-1
339  xf(1,j)=sum(x(imar-1,j-1:j+1)*weightpb(-1,:))
340  xf(1,j)=xf(1,j)+sum(x(1:2,j-1:j+1)*weightpb(0:1,-1:1))
341  xf(imar,j)=xf(1,j)
342  END DO
343  xf(:, 1)=xf(:, 2)
344  xf(:,jmar)=xf(:,jmar-1)
345  x(:,:)=xf(:,:)
346 
347 END SUBROUTINE mva9
348 !
349 !-------------------------------------------------------------------------------
350 
351 
352 END MODULE grid_noro_m
353 
grid_noro_m::mva9
subroutine mva9(x)
Definition: grid_noro_m.F90:318
assert_eq_m
Definition: assert_eq_m.F90:2
grid_noro_m::grid_noro
subroutine, public grid_noro(xd, yd, zd, x, y, zphi, zmea, zstd, zsig, zgam, zthe, zpic, zval, mask)
Definition: grid_noro_m.F90:15
assert_eq_m::assert_eq
Definition: assert_eq_m.F90:6
i
!$Id klon initialisation mois suivants day_rain itap ENDIF!Calcul fin de nday_rain calcul nday_rain itap DO i
Definition: calcul_divers.h:24
print_control_mod
Definition: print_control_mod.F90:2
lunout
!$Header!gestion des impressions de sorties et de débogage la sortie standard prt_level COMMON comprint lunout
Definition: iniprint.h:7
grid_noro_m
Definition: grid_noro_m.F90:1