Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

obs_averg_h2d.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago
Mixed precision version, tested up to 30 years on ORCA2.
Property svn:keywords set to `Id`
File size: 36.1 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	!! * Substitutions
45	# include "single_precision_substitute.h90"
46
47	!!----------------------------------------------------------------------
48	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
49	!! $Id$
50	!! Software governed by the CeCILL license (see ./LICENSE)
51	!!----------------------------------------------------------------------
52
53	CONTAINS
54	SUBROUTINE obs_avg_h2d_init( kpk, kpk2, kmaxifp, kmaxjfp, k2dint, plam, pphi, &
55	& pglam, pgphi, pglamf, pgphif, pmask, plamscl, pphiscl, lindegrees, &
56	& pweig, iminpoints )
57	!!-----------------------------------------------------------------------
58	!!
59	!! * ROUTINE obs_avg_h2d_init *
60	!!
61	!! ** Purpose : Computes weights for horizontal averaging to the
62	!! observation footprint.
63	!!
64	!! ** Method : Horizontal averaging to the observation footprint using
65	!! model values at a defined area.
66	!!
67	!! Averaging schemes :
68	!!
69	!! Two horizontal averaging schemes are available:
70	!! - weighted radial footprint (k2dint = 5)
71	!! - weighted rectangular footprint (k2dint = 6)
72	!!
73	!! History :
74	!! ! 13-10 (M. Martin)
75	!!-----------------------------------------------------------------------
76	!! * Modules used
77	!! * Arguments
78	INTEGER, INTENT(IN) :: &
79	& kpk, & ! Parameter values for automatic arrays
80	& kpk2, &
81	& kmaxifp, & ! Max size of model points in i-direction for obs footprint
82	& kmaxjfp, & ! Max size of model points in j-direction for obs footprint
83	& k2dint ! Averaging scheme options - see header
84	REAL(KIND=dp), INTENT(INOUT) :: &
85	& plam, & ! Geographical (lat,lon) coordinates of
86	& pphi ! observation
87	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
88	& pglam, & ! Model variable lon
89	& pgphi ! Model variable lat
90	REAL(KIND=dp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
91	& pglamf, & ! Model variable lon at corners of grid-boxes
92	& pgphif ! Model variable lat at corners of grid-boxes
93	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
94	& pmask ! Model variable mask
95	REAL(KIND=dp), INTENT(IN) :: &
96	& plamscl, & ! Diameter (lat,lon) of obs footprint in metres
97	& pphiscl ! This is the full width (rather than half-width)
98	LOGICAL, INTENT(IN) :: &
99	& lindegrees ! T=> obs footprint specified in degrees, F=> in metres
100	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) :: &
101	& pweig ! Weights for averaging
102	INTEGER, INTENT(IN), OPTIONAL :: &
103	& iminpoints ! Reject point which is not surrounded
104	! by at least iminpoints sea points
105
106	!! * Local declarations
107	INTEGER :: &
108	& jk
109	INTEGER :: &
110	& ikmax
111
112
113	!------------------------------------------------------------------------
114	!
115	!------------------------------------------------------------------------
116
117	!------------------------------------------------------------------------
118	! Initialize number of levels
119	!------------------------------------------------------------------------
120	IF ( kpk2 == 1 ) THEN
121	ikmax = 1
122	ELSEIF ( kpk2 == kpk) THEN
123	ikmax = kpk-1
124	ENDIF
125
126
127	SELECT CASE (k2dint)
128	CASE(5)
129	CALL obs_avg_h2d_rad( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
130	& plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
131	CASE(6)
132	CALL obs_avg_h2d_rec( kpk2, ikmax, kmaxifp, kmaxjfp, plam, pphi, &
133	& plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
134	END SELECT
135
136
137	END SUBROUTINE obs_avg_h2d_init
138
139
140	SUBROUTINE obs_avg_h2d_rad( kpk2, kmax, kmaxifp, kmaxjfp, plam, pphi, &
141	& plamscl, pphiscl, lindegrees, pmask, pglam, pgphi, pglamf, pgphif, pweig )
142	!!-----------------------------------------------------------------------
143	!!
144	!! * ROUTINE obs_avg_h2d_rad *
145	!!
146	!! ** Purpose : Computes weights for horizontal averaging to the
147	!! observation using a radial footprint.
148	!!
149	!! ** Method : Calculate whether each grid box is completely or
150	!! partially within the observation footprint.
151	!! If it is partially in the footprint then calculate
152	!! the ratio of the area inside the footprint to the total
153	!! area of the grid box.
154	!!
155	!! History :
156	!! ! 14-01 (M. Martin)
157	!!-----------------------------------------------------------------------
158	!! * Modules used
159	USE phycst, ONLY : & ! Physical constants
160	& ra, &
161	& rpi
162
163	!! * Arguments
164	INTEGER, INTENT(IN) :: &
165	& kpk2, & ! Parameter values for automatic arrays
166	& kmax
167
168	INTEGER, INTENT(IN) :: &
169	& kmaxifp, & ! Max size of model points in i-direction for obs footprint
170	& kmaxjfp ! Max size of model points in j-direction for obs footprint
171
172	REAL(KIND=dp), INTENT(IN) :: plam, pphi
173
174
175
176	! observation
177	REAL(KIND=dp), 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=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
183	& pmask ! Model variable mask
184	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
185	& pglam, & ! Model variable lon
186	& pgphi ! Model variable lat
187	REAL(KIND=dp), 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=dp), 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=dp), 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=dp), 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=dp) :: &
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=dp) :: &
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=dp), INTENT(IN) :: plam, pphi ! Geographical (lat,lon) coordinates of
357
358
359
360	! observation
361	REAL(KIND=dp), INTENT(IN) :: plamscl, pphiscl ! Width in x/y directions of obs footprint in metres
362
363
364
365	! This is the full width (rather than half-width)
366	LOGICAL, INTENT(IN) :: &
367	& lindegrees !T=> scales specified in degrees, F=> in metres
368	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
369	& pmask ! Model variable mask
370	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp), INTENT(IN) :: &
371	& pglam, & ! Model variable lat at centre of grid boxes
372	& pgphi ! Model variable lon at centre of grid boxes
373	REAL(KIND=dp), DIMENSION(kmaxifp+1,kmaxjfp+1), INTENT(IN) :: &
374	& pglamf, & ! Model variable lat at corners of grid boxes
375	& pgphif ! Model variable lon at corners of grid boxes
376	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(OUT) :: &
377	& pweig ! Weights for interpolation
378
379	!! Local declarations
380	INTEGER :: ji, jj, jk
381	INTEGER :: jvert
382	INTEGER, DIMENSION(4) :: &
383	& jnumvert
384	REAL(KIND=dp), DIMENSION(4) :: &
385	& zxvert, zyvert ! The lon/lat of the vertices(corners) of the grid box in m relative to centre of obs fp
386	REAL(KIND=dp), DIMENSION(4) :: &
387	& zdist ! Distance of each vertex to the centre of the obs footprint
388	REAL(KIND=dp), DIMENSION(4) :: &
389	& zxgrid, zygrid, & ! Distance of each vertex of grid box to the centre of the grid box in x/y directions
390	& zdgrid
391	REAL(KIND=dp) :: &
392	& zareabox, & ! Total area of grid box
393	& zarea_fp, & ! Total area of obs footprint
394	& zarea_intersect ! Area of the intersection between the grid box and the obs footprint
395	REAL(KIND=dp) :: &
396	& zlamscl_m, &
397	& zphiscl_m ! Total width (lat,lon) of obs footprint in metres
398	REAL(KIND=dp) :: &
399	& z_awidth, z_aheight, & ! Width and height of model grid box
400	& z_cwidth, z_cheight ! Width and height of union of model grid box and obs footprint
401	REAL(KIND=dp) :: &
402	& zleft, zright, & ! Distance (metres) of corners area of intersection
403	& ztop, zbottom ! between grid box and obs footprint
404
405	!-----------------------------------------------------------------------
406
407	!Initialise weights to zero
408	pweig(:,:,:) = 0.0_wp
409
410	!Loop over the grid boxes which have been identified as potentially being within the
411	!observation footprint
412	DO jk = 1, kmax
413	DO ji = 1, kmaxifp
414	DO jj = 1, kmaxjfp
415
416	IF ( pmask(ji,jj,jk) == 1.0_wp ) THEN
417
418
419	IF ( lindegrees ) THEN
420	!If the scales are specified in degrees, work out the
421	!scales (metres) in x/y directions
422	CALL obs_deg2dist( 1, 1, pglam(ji,jj), pgphi(ji,jj), &
423	& plamscl, pphiscl, zlamscl_m, zphiscl_m )
424	ELSE
425	zlamscl_m = plamscl
426	zphiscl_m = pphiscl
427	ENDIF
428
429	! Work out the area of the grid box using distance of corners relative to centre of grid box
430	CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
431	& pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
432	& pglam(ji,jj), pgphi(ji,jj), zxgrid, zygrid, zdgrid)
433
434	!Calculate width and height of model grid box
435	z_awidth = ABS( zxgrid(1) - zxgrid(2) )
436	z_aheight = ABS( zygrid(1) - zygrid(4) )
437	zareabox = z_awidth * z_aheight
438
439	! Work out area of the observation footprint
440	zarea_fp = zlamscl_m * zphiscl_m
441
442	! For each corner of the grid box, calculate its location and distance relative
443	! to the centre of the observation footprint
444	CALL obs_dist2corners(pglamf(ji,jj), pglamf(ji+1,jj), pglamf(ji,jj+1), pglamf(ji+1,jj+1), &
445	& pgphif(ji,jj), pgphif(ji+1,jj), pgphif(ji,jj+1), pgphif(ji+1,jj+1), &
446	& plam, pphi, zxvert, zyvert, zdist)
447
448	!Work out maximum width and height of rectangle covered by corners of obs fp and grid box
449	z_cwidth = MAX( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp ) - &
450	& MIN( zxvert(1), zxvert(2), -zlamscl_m/2.0_wp, zlamscl_m/2.0_wp )
451
452	z_cheight = MAX( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp ) - &
453	& MIN( zyvert(1), zyvert(4), zphiscl_m/2.0_wp, -zphiscl_m/2.0_wp )
454
455	IF ( ( z_cwidth >= z_awidth + zlamscl_m ) .OR. &
456	& ( z_cheight >= z_aheight + zphiscl_m ) ) THEN
457	!The obs footprint and the model grid box don't overlap so set weight to zero
458	pweig(ji,jj,jk) = 0.0_wp
459	ELSE IF ( ( z_cwidth == zlamscl_m ) .AND. &
460	& ( z_cheight == zphiscl_m ) ) THEN
461	!The grid box is totally contained within the obs footprint so set weight to one
462	pweig(ji,jj,jk) = 1.0_wp
463	ELSE IF ( ( z_cwidth == z_awidth ) .AND. &
464	& ( z_cheight == z_aheight ) ) THEN
465	!The obs footprint is totally contained within the grid box so set weight as ratio of the two
466	pweig(ji,jj,jk) = zarea_fp / zareabox
467	ELSE
468	!The obs footprint and the grid box overlap so calculate the area of the intersection of the two
469	zleft = max(zxvert(1), -zlamscl_m/2.0_wp)
470	zright = min(zxvert(2), zlamscl_m/2.0_wp)
471	zbottom = max(zyvert(4), -zphiscl_m/2.0_wp)
472	ztop = min(zyvert(1), zphiscl_m/2.0_wp)
473
474	IF ( ( zleft < zright ) .AND. ( zbottom < ztop ) ) THEN
475	zarea_intersect = ( zright - zleft ) * ( ztop - zbottom )
476	pweig(ji,jj,jk) = zarea_intersect / zareabox
477	ENDIF
478	ENDIF
479
480	END IF !pmask
481	END DO
482	END DO
483	END DO
484
485	END SUBROUTINE obs_avg_h2d_rec
486
487	SUBROUTINE obs_avg_h2d( kpk, kpk2, kmaxifp, kmaxjfp, pweig, pmod, pobsk )
488
489	!!-----------------------------------------------------------------------
490	!!
491	!! * ROUTINE obs_int_h2d *
492	!!
493	!! ** Purpose : Horizontal averaging to the observation footprint.
494	!!
495	!! ** Method : Average the model points based on the weights already calculated.
496	!!
497	!! ** Action :
498	!!
499	!! References :
500	!!
501	!! History :
502	!! ! 13/10. M. Martin.
503	!!-----------------------------------------------------------------------
504	!! * Modules used
505	!! * Arguments
506	INTEGER, INTENT(IN) :: &
507	& kpk, & ! Parameter values for automatic arrays
508	& kpk2
509	INTEGER, INTENT(IN) :: &
510	& kmaxifp, & ! Max size of model points in i-direction for obs footprint
511	& kmaxjfp ! Max size of model points in j-direction for obs footprint
512	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
513	& pweig ! Interpolation weights
514	REAL(KIND=dp), DIMENSION(kmaxifp,kmaxjfp,kpk2), INTENT(IN) :: &
515	& pmod ! Model variable to interpolate
516	REAL(KIND=dp), DIMENSION(kpk2), INTENT(OUT) :: &
517	& pobsk ! Model profile interpolated to obs (i,j) pt
518
519	INTEGER :: &
520	& jk
521	INTEGER :: &
522	& ikmax
523	REAL(KIND=dp) :: zsum
524
525
526
527	!------------------------------------------------------------------------
528	! Initialize number of levels
529	!------------------------------------------------------------------------
530	IF ( kpk2 == 1 ) THEN
531	ikmax = 1
532	ELSEIF ( kpk2 == kpk) THEN
533	ikmax = kpk-1
534	ENDIF
535
536	!------------------------------------------------------------------------
537	! Average model values to the observation footprint
538	!------------------------------------------------------------------------
539	pobsk = obfillflt
540
541	DO jk = 1, ikmax
542
543	zsum = SUM( pweig(:,:,jk) )
544
545	IF ( zsum /= 0.0_wp ) THEN
546	pobsk(jk) = SUM ( pweig(:,:,jk) * pmod(:,:,jk), Mask=pweig(:,:,jk) > 0.0_wp )
547	pobsk(jk) = pobsk(jk) / zsum
548	END IF
549
550	END DO
551
552	END SUBROUTINE obs_avg_h2d
553
554	SUBROUTINE obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, pmask, kmaxifp, kmaxjfp )
555	!!-----------------------------------------------------------------------
556	!!
557	!! * ROUTINE obs_max_fpsize *
558	!!
559	!! ** Purpose : Calculate maximum number of grid points which may
560	!! need to be used in the averaging in the global domain.
561	!!
562	!!
563	!! ** Method : Work out the minimum grid size and work out
564	!! how many of the smallest grid points would be needed
565	!! to cover the scale of the observation footprint.
566	!! This needs to be done using the max/min of the global domain
567	!! as the obs can be distributed from other parts of the grid.
568	!!
569	!! ** Action :
570	!!
571	!! References :
572	!!
573	!! History :
574	!! ! 14/01. M. Martin.
575	!!-----------------------------------------------------------------------
576	!! * Modules used
577	!! * Arguments
578	INTEGER , INTENT(IN) :: &
579	& k2dint !Type of interpolation/averaging used
580	REAL(KIND=dp), INTENT(IN) :: &
581	& plamscl, & !Total width/radius in metres of the observation footprint
582	& pphiscl !
583	LOGICAL, INTENT(IN) :: &
584	& lindegrees !T=> plamscl and pphiscl are specified in degrees
585	REAL(KIND=dp), DIMENSION(jpi,jpj), INTENT(IN) :: &
586	& pmask !Land/sea mask
587	!F=> plamscl and pphiscl are specified in metres
588	INTEGER, INTENT(OUT) :: &
589	& kmaxifp, & !Max number of grid points in i,j directions to use in averaging
590	& kmaxjfp !these have to be even so that the footprint is centred
591
592	!! * Local variables
593	REAL(KIND=dp) :: &
594	& ze1min, & !Minimum global grid-size in i,j directions
595	& ze2min
596	REAL(KIND=dp) :: &
597	& zphiscl_m, &
598	& zlamscl_m
599	!------------------------------------------------------------------------
600
601	IF ( k2dint <= 4 ) THEN
602	!If interpolation is being used then only need to use a 2x2 footprint
603	kmaxifp = 2
604	kmaxjfp = 2
605
606	ELSE
607
608	IF ( lindegrees ) THEN
609	!If the scales are specified in degrees, work out the max
610	!distance (metres) in x/y directions
611	CALL obs_deg2dist( jpi, jpj, glamt, gphit, &
612	& plamscl, pphiscl, zlamscl_m, zphiscl_m )
613	ELSE
614	zlamscl_m = plamscl
615	zphiscl_m = pphiscl
616	ENDIF
617
618	ze1min = MINVAL( e1t(:,:), mask = pmask(:,:) == 1._wp )
619	ze2min = MINVAL( e2t(:,:), mask = pmask(:,:) == 1._wp )
620
621	IF(lk_mpp) THEN
622	CALL mpp_min( 'obs_averg_h2d', ze1min )
623	CALL mpp_min( 'obs_averg_h2d', ze2min )
624	ENDIF
625
626	kmaxifp = ceiling(zlamscl_m/ze1min) + 1
627	kmaxjfp = ceiling(zphiscl_m/ze2min) + 1
628
629	!Ensure that these numbers are even
630	kmaxifp = kmaxifp + MOD(kmaxifp,2)
631	kmaxjfp = kmaxjfp + MOD(kmaxjfp,2)
632
633
634	ENDIF
635
636	END SUBROUTINE obs_max_fpsize
637
638	SUBROUTINE obs_deg2dist( ki, kj, pglam, pgphi, plamscl_deg, pphiscl_deg, &
639	& plamscl_max, pphiscl_max )
640	!!-----------------------------------------------------------------------
641	!!
642	!! * ROUTINE obs_deg2dist *
643	!!
644	!! ** Purpose : Calculate the maximum distance in m of the length scale
645	!! in degrees.
646	!!
647	!! ** Method : At each lon/lat point, work out the distances in the
648	!! zonal and meridional directions.
649	!!
650	!! ** Action :
651	!!
652	!! References :
653	!!
654	!! History :
655	!! ! 14/01. M. Martin.
656	!!-----------------------------------------------------------------------
657	!! * Modules used
658	!! * Arguments
659	INTEGER , INTENT(IN) :: &
660	& ki, kj !x/y dimensions of input lat/lon variables
661	REAL(KIND=dp), INTENT(IN), DIMENSION(ki,kj) :: &
662	& pglam, pgphi !Longitude and latitudes of grid points
663	REAL(KIND=dp), INTENT(IN) :: &
664	& plamscl_deg, & !Size in degrees of the observation footprint
665	& pphiscl_deg !
666	REAL(KIND=dp), INTENT(OUT) :: &
667	& plamscl_max, & !Maximum size of obs footprint in metres
668	& pphiscl_max
669
670	!! * Local declarations
671	INTEGER :: &
672	& ji, jj !Counters
673	REAL(KIND=dp) :: &
674	& zlon1, zlon2, & !Lon values surrounding centre of grid point
675	& zlat1, zlat2, & !Lat values surrounding centre of grid point
676	& zdlat, zdlon !Distance in radians in lat/lon directions
677	REAL(KIND=dp) :: &
678	& za1, za2, za, zc, zd
679
680	plamscl_max = -1.0_wp
681	pphiscl_max = -1.0_wp
682
683	DO ji = 1, ki
684	DO jj = 1, kj
685
686	!Calculate distance in metres in zonal(x) direction
687
688	zlon1 = rad * ( pglam(ji,jj) + ( 0.5_wp * plamscl_deg ) )
689	zlon2 = rad * ( pglam(ji,jj) - ( 0.5_wp * plamscl_deg ) )
690	zlat1 = rad * pgphi(ji,jj)
691	zlat2 = rad * pgphi(ji,jj)
692	zdlon = zlon2 - zlon1
693	zdlat = zlat2 - zlat1
694
695	za1 = sin( zdlat/2.0_wp )
696	za2 = sin( zdlon/2.0_wp )
697	za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
698	zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
699	zd = ra * zc
700
701	IF ( zd > plamscl_max ) plamscl_max = zd
702
703	!Calculate distance in metres in meridional(y) direction
704
705	zlon1 = rad * pglam(ji,jj)
706	zlon2 = rad * pglam(ji,jj)
707	zlat1 = rad * ( pgphi(ji,jj) + ( 0.5_wp * pphiscl_deg ) )
708	zlat2 = rad * ( pgphi(ji,jj) - ( 0.5_wp * pphiscl_deg ) )
709	zdlon = zlon2 - zlon1
710	zdlat = zlat2 - zlat1
711
712	za1 = sin( zdlat/2.0_wp )
713	za2 = sin( zdlon/2.0_wp )
714	za = ( za1 * za1 ) + ( COS( zlat1 ) * COS( zlat2 ) * ( za2 * za2 ) )
715	zc = 2.0_wp * atan2( SQRT( za ), SQRT( 1.0_wp-za ) )
716	zd = ra * zc
717
718	IF ( zd > pphiscl_max ) pphiscl_max = zd
719
720	END DO
721	END DO
722
723	END SUBROUTINE obs_deg2dist
724
725	SUBROUTINE obs_dist2corners(pglam_bl, pglam_br, pglam_tl, pglam_tr, &
726	& pgphi_bl, pgphi_br, pgphi_tl, pgphi_tr, &
727	& plam, pphi, pxvert, pyvert, pdist)
728	!!-----------------------------------------------------------------------
729	!!
730	!! * ROUTINE obs_dist2corners *
731	!!
732	!! ** Purpose : Calculate distance from centre of obs footprint to the corners of a grid box
733	!!
734	!! ** Method : Use great circle distance formulae.
735	!! Order of corners is 1=topleft, 2=topright, 3=bottomright, 4=bottomleft
736	!!
737	!! ** Action :
738	!!
739	!! References :
740	!!
741	!! History :
742	!! ! 14/01. M. Martin.
743	!!-----------------------------------------------------------------------
744	!! * Modules used
745	!! * Arguments
746	REAL(KIND=dp), INTENT(IN) :: &
747	& pglam_bl, pglam_br, & !lon at corners of grid box
748	& pglam_tl, pglam_tr
749	REAL(KIND=dp), INTENT(IN) :: &
750	& pgphi_bl, pgphi_br, & !lat at corners of grid box
751	& pgphi_tl, pgphi_tr
752	REAL(KIND=dp), INTENT(IN) :: &
753	& pphi, plam !lat/lon of centre of obs footprint
754	REAL(KIND=dp), DIMENSION(4), INTENT(OUT) :: &
755	& pxvert, pyvert !x/y location (in metres relative to centre of obs footprint) of corners
756	REAL(KIND=dp), DIMENSION(4), INTENT(OUT) :: &
757	& pdist !distance (in metres) of each corner relative to centre of obs footprint
758
759	!! * Local variables
760	INTEGER :: &
761	& jvert !Counter for corners
762	REAL(KIND=dp) :: &
763	& zphi, zlam !Local values for lon/lat of corners
764	REAL(KIND=dp) :: &
765	& za1, za2, & !For great circle distance calculations
766	& zb1, zb2, &
767	& zc1, zc2
768	REAL(KIND=dp) :: &
769	& zdist_centre_lat, & !Distances in lat/lon directions (in metres)
770	& zdist_centre_lon
771
772	!!-----------------------------------------------------------------------
773
774	! Work out latitudinal and longitudinal distance from centre of
775	! obs fp to corners of grid box
776	DO jvert = 1, 4
777	SELECT CASE(jvert)
778	CASE(1)
779	zphi = pgphi_tl
780	zlam = pglam_tl
781	CASE(2)
782	zphi = pgphi_tr
783	zlam = pglam_tr
784	CASE(3)
785	zphi = pgphi_br
786	zlam = pglam_br
787	CASE(4)
788	zphi = pgphi_bl
789	zlam = pglam_bl
790	END SELECT
791
792	IF (zlam == plam ) THEN
793	pxvert(jvert) = 0.0_wp
794	ELSE
795	za1 = SIN( zphi * rad )
796	za2 = SIN( zphi * rad )
797	zb1 = COS( zphi * rad ) * COS( zlam * rad )
798	zb2 = COS( zphi * rad ) * COS( plam * rad )
799	zc1 = COS( zphi * rad ) * SIN( zlam * rad )
800	zc2 = COS( zphi * rad ) * SIN( plam * rad )
801	pxvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
802	pxvert(jvert) = ra * pxvert(jvert)
803	IF ( zlam < plam ) pxvert(jvert) = - pxvert(jvert)
804	ENDIF
805
806	IF ( zphi == pphi ) THEN
807	pyvert(jvert) = 0.0_wp
808	ELSE
809	za1 = SIN( zphi * rad )
810	za2 = SIN( pphi * rad )
811	zb1 = COS( zphi * rad ) * COS( zlam * rad )
812	zb2 = COS( pphi * rad ) * COS( zlam * rad )
813	zc1 = COS( zphi * rad ) * SIN( zlam * rad )
814	zc2 = COS( pphi * rad ) * SIN( zlam * rad )
815	pyvert(jvert) = grt_cir_dis( za1, za2, zb1, zb2, zc1, zc2 )
816	pyvert(jvert) = ra * pyvert(jvert)
817	IF ( zphi < pphi ) pyvert(jvert) = - pyvert(jvert)
818	ENDIF
819
820	!Calculate the distance of each vertex relative to centre of obs footprint
821	pdist(jvert) = SQRT( ( pxvert(jvert) * pxvert(jvert) ) + &
822	& ( pyvert(jvert) * pyvert(jvert) ) )
823
824	END DO
825
826	END SUBROUTINE obs_dist2corners
827
828	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/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_averg_h2d.F90 @ 15540

Download in other formats: