Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

obs_averg_h2d.F90

Last change on this file was 14275, checked in by smasson, 3 years ago
trunk: suppress nproc ( = mpprank = narea-1)
Property svn:keywords set to `Id`
File size: 36.2 KB

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

Note: See TracBrowser for help on using the repository browser.

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

Context Navigation

source: NEMO/trunk/src/OCE/OBS/obs_averg_h2d.F90

Download in other formats: