New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_averg_h2d.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/OBS – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/OBS/obs_averg_h2d.F90 @ 12182

Last change on this file since 12182 was 12182, checked in by davestorkey, 4 years ago

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