source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/OBS/obs_averg_h2d.F90 @ 10297

Last change on this file since 10297 was 10297, checked in by smasson, 23 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of mppmin/max/sum, see #2133

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