source: LMDZ5/trunk/libf/phylmd/grid_noro_m.F90 @ 5444

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