Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obsaverg_h2d.h90 @ 8646

Last change on this file since 8646 was 7992, checked in by jwhile, 7 years ago
This version of the code seems to work correctly after some bug fixes
File size: 34.6 KB

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

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

Download in other formats:

Original Format

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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_medusa/NEMOGCM/NEMO/OPA_SRC/OBS/obsaverg_h2d.h90 @ 8646

Download in other formats: