source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obsaverg_h2d.h90 @ 11202

Last change on this file since 11202 was 11202, checked in by jcastill, 15 months ago

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

File size: 34.6 KB
Line 
1   !!----------------------------------------------------------------------
2   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
3   !! $Id$
4   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
5   !!----------------------------------------------------------------------
6
7   SUBROUTINE obs_avg_h2d_init( kpk, kpk2, kmaxifp, kmaxjfp, k2dint, plam,  pphi, &
8      &                         pglam, pgphi, pglamf, pgphif, pmask, plamscl, pphiscl, lindegrees, &
9      &                         pweig, pobsmask, iminpoints )
10      !!-----------------------------------------------------------------------
11      !!
12      !!                     ***  ROUTINE obs_avg_h2d_init  ***
13      !!
14      !! ** Purpose : Computes weights for horizontal averaging to the
15      !!              observation footprint.
16      !!
17      !! ** Method  : Horizontal averaging to the observation footprint using
18      !!              model values at a defined area.
19      !!
20      !!    Averaging schemes :
21      !!
22      !!    Two horizontal averaging schemes are available:
23      !!        - weighted radial footprint        (k2dint = 5)
24      !!        - weighted rectangular footprint   (k2dint = 6)
25      !!
26      !! History :
27      !!        ! 13-10 (M. Martin)
28      !!-----------------------------------------------------------------------
29      !! * Modules used
30      !! * Arguments
31      INTEGER, INTENT(IN) :: &
32         & kpk,   &             ! Parameter values for automatic arrays
33         & kpk2,  &
34         & kmaxifp,  &          ! Max size of model points in i-direction for obs footprint
35         & kmaxjfp,  &          ! Max size of model points in j-direction for obs footprint
36         & k2dint               ! Averaging scheme options - see header
37      REAL(KIND=wp), INTENT(INOUT) :: &
38         & plam, &              ! Geographical (lat,lon) coordinates of
39         & pphi                 ! observation
40      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
41         & pglam, &             ! Model variable lon
42         & pgphi                ! Model variable lat
43      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
44         & pglamf, &            ! Model variable lon at corners of grid-boxes
45         & pgphif               ! Model variable lat at corners of grid-boxes
46      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
47         & pmask                ! Model variable mask
48      REAL(KIND=wp), INTENT(IN) :: &
49         & plamscl, &           ! Diameter (lat,lon) of obs footprint in metres
50         & pphiscl              ! This is the full width (rather than half-width)
51      LOGICAL, INTENT(IN) :: &
52         & lindegrees           ! T=> obs footprint specified in degrees, F=> in metres
53      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
54         & pweig                ! Weights for averaging
55      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
56         & pobsmask             ! Vertical mask for observations
57      INTEGER, INTENT(IN), OPTIONAL :: &
58         & iminpoints           ! Reject point which is not surrounded
59                                ! by at least iminpoints sea points
60
61      !! * Local declarations
62      INTEGER :: &
63         & jk
64      INTEGER :: &
65         & ikmax
66
67
68      !------------------------------------------------------------------------
69      !
70      !------------------------------------------------------------------------
71
72      !------------------------------------------------------------------------
73      ! Initialize number of levels
74      !------------------------------------------------------------------------
75      IF ( kpk2 == 1 ) THEN
76         ikmax = 1
77      ELSEIF ( kpk2 == kpk) THEN
78         ikmax = kpk-1
79      ENDIF
80
81
82         SELECT CASE (k2dint)
83         CASE(5)
84            CALL obs_avg_h2d_rad( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
85      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
86         CASE(6)
87            CALL obs_avg_h2d_rec( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
88      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
89         END SELECT
90
91
92      END SUBROUTINE obs_avg_h2d_init
93
94
95      SUBROUTINE obs_avg_h2d_rad( kpk2, kmax, kmaxifp, kmaxjfp, plam, pphi, &
96      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
97      !!-----------------------------------------------------------------------
98      !!
99      !!                     ***  ROUTINE obs_avg_h2d_rad  ***
100      !!
101      !! ** Purpose : Computes weights for horizontal averaging to the
102      !!              observation using a radial footprint.
103      !!
104      !! ** Method  : Calculate whether each grid box is completely or
105      !!              partially within the observation footprint.
106      !!              If it is partially in the footprint then calculate
107      !!              the ratio of the area inside the footprint to the total
108      !!              area of the grid box.
109      !!
110      !! History :
111      !!        ! 14-01 (M. Martin)
112      !!-----------------------------------------------------------------------
113      !! * Modules used
114      USE phycst,   ONLY : &  ! Physical constants
115         & ra,  &
116         & rpi
117
118      !! * Arguments
119      INTEGER, INTENT(IN) :: &
120         & kpk2, &             ! Parameter values for automatic arrays
121         & kmax
122
123      INTEGER, INTENT(IN) :: &
124         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
125         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
126
127      REAL(KIND=wp), INTENT(IN) :: &
128         & plam, &
129         & pphi                 ! Geographical (lat,lon) coordinates of
130                                ! observation
131      REAL(KIND=wp), INTENT(IN) :: &
132         & plamscl, &           ! Diameter (lat,lon) of obs footprint in metres or degrees (see below)
133         & pphiscl              ! This is the full width (rather than half-width)
134      LOGICAL, INTENT(IN) :: &
135         & lindegrees           ! T=>scales specified in degrees, F=> in metres
136      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
137         & pmask                ! Model variable mask
138      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
139         & pglam, &             ! Model variable lon
140         & pgphi                ! Model variable lat
141      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
142         & pglamf, &             ! Model variable lon at corners of grid boxes
143         & pgphif                ! Model variable lat at corners of grid boxes
144      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
145         & pweig                ! Weights for interpolation
146
147      !! Local declarations
148      INTEGER :: ji, jj, jk
149      INTEGER :: jvert, jis, jjs
150      INTEGER :: jnumvert, jnumvertbig
151      INTEGER, PARAMETER :: &
152         & jnumsubgrid = 20     ! The number of sub grid-boxes (in x and y directions) used to approximate area of obs fp
153
154      REAL(KIND=wp), DIMENSION(4) :: &
155         & zxvert, zyvert, &    ! The lon/lat of the vertices(corners) of the grid box in m relative to centre of obs fp
156         & zdist                ! Distance of each vertex to the centre of the obs footprint
157      REAL(KIND=wp), DIMENSION(4) :: &
158         & zxgrid, zygrid, &      ! Distance of each vertex of grid box to the centre of the grid box in x/y directions
159         & zdgrid
160      REAL(KIND=wp) :: &
161         & zdx, zdy, &          ! The sub grid-box sizes (in metres)
162         & zarea_subbox, &      ! The area of each sub grid-box (in metres squared)
163         & zxpos, zypos, &      ! The x,y position (relative to centre of obs footprint) of the centre of each sub grid-box
164         & zsubdist, &          ! The distance of the centre of each sub grid-box from the centre of the obs footprint
165         & zarea_fp, &          ! Total area of obs footprint within the grid box
166         & zareabox             ! Total area of the grid box
167      REAL(KIND=wp) :: &
168         & zphiscl_m, &         ! Diameter of obs footprint in metres
169         & zlamscl_m
170      !---------------------------------------------------------------------------------------------------
171      !Initialise weights to zero.
172      pweig(:,:,:) = 0.0_wp
173     
174      !Two footprint sizes can be specified in the namelist but this routine assumes a circular footprint.
175      !If the two sizes are different then write out a warning.
176      IF ( pphiscl /= plamscl ) THEN
177               CALL ctl_warn( 'obs_avg_h2d_rad:',   &
178                  &           'The two components of the obs footprint size are not equal', &
179                  &           'yet the radial option has been selected - using pphiscl here' )
180      ENDIF
181     
182      DO jk = 1, kmax
183         DO ji = 1, kmaxifp
184            DO jj = 1, kmaxjfp
185           
186               IF ( pmask(ji,jj,jk) == 1.0_wp ) THEN
187
188                  IF ( lindegrees ) THEN
189                     !If the scales are specified in degrees, work out the
190                     !scales (metres) in x/y directions
191                     CALL obs_deg2dist( 1, 1, pglam(ji,jj), pgphi(ji,jj), &
192                        &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
193                  ELSE
194                     zphiscl_m = pphiscl
195                  ENDIF
196
197
198                  ! Work out the area of the grid box using distance of corners relative to centre of grid box
199                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
200                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
201                     &                  pglam(ji,jj), pgphi(ji,jj), zxgrid, zygrid, zdgrid)
202                  zareabox = ABS( zxgrid(1) - zxgrid(2) ) * ABS( zygrid(1) - zygrid(4) )
203
204                  !1. Determine how many of the vertices of the grid box lie within the circle
205                 
206                  !For each vertex, calculate its location and distance relative
207                  !to the centre of the observation footprint
208                 
209                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
210                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
211                     &                  plam, pphi, zxvert, zyvert, zdist)
212
213                  jnumvert = 0
214                  jnumvertbig = 0
215                  DO jvert = 1, 4
216
217                     !If the distance to the center to the observation footprint is less
218                     !than the radius of the footprint (half the diameter) then this
219                     !vertex is within the observation footprint
220                     IF ( zdist(jvert) <= ( zphiscl_m / 2.0_wp ) ) jnumvert = jnumvert + 1
221                     
222                     !For expediency, check if the vertices are "nearly" within the obs
223                     !footprint as if none of them are close to the edge of the footprint
224                     !then the footprint is unlikely to be intersecting the grid box
225                     IF ( zdist(jvert) - ( 0.5_wp * zareabox ) <= ( zphiscl_m / 2.0 ) ) &
226                        & jnumvertbig = jnumvertbig + 1
227                     
228                  END DO
229                 
230                  !2. If none of the vertices are even close to the edge of the obs
231                  !footprint then leave weight as zero and cycle to next grid box.
232                  IF ( jnumvertbig == 0 ) CYCLE
233
234                  !3. If all the vertices of the box are within the observation footprint then the
235                  !      whole grid box is within the footprint so set the weight to one and
236                  !      move to the next grid box.
237                  IF ( jnumvert == 4 ) THEN
238                     pweig(ji,jj,jk) = 1.0_wp
239                     CYCLE
240                  ENDIF
241
242
243                  !4. Use a brute force technique for calculating the area within
244                  !   the grid box covered by the obs footprint.
245                  !  (alternative could be to use formulae on
246                  !         http://mathworld.wolfram.com/Circle-LineIntersection.html)
247                  !   For now split the grid box into a specified number of smaller
248                  !   boxes and add up the area of those whose centre is within the obs footprint.
249                  !   Order of vertices is 1=topleft, 2=topright, 3=bottomright, 4=bottomleft
250                  zdx = ABS( zxvert(3) - zxvert(4) ) / REAL(jnumsubgrid, wp)
251                  zdy = ABS( zyvert(1) - zyvert(4) ) / REAL(jnumsubgrid, wp)
252                  zarea_subbox = zdx * zdy
253
254                  zarea_fp = 0.0_wp
255                  DO jis = 1, jnumsubgrid
256                     zxpos = zxvert(4) + ( REAL(jis, wp) * zdx ) - (0.5_wp * zdx )
257                     DO jjs = 1, jnumsubgrid
258                        !Find the distance of the centre of this sub grid box to the
259                        !centre of the obs footprint
260                        zypos = zyvert(4) + ( REAL(jjs, wp) * zdy ) - ( 0.5_wp * zdy )
261                        zsubdist = SQRT( (zxpos * zxpos) + (zypos * zypos) )
262                        IF ( zsubdist < ( zphiscl_m / 2.0_wp ) ) &
263                           &  zarea_fp = zarea_fp + zarea_subbox
264                     END DO
265                  END DO
266
267                  !6. Calculate the ratio of the area of the footprint within the box
268                  !   to the total area of the grid box and use this fraction to weight
269                  !   the model value in this grid box.
270                  pweig(ji,jj,jk) = MIN( zarea_fp / zareabox, 1.0_wp )
271
272               END IF  !pmask
273            END DO
274         END DO
275      END DO
276     
277      END SUBROUTINE obs_avg_h2d_rad
278
279
280      SUBROUTINE obs_avg_h2d_rec( kpk2, kmax, kmaxifp, kmaxjfp, plam, pphi, &
281      &                           plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
282      !!-----------------------------------------------------------------------
283      !!
284      !!                     ***  ROUTINE obs_avg_h2d_rec  ***
285      !!
286      !! ** Purpose : Computes weights for horizontal averaging to the
287      !!              observation using a rectangular footprint which
288      !!              is aligned with lines of lat/lon.
289      !!
290      !! ** Method  : Horizontal averaging to the observation footprint using
291      !!              model values at a defined area.
292      !!
293      !! History :
294      !!        ! 14-01 (M. Martin)
295      !!-----------------------------------------------------------------------
296      !! * Modules used
297      USE phycst,   ONLY : &    ! Physical constants
298         & ra,  &
299         & rpi
300
301      !! * Arguments
302      INTEGER, INTENT(IN) :: &
303         & kpk2, &              ! Parameter values for automatic arrays
304         & kmax
305
306      INTEGER, INTENT(IN) :: &
307         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
308         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
309
310      REAL(KIND=wp), INTENT(IN) :: &
311         & plam, &
312         & pphi                 ! Geographical (lat,lon) coordinates of
313                                ! observation
314      REAL(KIND=wp), INTENT(IN) :: &
315         & plamscl, &
316         & pphiscl              ! Width in x/y directions of obs footprint in metres
317                                ! This is the full width (rather than half-width)
318      LOGICAL, INTENT(IN) :: &
319         & lindegrees           !T=> scales specified in degrees, F=> in metres
320      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
321         & pmask                ! Model variable mask
322      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
323         & pglam, &             ! Model variable lat at centre of grid boxes
324         & pgphi                ! Model variable lon at centre of grid boxes
325      REAL(KIND=wp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
326         & pglamf, &             ! Model variable lat at corners of grid boxes
327         & pgphif                ! Model variable lon at corners of grid boxes
328      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) ::  &
329         & pweig                ! Weights for interpolation
330
331      !! Local declarations
332      INTEGER :: ji, jj, jk
333      INTEGER :: jvert
334      INTEGER, DIMENSION(4) :: &
335         & jnumvert
336      REAL(KIND=wp), DIMENSION(4) :: &
337         & zxvert, zyvert         ! The lon/lat of the vertices(corners) of the grid box in m relative to centre of obs fp
338      REAL(KIND=wp), DIMENSION(4) :: &
339         & zdist                  ! Distance of each vertex to the centre of the obs footprint
340      REAL(KIND=wp), DIMENSION(4) :: &
341         & zxgrid, zygrid, &      ! Distance of each vertex of grid box to the centre of the grid box in x/y directions
342         & zdgrid
343      REAL(KIND=wp) :: &
344         & zareabox, &            ! Total area of grid box
345         & zarea_fp, &            ! Total area of obs footprint
346         & zarea_intersect        ! Area of the intersection between the grid box and the obs footprint
347      REAL(KIND=wp) :: &
348         & zlamscl_m, &
349         & zphiscl_m              ! Total width (lat,lon) of obs footprint in metres
350      REAL(KIND=wp) :: &
351         & z_awidth, z_aheight, & ! Width and height of model grid box
352         & z_cwidth, z_cheight    ! Width and height of union of model grid box and obs footprint
353      REAL(KIND=wp) :: &
354         & zleft, zright, &       ! Distance (metres) of corners area of intersection
355         & ztop, zbottom          ! between grid box and obs footprint
356
357      !-----------------------------------------------------------------------
358
359      !Initialise weights to zero
360      pweig(:,:,:) = 0.0_wp
361     
362      !Loop over the grid boxes which have been identified as potentially being within the
363      !observation footprint
364      DO jk = 1, kmax
365         DO ji = 1, kmaxifp
366            DO jj = 1, kmaxjfp
367
368               IF ( pmask(ji,jj,jk) == 1.0_wp ) THEN
369
370
371                  IF ( lindegrees ) THEN
372                     !If the scales are specified in degrees, work out the
373                     !scales (metres) in x/y directions
374                     CALL obs_deg2dist( 1, 1, pglam(ji,jj), pgphi(ji,jj), &
375                        &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
376                  ELSE
377                     zlamscl_m = plamscl
378                     zphiscl_m = pphiscl
379                  ENDIF
380
381                  ! Work out the area of the grid box using distance of corners relative to centre of grid box
382                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
383                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
384                     &                  pglam(ji,jj), pgphi(ji,jj), zxgrid, zygrid, zdgrid)
385
386                  !Calculate width and height of model grid box
387                  z_awidth  = ABS( zxgrid(1) - zxgrid(2) )
388                  z_aheight = ABS( zygrid(1) - zygrid(4) )
389                  zareabox = z_awidth * z_aheight
390
391                  ! Work out area of the observation footprint
392                  zarea_fp = zlamscl_m * zphiscl_m
393
394                  ! For each corner of the grid box, calculate its location and distance relative
395                  ! to the centre of the observation footprint
396                  CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
397                     &                  pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
398                     &                  plam, pphi, zxvert, zyvert, zdist)
399
400                  !Work out maximum width and height of rectangle covered by corners of obs fp and grid box
401                  z_cwidth  = MAX( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp ) - &
402                     &       MIN( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp )
403                     
404                  z_cheight = MAX( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp ) - &
405                     &       MIN( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp )
406                 
407                  IF ( ( z_cwidth  >= z_awidth  + zlamscl_m  ) .OR. &
408                     & ( z_cheight >= z_aheight + zphiscl_m ) ) THEN
409                     !The obs footprint and the model grid box don't overlap so set weight to zero
410                     pweig(ji,jj,jk) = 0.0_wp
411                  ELSE IF ( ( z_cwidth  == zlamscl_m ) .AND. &
412                     &      ( z_cheight == zphiscl_m ) ) THEN
413                     !The grid box is totally contained within the obs footprint so set weight to one
414                     pweig(ji,jj,jk) = 1.0_wp
415                  ELSE IF ( ( z_cwidth  == z_awidth  ) .AND. &
416                     &      ( z_cheight == z_aheight ) ) THEN
417                     !The obs footprint is totally contained within the grid box so set weight as ratio of the two
418                     pweig(ji,jj,jk) = zarea_fp / zareabox
419                  ELSE
420                     !The obs footprint and the grid box overlap so calculate the area of the intersection of the two
421                     zleft   = max(zxvert(1), -zlamscl_m/2.0_wp)
422                     zright  = min(zxvert(2),  zlamscl_m/2.0_wp)
423                     zbottom = max(zyvert(4), -zphiscl_m/2.0_wp)
424                     ztop    = min(zyvert(1),  zphiscl_m/2.0_wp)
425                     
426                     IF (  ( zleft < zright ) .AND. ( zbottom < ztop ) ) THEN
427                        zarea_intersect = ( zright - zleft ) * ( ztop - zbottom )
428                        pweig(ji,jj,jk) = zarea_intersect / zareabox
429                     ENDIF
430                  ENDIF
431
432               END IF !pmask
433            END DO
434         END DO
435      END DO
436
437   END SUBROUTINE obs_avg_h2d_rec
438
439   SUBROUTINE obs_avg_h2d( kpk, kpk2, kmaxifp, kmaxjfp, pweig, pmod, pobsk )
440
441      !!-----------------------------------------------------------------------
442      !!
443      !!                     ***  ROUTINE obs_int_h2d  ***
444      !!
445      !! ** Purpose : Horizontal averaging to the observation footprint.
446      !!
447      !! ** Method  : Average the model points based on the weights already calculated.
448      !!
449      !! ** Action  :
450      !!                   
451      !! References :
452      !!
453      !! History :
454      !!        ! 13/10. M. Martin.
455      !!-----------------------------------------------------------------------
456      !! * Modules used
457      !! * Arguments
458      INTEGER, INTENT(IN) :: &
459         & kpk,   &             ! Parameter values for automatic arrays
460         & kpk2
461      INTEGER, INTENT(IN) :: &
462         & kmaxifp,   &         ! Max size of model points in i-direction for obs footprint
463         & kmaxjfp              ! Max size of model points in j-direction for obs footprint
464      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
465         & pweig                ! Interpolation weights
466      REAL(KIND=wp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
467         & pmod                 ! Model variable to interpolate
468      REAL(KIND=wp), DIMENSION(kpk2), INTENT(OUT) ::  &
469         & pobsk                ! Model profile interpolated to obs (i,j) pt
470
471      INTEGER :: &
472         & jk
473      INTEGER :: &
474         & ikmax
475      REAL(KIND=wp) :: &
476         & zsum
477
478      !------------------------------------------------------------------------
479      ! Initialize number of levels
480      !------------------------------------------------------------------------
481      IF ( kpk2 == 1 ) THEN
482         ikmax = 1
483      ELSEIF ( kpk2 == kpk) THEN
484         ikmax = kpk-1
485      ENDIF
486
487      !------------------------------------------------------------------------
488      ! Average model values to the observation footprint
489      !------------------------------------------------------------------------
490      pobsk = obfillflt
491
492      DO jk = 1, ikmax
493
494         zsum = SUM( pweig(:,:,jk) )
495
496         IF ( zsum /= 0.0_wp ) THEN
497            pobsk(jk) = SUM ( pweig(:,:,jk) * pmod(:,:,jk), Mask=pweig(:,:,jk) > 0.0_wp )
498            pobsk(jk) = pobsk(jk) / zsum
499         END IF
500
501      END DO
502
503   END SUBROUTINE obs_avg_h2d
504
505   SUBROUTINE obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, pmask, kmaxifp, kmaxjfp )
506      !!-----------------------------------------------------------------------
507      !!
508      !!                     ***  ROUTINE obs_max_fpsize  ***
509      !!
510      !! ** Purpose : Calculate maximum number of grid points which may
511      !!              need to be used in the averaging in the global domain.
512      !!
513      !!
514      !! ** Method  : Work out the minimum grid size and work out
515      !!              how many of the smallest grid points would be needed
516      !!              to cover the scale of the observation footprint.
517      !!              This needs to be done using the max/min of the global domain
518      !!              as the obs can be distributed from other parts of the grid.
519      !!
520      !! ** Action  :
521      !!                   
522      !! References :
523      !!
524      !! History :
525      !!        ! 14/01. M. Martin.
526      !!-----------------------------------------------------------------------
527      !! * Modules used
528      !! * Arguments
529      INTEGER , INTENT(IN) :: &
530         & k2dint                   !Type of interpolation/averaging used
531      REAL(KIND=wp), INTENT(IN) :: &
532         & plamscl,   &             !Total width/radius in metres of the observation footprint
533         & pphiscl                  !
534      LOGICAL, INTENT(IN) :: &
535         & lindegrees               !T=> plamscl and pphiscl are specified in degrees
536      REAL(KIND=wp), DIMENSION(jpi,jpj), INTENT(IN) :: &
537         & pmask                    !Land/sea mask
538                                    !F=> plamscl and pphiscl are specified in metres
539      INTEGER, INTENT(OUT)  :: &
540         & kmaxifp,   &             !Max number of grid points in i,j directions to use in averaging
541         & kmaxjfp                  !these have to be even so that the footprint is centred
542
543      !! * Local variables
544      REAL(KIND=wp) :: &
545         & ze1min,     &            !Minimum global grid-size in i,j directions
546         & ze2min
547      REAL(KIND=wp) :: &
548         & zphiscl_m, &
549         & zlamscl_m
550      !------------------------------------------------------------------------
551
552      IF ( k2dint <= 4 ) THEN
553         !If interpolation is being used then only need to use a 2x2 footprint
554         kmaxifp = 2
555         kmaxjfp = 2
556
557      ELSE
558
559         IF ( lindegrees ) THEN
560            !If the scales are specified in degrees, work out the max
561            !distance (metres) in x/y directions
562            CALL obs_deg2dist( jpi, jpj, glamt, gphit, &
563               &               plamscl, pphiscl, zlamscl_m, zphiscl_m )
564         ELSE
565            zlamscl_m = plamscl
566            zphiscl_m = pphiscl
567         ENDIF
568
569         ze1min = MINVAL( e1t(:,:), mask = pmask(:,:) == 1._wp )
570         ze2min = MINVAL( e2t(:,:), mask = pmask(:,:) == 1._wp )
571         
572         IF(lk_mpp) THEN
573            CALL mpp_min( ze1min )
574            CALL mpp_min( ze2min )
575         ENDIF
576
577         kmaxifp = ceiling(zlamscl_m/ze1min) + 1
578         kmaxjfp = ceiling(zphiscl_m/ze2min) + 1
579         
580         !Ensure that these numbers are even
581         kmaxifp = kmaxifp + MOD(kmaxifp,2)
582         kmaxjfp = kmaxjfp + MOD(kmaxjfp,2)
583         
584
585      ENDIF
586
587   END SUBROUTINE obs_max_fpsize
588
589   SUBROUTINE obs_deg2dist( ki, kj, pglam, pgphi, plamscl_deg, pphiscl_deg, &
590      &                     plamscl_max, pphiscl_max )
591      !!-----------------------------------------------------------------------
592      !!
593      !!                     ***  ROUTINE obs_deg2dist  ***
594      !!
595      !! ** Purpose : Calculate the maximum distance in m of the length scale
596      !!              in degrees.
597      !!
598      !! ** Method  : At each lon/lat point, work out the distances in the
599      !!              zonal and meridional directions.
600      !!
601      !! ** Action  :
602      !!                   
603      !! References :
604      !!
605      !! History :
606      !!        ! 14/01. M. Martin.
607      !!-----------------------------------------------------------------------
608      !! * Modules used
609      !! * Arguments
610      INTEGER , INTENT(IN) :: &
611         & ki, kj                   !x/y dimensions of input lat/lon variables
612      REAL(KIND=wp), INTENT(IN), DIMENSION(ki,kj) :: &
613         & pglam, pgphi             !Longitude and latitudes of grid points
614      REAL(KIND=wp), INTENT(IN) :: &
615         & plamscl_deg,   &         !Size in degrees of the observation footprint
616         & pphiscl_deg              !
617      REAL(KIND=wp), INTENT(OUT) :: &
618         & plamscl_max, &           !Maximum size of obs footprint in metres
619         & pphiscl_max
620     
621      !! * Local declarations
622      INTEGER :: &
623         & ji, jj                   !Counters
624      REAL(KIND=wp) :: &
625         & zlon1, zlon2, &          !Lon values surrounding centre of grid point
626         & zlat1, zlat2, &          !Lat values surrounding centre of grid point
627         & zdlat, zdlon             !Distance in radians in lat/lon directions
628      REAL(KIND=wp) :: &
629         & za1, za2, za, zc, zd
630     
631      plamscl_max = -1.0_wp
632      pphiscl_max = -1.0_wp
633
634      DO ji = 1, ki
635         DO jj = 1, kj
636
637            !Calculate distance in metres in zonal(x) direction
638
639            zlon1 = rad * ( pglam(ji,jj) + ( 0.5_wp * plamscl_deg ) )
640            zlon2 = rad * ( pglam(ji,jj) - ( 0.5_wp * plamscl_deg ) )
641            zlat1 = rad * pgphi(ji,jj)
642            zlat2 = rad * pgphi(ji,jj)
643            zdlon = zlon2 - zlon1
644            zdlat = zlat2 - zlat1
645
646            za1 = sin( zdlat/2.0_wp )
647            za2 = sin( zdlon/2.0_wp )
648            za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
649            zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
650            zd = ra * zc
651
652            IF ( zd > plamscl_max ) plamscl_max = zd
653
654            !Calculate distance in metres in meridional(y) direction
655
656            zlon1 = rad * pglam(ji,jj)
657            zlon2 = rad * pglam(ji,jj)
658            zlat1 = rad * ( pgphi(ji,jj) + ( 0.5_wp * pphiscl_deg ) )
659            zlat2 = rad * ( pgphi(ji,jj) - ( 0.5_wp * pphiscl_deg ) )
660            zdlon = zlon2 - zlon1
661            zdlat = zlat2 - zlat1
662
663            za1 = sin( zdlat/2.0_wp )
664            za2 = sin( zdlon/2.0_wp )
665            za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
666            zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
667            zd = ra * zc
668
669            IF ( zd > pphiscl_max ) pphiscl_max = zd
670
671         END DO
672      END DO
673         
674   END SUBROUTINE obs_deg2dist
675
676   SUBROUTINE obs_dist2corners(pglam_bl, pglam_br, pglam_tl, pglam_tr, &
677      &                        pgphi_bl, pgphi_br, pgphi_tl, pgphi_tr, &
678      &                        plam, pphi, pxvert, pyvert, pdist)
679      !!-----------------------------------------------------------------------
680      !!
681      !!                     ***  ROUTINE obs_dist2corners  ***
682      !!
683      !! ** Purpose : Calculate distance from centre of obs footprint to the corners of a grid box
684      !!
685      !! ** Method  : Use great circle distance formulae.
686      !!              Order of corners is 1=topleft, 2=topright, 3=bottomright, 4=bottomleft
687      !!
688      !! ** Action  :
689      !!                   
690      !! References :
691      !!
692      !! History :
693      !!        ! 14/01. M. Martin.
694      !!-----------------------------------------------------------------------
695      !! * Modules used
696      !! * Arguments
697         REAL(KIND=wp), INTENT(IN) :: &
698            &  pglam_bl, pglam_br, &       !lon at corners of grid box
699            &  pglam_tl, pglam_tr
700         REAL(KIND=wp), INTENT(IN) :: &
701            &  pgphi_bl, pgphi_br, &       !lat at corners of grid box
702            &  pgphi_tl, pgphi_tr
703         REAL(KIND=wp), INTENT(IN) :: &
704            &  pphi, plam                  !lat/lon of centre of obs footprint
705         REAL(KIND=wp), DIMENSION(4), INTENT(OUT) :: &
706            &  pxvert, pyvert              !x/y location (in metres relative to centre of obs footprint) of corners
707         REAL(KIND=wp), DIMENSION(4), INTENT(OUT) :: &
708            &  pdist                       !distance (in metres) of each corner relative to centre of obs footprint
709
710      !! * Local variables
711         INTEGER :: &
712            &  jvert                       !Counter for corners
713         REAL(KIND=wp) :: &
714            &  zphi, zlam                  !Local values for lon/lat of corners
715         REAL(KIND=wp) :: &
716            &  za1, za2,  &                !For great circle distance calculations
717            &  zb1, zb2,  &
718            &  zc1, zc2
719         REAL(KIND=wp) :: &
720            &  zdist_centre_lat, &         !Distances in lat/lon directions (in metres)
721            &  zdist_centre_lon
722
723      !!-----------------------------------------------------------------------
724
725         ! Work out latitudinal and longitudinal distance from centre of
726         ! obs fp to corners of grid box
727         DO jvert = 1, 4
728            SELECT CASE(jvert)
729            CASE(1)
730               zphi = pgphi_tl
731               zlam = pglam_tl
732            CASE(2)
733               zphi = pgphi_tr
734               zlam = pglam_tr
735            CASE(3)
736               zphi = pgphi_br
737               zlam = pglam_br
738            CASE(4)
739               zphi = pgphi_bl
740               zlam = pglam_bl
741            END SELECT
742           
743            IF (zlam == plam ) THEN
744               pxvert(jvert) = 0.0_wp
745            ELSE
746               za1 = SIN( zphi * rad )
747               za2 = SIN( zphi * rad )
748               zb1 = COS( zphi * rad ) * COS( zlam * rad )
749               zb2 = COS( zphi * rad ) * COS( plam  * rad )
750               zc1 = COS( zphi * rad ) * SIN( zlam * rad )
751               zc2 = COS( zphi * rad ) * SIN( plam  * rad )
752               pxvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
753               pxvert(jvert) =  ra * pxvert(jvert)
754               IF ( zlam < plam ) pxvert(jvert) = - pxvert(jvert)
755            ENDIF
756           
757            IF ( zphi == pphi ) THEN
758               pyvert(jvert) = 0.0_wp
759            ELSE
760               za1 = SIN( zphi * rad )
761               za2 = SIN( pphi * rad )
762               zb1 = COS( zphi * rad ) * COS( zlam * rad )
763               zb2 = COS( pphi * rad ) * COS( zlam * rad )
764               zc1 = COS( zphi * rad ) * SIN( zlam * rad )
765               zc2 = COS( pphi * rad ) * SIN( zlam * rad )
766               pyvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
767               pyvert(jvert) =  ra * pyvert(jvert)
768               IF ( zphi < pphi ) pyvert(jvert) = - pyvert(jvert)
769            ENDIF
770
771            !Calculate the distance of each vertex relative to centre of obs footprint
772            pdist(jvert) = SQRT( ( pxvert(jvert) * pxvert(jvert) ) + &
773               &                 ( pyvert(jvert) * pyvert(jvert) ) )
774
775         END DO
776
777      END SUBROUTINE obs_dist2corners
Note: See TracBrowser for help on using the repository browser.