source: LMDZ6/branches/blowing_snow/libf/phylmd/grid_noro_m.F90 @ 5424

Last change on this file since 5424 was 3435, checked in by Laurent Fairhead, 6 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

  • as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

  • due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, YM, LF
File size: 23.5 KB
Line 
1!
2! $Id$
3!
4MODULE 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
14CONTAINS
15
16
17!-------------------------------------------------------------------------------
18!
19SUBROUTINE 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
298END SUBROUTINE grid_noro
299!
300!-------------------------------------------------------------------------------
301
302
303!-------------------------------------------------------------------------------
304!
305SUBROUTINE 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
425END SUBROUTINE grid_noro0
426!
427!-------------------------------------------------------------------------------
428
429
430!-------------------------------------------------------------------------------
431!
432SUBROUTINE 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
494CONTAINS
495
496
497SUBROUTINE 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,:)
503END SUBROUTINE get_fld
504
505SUBROUTINE 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
520END SUBROUTINE check_dim
521
522SUBROUTINE 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
533END SUBROUTINE ncerr
534
535END SUBROUTINE read_noro
536!
537!-------------------------------------------------------------------------------
538
539
540!-------------------------------------------------------------------------------
541!
542SUBROUTINE 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
572END SUBROUTINE MVA9
573!
574!-------------------------------------------------------------------------------
575
576
577END MODULE grid_noro_m
578
579
Note: See TracBrowser for help on using the repository browser.