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/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS/obs_averg_h2d.F90 @ 14219

Last change on this file since 14219 was 14219, checked in by mcastril, 4 years ago

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