-
Notifications
You must be signed in to change notification settings - Fork 145
/
location_mod.f90
2430 lines (1961 loc) · 94 KB
/
location_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!> Implements location interfaces for a three dimensional spherical shell
!> with a choice of vertical coordinates.
!> Horizontal coordinates are always latitude and longitude.
!> Vertical coordinate choices include pressure, height, model level,
!> scale height, surface, and non-specific (column-integrated values, or
!> with no logically defined vertical location, e.g. hurricane vortex center)
!> The internal representation of the location is stored as
!> radians from 0 to 2 PI for longitude and -PI/2 to PI/2 for latitude to
!> minimize computational cost for distances. However, the external
!> representation is longitude in degrees from 0 to 360 and latitude
!> from -90 to 90 for consistency with most applications in the field.
!>
!> This version supports multiple cutoff distances in an efficient manner.
!> Smaller cutoff values will do less searching than larger ones. (This was
!> not true in earlier implementations of this code.)
!>
module location_mod
use types_mod, only : r8, MISSING_R8, MISSING_I, PI, RAD2DEG, DEG2RAD, OBSTYPELENGTH, i8
use utilities_mod, only : error_handler, E_ERR, ascii_file_format, &
E_MSG, open_file, close_file, set_output, &
logfileunit, nmlfileunit, find_namelist_in_file, &
check_namelist_read, do_output, do_nml_file, &
do_nml_term, is_longitude_between
use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform
use obs_kind_mod, only : get_num_types_of_obs, get_name_for_type_of_obs, get_index_for_type_of_obs
use mpi_utilities_mod, only : my_task_id, task_count
use ensemble_manager_mod, only : ensemble_type
implicit none
private
public :: location_type, get_location, set_location, &
set_location_missing, is_location_in_region, get_maxdist, &
write_location, read_location, interactive_location, query_location, &
LocationDims, LocationName, LocationLName, LocationStorageOrder, LocationUnits, &
get_close_type, get_close_init, get_close_obs, get_close_state, get_close_destroy, &
operator(==), operator(/=), get_dist, has_vertical_choice, vertical_localization_on, &
set_vertical, is_vertical, get_vertical_localization_coord, get_close, &
set_vertical_localization_coord, convert_vertical_obs, convert_vertical_state, &
VERTISUNDEF, VERTISSURFACE, VERTISLEVEL, VERTISPRESSURE, &
VERTISHEIGHT, VERTISSCALEHEIGHT, print_get_close_type
character(len=*), parameter :: source = 'threed_sphere/location_mod.f90'
integer, parameter :: LocationDims = 3
character(len=*), parameter :: LocationName = "loc3Dsphere"
character(len=*), parameter :: LocationLName = &
"threed sphere locations: lon, lat, vertical"
character(len=*), parameter :: LocationStorageOrder = "Lon Lat Vertical"
character(len=*), parameter :: LocationUnits = "degrees degrees which_vert"
! The possible numeric values for the location_type%which_vert component.
! The numeric values are PRIVATE to this module. The parameter names are PUBLIC.
integer, parameter :: VERTISUNDEF = -2 ! has no specific vertical location (undefined)
integer, parameter :: VERTISSURFACE = -1 ! surface value (value is surface elevation in m)
integer, parameter :: VERTISLEVEL = 1 ! by level
integer, parameter :: VERTISPRESSURE = 2 ! by pressure (in pascals)
integer, parameter :: VERTISHEIGHT = 3 ! by height (in meters)
integer, parameter :: VERTISSCALEHEIGHT = 4 ! by scale height (unitless)
type location_type
private
real(r8) :: lon, lat ! lon, lat are stored in radians
real(r8) :: vloc ! units vary based on value of which_vert
integer :: which_vert ! determines if vert is level, height, pressure, ...
end type location_type
! Derived type to facilitate efficient computation of locations close to a given observation.
type get_close_type_by_type
private
integer :: num
real(r8) :: maxdist ! furthest separation between "close" locations
integer, allocatable :: lon_offset(:, :) ! (nlat, nlat), lon box indices searched (nlat 2x IS correct)
integer, allocatable :: loc_box(:) ! (nloc), List of loc indices in boxes
integer, allocatable :: count(:, :) ! (nlon, nlat), # of loc in each box
integer, allocatable :: start(:, :) ! (nlon, nlat), Start of list of loc in this box
real(r8) :: bot_lat, top_lat ! Bottom and top latitudes of latitude boxes
real(r8) :: bot_lon, top_lon ! Bottom and top longitudes of longitude boxes
real(r8) :: lon_width, lat_width ! Width of boxes in lon and lat
logical :: lon_cyclic ! Do boxes wraparound in longitude?
end type get_close_type_by_type
! Support more than a single cutoff distance. nt is the count of
! distinct cutoffs, which are selected per specific observation type.
! The map associates the incoming location type with the
! corresponding gtt index. There are only as many close_types
! as there are distinct cutoff distances; if more than one specific
! type has the same cutoff distances they share the type.
type get_close_type
integer :: nt ! number of distinct cutoffs
integer, allocatable :: type_to_cutoff_map(:) ! mapping of types to index
type(get_close_type_by_type), allocatable :: gtt(:) ! array of close_types by type
end type get_close_type
! Horizontal localization/cutoff values are passed in by the caller.
! The Vertical normalization values are globals; are specified by namelist
! here, and apply to all specific types unless a 'special list' is also specified
! that overrides the default values.
logical :: has_special_vertical_norms = .false.
integer :: num_special_vert_norms = 0
integer :: location_vertical_localization_coord = 0
! Some calls include information about the type or kind of the location.
! If the location refers to an observation it is possible to have both
! a specific type and a generic kind. If the location refers to a
! state vector item, it only has a generic kind. Some routines have
! a specific type and a generic kind as arguments; look carefully at
! the argument names before using.
type(random_seq_type) :: ran_seq
logical :: ran_seq_init = .false.
logical, save :: module_initialized = .false.
character(len = 512) :: msgstring, msgstring1, msgstring2
! Global storage for vertical distance normalization factors
! The 4 below is for the 4 vertical units (pressure, level, height,
! scale height). undefined and surface don't need vert distances.
! NOTE: Code that uses VERT_TYPE_COUNT depends on pressure, level,
! height, and scale height having actual values between 1 and 4, or
! this code will break.
integer, parameter :: VERT_TYPE_COUNT = 4
real(r8) :: vert_normalization(VERT_TYPE_COUNT)
real(r8), allocatable :: per_type_vert_norm(:,:) ! if doing per-type
! Global storage for fast approximate sin and cosine lookups
! PAR For efficiency for small cases might want to fill tables as needed
! Also these could be larger for more accurate computations, if needed.
! 630 is 2 * PI rounded up, times 100.
real(r8), parameter :: SINCOS_DELTA = 100.0_r8
integer, parameter :: SINCOS_LIMIT = 630
real(r8), parameter :: ACOS_DELTA = 1000.0_r8
integer, parameter :: ACOS_LIMIT = 1000
real(r8) :: my_sin(-SINCOS_LIMIT:SINCOS_LIMIT)
real(r8) :: my_cos(-SINCOS_LIMIT:SINCOS_LIMIT)
real(r8) :: my_acos(-ACOS_LIMIT:ACOS_LIMIT)
! Tolerance for the top latitude boundary test. All locations which
! are located on a box boundary are added to the bin on the larger
! side of the boundary (e.g. for north-south latitudes, it rounds
! towards the north). But when a value falls exactly on the edge
! of the last box, technically it is inside the region but would be
! rounded up and outside the region unless handled specially.
! this tolerance below is used to determine if a value is within
! the range of the last box boundary and if so, the location is
! included in the last box. In particular, for global grids this
! preserves locations which are at exactly 90.0 degrees latitude.
real(r8), parameter :: EDGE_TOLERANCE = 100.0_r8 * epsilon(0.0_r8)
! Option for verification using exhaustive search
logical :: COMPARE_TO_CORRECT = .false. ! normally false
!-----------------------------------------------------------------
! Namelist with default values
! horiz_dist_only == .true. -> Only the great circle horizontal distance is
! computed in get_dist.
! horiz_dist_only == .false. -> Square root of sum of squared horizontal and
! normalized vertical dist computed in get_dist
! vert_normalization_pressure -> Number pascals that give a distance equivalent
! to one radian in horizontal
! vert_normalization_height -> Number meters that give a distance equivalent
! to one radian in horizontal
! vert_normalization_level -> Number levels that give a distance equivalent
! to one radian in horizontal
! vert_normalization_scale_height -> Number scale heights that give a distance
! equivalent to one radian in horizontal
! approximate_distance -> Use a faster table lookup for the trig math.
! Works well for global models and large areas,
! and improves performance. For smaller regions
! might introduce banding, so leave .false.
! nlon -> Number longitude boxes for get_close
! nlon MUST BE ODD
! nlat -> Number latitude boxes for get_close
! output_box_info -> Useful for debugging performance problems.
! print_box_level -> How much data to print out.
! special_vert_normalization_obs_types -> Which obs types to modify the default vert
! normalization values
! special_vert_normalization_pressure -> must give all 4 values for each type listed
! special_vert_normalization_heights
! special_vert_normalization_levels
! special_vert_normalization_scale_heights
logical :: horiz_dist_only = .true.
real(r8) :: vert_normalization_pressure = 100000.0_r8
real(r8) :: vert_normalization_height = 10000.0_r8
real(r8) :: vert_normalization_level = 20.0_r8
real(r8) :: vert_normalization_scale_height = 5.0_r8
logical :: approximate_distance = .false.
integer :: nlon = 71
integer :: nlat = 36
logical :: output_box_info = .false.
integer :: print_box_level = 0
!>@todo make these allocatable instead of fixed length
integer, parameter :: MAX_ITEMS = 500
character(len=OBSTYPELENGTH) :: special_vert_normalization_obs_types(MAX_ITEMS)
real(r8) :: special_vert_normalization_pressures(MAX_ITEMS)
real(r8) :: special_vert_normalization_heights(MAX_ITEMS)
real(r8) :: special_vert_normalization_levels(MAX_ITEMS)
real(r8) :: special_vert_normalization_scale_heights(MAX_ITEMS)
namelist /location_nml/ horiz_dist_only, vert_normalization_pressure, &
vert_normalization_height, vert_normalization_level, &
vert_normalization_scale_height, approximate_distance, nlon, nlat, &
output_box_info, print_box_level, &
special_vert_normalization_obs_types, special_vert_normalization_pressures, &
special_vert_normalization_heights, special_vert_normalization_levels, &
special_vert_normalization_scale_heights
!-----------------------------------------------------------------
interface operator(==); module procedure loc_eq; end interface
interface operator(/=); module procedure loc_ne; end interface
interface set_location
module procedure set_location_single
module procedure set_location_array
end interface set_location
contains
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
! basic location routines
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
! things which need doing exactly once.
subroutine initialize_module()
integer :: iunit, io, i, k, typecount, type_index
if (module_initialized) return
module_initialized = .true.
! give these initial values before reading them from the namelist.
special_vert_normalization_obs_types(:) = 'null'
special_vert_normalization_pressures(:) = missing_r8
special_vert_normalization_heights(:) = missing_r8
special_vert_normalization_levels(:) = missing_r8
special_vert_normalization_scale_heights(:) = missing_r8
! Read the namelist entry
call find_namelist_in_file("input.nml", "location_nml", iunit)
read(iunit, nml = location_nml, iostat = io)
call check_namelist_read(iunit, io, "location_nml")
! Write the namelist values to the log file
if(do_nml_file()) write(nmlfileunit, nml=location_nml)
if(do_nml_term()) write( * , nml=location_nml)
call error_handler(E_MSG, 'location_mod:', 'using code with optimized cutoffs', &
source)
! Simple error checking for nlon, nlat to be sure we can use them for allocations.
if (nlon < 1 .or. nlat < 1) then
write(msgstring, '(A,I5,A,I5)') 'from the namelist, nlon is ', nlon, '; nlat is ', nlat
write(msgstring2, '(A)') 'both must be >= 1, and nlon must be odd'
call error_handler(E_ERR, 'location_mod', msgstring, source, text2=msgstring2)
endif
! One additional restriction; number of longitudes, nlon, for get_close
if((nlon / 2) * 2 == nlon) then
write(msgstring, '(A,I5,A)') 'nlon is ', nlon, '. Must be odd'
call error_handler(E_ERR, 'location_mod', msgstring, source)
endif
! Copy the normalization factors in the vertical into an array
vert_normalization(VERTISLEVEL) = vert_normalization_level
vert_normalization(VERTISPRESSURE) = vert_normalization_pressure
vert_normalization(VERTISHEIGHT) = vert_normalization_height
vert_normalization(VERTISSCALEHEIGHT) = vert_normalization_scale_height
! If the user has set a namelist to make different vertical normalization factors
! when computing the localization distances, allocate and set that up here.
! It overrides the defaults above.
if (special_vert_normalization_obs_types(1) /= 'null' .or. &
special_vert_normalization_pressures(1) /= missing_r8 .or. &
special_vert_normalization_heights(1) /= missing_r8 .or. &
special_vert_normalization_levels(1) /= missing_r8 .or. &
special_vert_normalization_scale_heights(1) /= missing_r8) then
! FIXME: add code to check for mismatched length lists. are we going to force
! users to specify all 4 values for any obs type that is not using the defaults?
typecount = get_num_types_of_obs() ! ignore function name, this is specific type count
allocate(per_type_vert_norm(VERT_TYPE_COUNT, typecount))
! Set the defaults for all specific types not listed in the special list
per_type_vert_norm(VERTISLEVEL, :) = vert_normalization_level
per_type_vert_norm(VERTISPRESSURE, :) = vert_normalization_pressure
per_type_vert_norm(VERTISHEIGHT, :) = vert_normalization_height
per_type_vert_norm(VERTISSCALEHEIGHT, :) = vert_normalization_scale_height
! Go through special-treatment observation kinds, if any.
num_special_vert_norms = 0
k = 0
do i = 1, MAX_ITEMS
if(special_vert_normalization_obs_types(i) == 'null') exit
k = k + 1
enddo
num_special_vert_norms = k
if (num_special_vert_norms > 0) has_special_vertical_norms = .true.
do i = 1, num_special_vert_norms
type_index = get_index_for_type_of_obs(special_vert_normalization_obs_types(i))
if (type_index < 0) then
write(msgstring, *) 'unrecognized TYPE_ in the special vertical normalization namelist:'
call error_handler(E_ERR,'location_mod:', msgstring, source, &
text2=trim(special_vert_normalization_obs_types(i)))
endif
per_type_vert_norm(VERTISLEVEL, type_index) = special_vert_normalization_levels(i)
per_type_vert_norm(VERTISPRESSURE, type_index) = special_vert_normalization_pressures(i)
per_type_vert_norm(VERTISHEIGHT, type_index) = special_vert_normalization_heights(i)
per_type_vert_norm(VERTISSCALEHEIGHT, type_index) = special_vert_normalization_scale_heights(i)
enddo
if (any(per_type_vert_norm == missing_r8)) then
write(msgstring, *) 'one or more special vertical normalization values is uninitialized.'
call error_handler(E_ERR,'location_mod:', &
'special vert normalization value namelist requires all 4 values per type', &
source, text2=trim(msgstring))
endif
endif
! if the namelist control says we are only basing the localization on
! distances in the horizontal, log that in the dart log file.
if (horiz_dist_only) then
call error_handler(E_MSG,'location_mod:', &
'Ignoring vertical separation when computing distances; horizontal distances only', &
source)
else
call error_handler(E_MSG,'location_mod:', &
'Including vertical separation when computing distances:', source)
write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', vert_normalization_pressure
call error_handler(E_MSG,'location_mod:',msgstring,source)
write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', vert_normalization_height
call error_handler(E_MSG,'location_mod:',msgstring,source)
write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', vert_normalization_level
call error_handler(E_MSG,'location_mod:',msgstring,source)
write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', vert_normalization_scale_height
call error_handler(E_MSG,'location_mod:',msgstring,source)
if (allocated(per_type_vert_norm)) then
typecount = get_num_types_of_obs() ! ignore function name, this is specific type count
do i = 1, typecount
if ((per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) .or. &
(per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) .or. &
(per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) .or. &
(per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height)) then
write(msgstring,'(2A)') 'Altering default vertical normalization for type ', trim(get_name_for_type_of_obs(i))
call error_handler(E_MSG,'location_mod:',msgstring,source)
if (per_type_vert_norm(VERTISPRESSURE, i) /= vert_normalization_pressure) then
write(msgstring,'(A,f17.5)') ' # pascals ~ 1 horiz radian: ', &
per_type_vert_norm(VERTISPRESSURE, i)
call error_handler(E_MSG,'location_mod:',msgstring,source)
endif
if (per_type_vert_norm(VERTISHEIGHT, i) /= vert_normalization_height) then
write(msgstring,'(A,f17.5)') ' # meters ~ 1 horiz radian: ', &
per_type_vert_norm(VERTISHEIGHT, i)
call error_handler(E_MSG,'location_mod:',msgstring,source)
endif
if (per_type_vert_norm(VERTISLEVEL, i) /= vert_normalization_level) then
write(msgstring,'(A,f17.5)') ' # model levels ~ 1 horiz radian: ', &
per_type_vert_norm(VERTISLEVEL, i)
call error_handler(E_MSG,'location_mod:',msgstring,source)
endif
if (per_type_vert_norm(VERTISSCALEHEIGHT, i) /= vert_normalization_scale_height) then
write(msgstring,'(A,f17.5)') ' # scale heights ~ 1 horiz radian: ', &
per_type_vert_norm(VERTISSCALEHEIGHT, i)
call error_handler(E_MSG,'location_mod:',msgstring,source)
endif
endif
enddo
endif
endif
! Set up a lookup table for cos and sin for approximate but fast distances
! Don't worry about rounding errors as long as one gives more weight
! Really only need tables half this size, too (sin from -pi/2 to pi/2, cos only +)
if(approximate_distance) then
do i = -SINCOS_LIMIT, SINCOS_LIMIT
my_cos(i) = cos(i / SINCOS_DELTA)
my_sin(i) = sin(i / SINCOS_DELTA)
end do
do i = -ACOS_LIMIT, ACOS_LIMIT
my_acos(i) = acos(i / ACOS_DELTA)
end do
call error_handler(E_MSG,'location_mod:', &
'Using table-lookup approximation for distance computations', source)
endif
end subroutine initialize_module
!----------------------------------------------------------------------------
function get_dist(loc1, loc2, type1, kind2, no_vert)
! returns the distance between 2 locations in units of radians.
! Distance depends on only horizontal distance or both horizontal and
! vertical distance.
! The choice is determined by horiz_dist_only and the which_vert of loc1.
! May want to allow return of some joint distance in the long run?
! Or just a distance that is a function of all 3 things.
! The namelist controls whether default computations use just horizontal distance.
! However, this behavior can be over-ridden by the no_vert optional argument.
! If set to false, this will always do full 3d distance if possible. If set to
! true it will never do the full 3d distance. At present asking to do a vertical
! distance computation for incompatible vertical location units results
! in a fatal error unless one of the vertical units is UNDEFINED.
! CHANGE from previous versions: the 3rd argument is now a specific type
! (e.g. RADIOSONDE_TEMPERATURE, AIRCRAFT_SPECIFIC_HUMIDITY) associated
! with loc1, while the 4th argument is a generic kind (QTY_TEMPERATURE,
! QTY_U_WIND_COMPONENT) associated with loc2.
! The type and kind are part of the interface in case user-code wants to do
! a more sophisticated distance calculation based on the base type or target
! kind. In the usual case this code still doesn't use the kind/type, but
! it does require at least type1 if using per-type vertical normalization.
type(location_type), intent(in) :: loc1, loc2
integer, optional, intent(in) :: type1, kind2
logical, optional, intent(in) :: no_vert
real(r8) :: get_dist
real(r8) :: lon_dif, vert_dist, rtemp
integer :: lat1_ind, lat2_ind, lon_ind, temp ! indexes into lookup tables
logical :: comp_h_only
if ( .not. module_initialized ) call initialize_module()
! Begin with the horizontal distance
! Compute great circle path shortest route between two points
lon_dif = loc1%lon - loc2%lon
if(approximate_distance) then
! Option 1: Use table lookup; faster but less accurate
lat1_ind = int(loc1%lat*SINCOS_DELTA)
lat2_ind = int(loc2%lat*SINCOS_DELTA)
lon_ind = int(lon_dif *SINCOS_DELTA)
temp = int(ACOS_DELTA * (my_sin(lat2_ind) * my_sin(lat1_ind) + &
my_cos(lat2_ind) * my_cos(lat1_ind) * my_cos(lon_ind)))
get_dist = my_acos(temp)
else
! Option 2: Use pre-defined trig functions: accurate but slow
! First 2 ifs avoids round-off error that can kill acos;
if(abs(loc1%lat) >= PI/2.0_r8 .or. abs(loc2%lat) >= PI/2.0_r8 .or. &
lon_dif == 0.0_r8) then
get_dist = abs(loc2%lat - loc1%lat)
else
! This test is for apparent roundoff error which may be a result of
! running r8 == r4.
rtemp = sin(loc2%lat) * sin(loc1%lat) + &
cos(loc2%lat) * cos(loc1%lat) * cos(lon_dif)
if (rtemp < -1.0_r8) then
get_dist = PI
else if (rtemp > 1.0_r8) then
get_dist = 0.0_r8
else
get_dist = acos(rtemp)
endif
endif
endif
! Now compute a vertical distance if requested. Highest priority is
! the optional no_vert argument, so test it first.
if(present(no_vert)) then
comp_h_only = no_vert
! Namelist horizontal only has second highest priority
else
comp_h_only = horiz_dist_only
endif
! If which_vert has no vertical definition for either location do only horizontal
if(loc1%which_vert == VERTISUNDEF .or. loc2%which_vert == VERTISUNDEF) comp_h_only = .true.
! If both verts are surface, do only horizontal
if(loc1%which_vert == VERTISSURFACE .and. loc2%which_vert == VERTISSURFACE) comp_h_only = .true.
! Add in vertical component if required
if(.not. comp_h_only) then
! Vert distance can only be done for like vertical locations units
if(loc1%which_vert /= loc2%which_vert) then
write(msgstring,*)'loc1%which_vert (',loc1%which_vert, &
') /= loc2%which_vert (',loc2%which_vert,')'
call error_handler(E_MSG, 'get_dist', msgstring, source)
call write_location(logfileunit,loc1)
call write_location(logfileunit,loc2)
call write_location(0,loc1)
call write_location(0,loc2)
call error_handler(E_ERR, 'get_dist', &
'Dont know how to compute vertical distance for unlike vertical coordinates', &
source)
endif
! Compute the difference and divide by the appropriate normalization factor
! Normalization factor computes relative distance in vertical compared to one radian
! This is new - if per-type localization distances given, use the specific type of loc1
! to determine the vertical mapping distance. it defaults to the 4 standard ones,
! but can be specified separately if desired.
! note that per-type vertical conversion factors are a global here, and were set
! by namelist values. they apply to all calls to get_dist() based on the obs type.
if (allocated(per_type_vert_norm)) then
if (.not. present(type1)) then
write(msgstring, *) 'obs type required in get_dist() if doing per-type vertical normalization'
call error_handler(E_MSG, 'get_dist', msgstring, source)
endif
vert_dist = abs(loc1%vloc - loc2%vloc) / per_type_vert_norm(loc1%which_vert, type1)
else
vert_dist = abs(loc1%vloc - loc2%vloc) / vert_normalization(loc1%which_vert)
endif
! Spherical distance shape is computed here, other flavors can be computed
get_dist = sqrt(get_dist**2 + vert_dist**2)
endif
end function get_dist
!---------------------------------------------------------------------------
function loc_eq(loc1,loc2)
! Interface operator used to compare two locations.
! Returns true only if all components are 'the same' to within machine
! precision. There is some debate as to whether or not the vertical
! locations need to be identical if 'VERTISUNDEF' ... hard to peruse
! the code tree to find where this may be affected.
type(location_type), intent(in) :: loc1, loc2
logical :: loc_eq
if ( .not. module_initialized ) call initialize_module()
loc_eq = .false.
! if ( loc1%which_vert /= loc2%which_vert ) return
if ( abs(loc1%lon - loc2%lon) > epsilon(loc1%lon) ) return
if ( abs(loc1%lat - loc2%lat) > epsilon(loc1%lat) ) return
!if ( loc1%which_vert /= VERTISUNDEF ) then
if ( abs(loc1%vloc - loc2%vloc) > epsilon(loc1%vloc) ) return
!endif
loc_eq = .true.
end function loc_eq
!---------------------------------------------------------------------------
function loc_ne(loc1,loc2)
! Interface operator used to compare two locations.
! Returns true if locations are not identical to machine precision.
type(location_type), intent(in) :: loc1, loc2
logical :: loc_ne
if ( .not. module_initialized ) call initialize_module()
loc_ne = (.not. loc_eq(loc1,loc2))
end function loc_ne
!---------------------------------------------------------------------------
function get_location(loc)
! Given a location type (in radians),
! return the longitude, latitude (in degrees) and vertical value
type(location_type), intent(in) :: loc
real(r8), dimension(3) :: get_location
if ( .not. module_initialized ) call initialize_module()
get_location(1) = loc%lon * RAD2DEG
get_location(2) = loc%lat * RAD2DEG
get_location(3) = loc%vloc
end function get_location
!---------------------------------------------------------------------------
function set_location_single(lon, lat, vert_loc, which_vert)
! Puts the given longitude, latitude, and vertical location
! into a location datatype. Arguments to this function are in degrees,
! but the values are stored as radians.
real(r8), intent(in) :: lon, lat
real(r8), intent(in) :: vert_loc
integer, intent(in) :: which_vert
type (location_type) :: set_location_single
if ( .not. module_initialized ) call initialize_module()
if(lon < 0.0_r8 .or. lon > 360.0_r8) then
write(msgstring,*)'longitude (',lon,') is not within range [0,360]'
call error_handler(E_ERR, 'set_location', msgstring, source)
endif
if(lat < -90.0_r8 .or. lat > 90.0_r8) then
write(msgstring,*)'latitude (',lat,') is not within range [-90,90]'
call error_handler(E_ERR, 'set_location', msgstring, source)
endif
set_location_single%lon = lon * DEG2RAD
set_location_single%lat = lat * DEG2RAD
if(which_vert < VERTISUNDEF .or. which_vert == 0 .or. which_vert > VERTISSCALEHEIGHT) then
write(msgstring,*)'which_vert (',which_vert,') must be one of -2, -1, 1, 2, 3, or 4'
call error_handler(E_ERR,'set_location', msgstring, source)
endif
set_location_single%vloc = vert_loc
set_location_single%which_vert = which_vert
end function set_location_single
!----------------------------------------------------------------------------
function set_location_array(list)
! location semi-independent interface routine
! given 4 float numbers, call the underlying set_location routine
real(r8), intent(in) :: list(:)
type (location_type) :: set_location_array
if ( .not. module_initialized ) call initialize_module()
if (size(list) < 4) then
write(msgstring,*)'requires 4 input values'
call error_handler(E_ERR, 'set_location', msgstring, source)
endif
set_location_array = set_location_single(list(1), list(2), list(3), nint(list(4)))
end function set_location_array
!----------------------------------------------------------------------------
function set_location_missing()
! Initialize a location type to indicate the contents are unset.
type (location_type) :: set_location_missing
if ( .not. module_initialized ) call initialize_module()
set_location_missing%lon = MISSING_R8
set_location_missing%lat = MISSING_R8
set_location_missing%vloc = MISSING_R8
set_location_missing%which_vert = MISSING_I
end function set_location_missing
!---------------------------------------------------------------------------
function query_location(loc, attr)
! Returns the value of the attribute
type(location_type), intent(in) :: loc
character(len=*), optional, intent(in) :: attr
real(r8) :: query_location
if ( .not. module_initialized ) call initialize_module()
! Workaround for apparent bug in mac osx intel 10.x fortran compiler.
! Previous code had a 16 byte local character variable which was
! apparently not getting deallocated after this function returned.
! Workaround for a second compiler bug. The PGI 6.1.x compiler
! refused to compile the select statement below when it was
! selectcase(adjustl(attr)). i'm rearranging the code again because
! this particular routine has been so troublesome. i'm removing
! the result(fval) construct, and setting the default to which_vert
! and then overwriting it the case we recognize another quantity.
! if we fall through to the end of the routine, the return value
! is which_vert. (This routine is rapidly becoming my favorite
! problem child. How can such a short piece of code be so troublesome?)
! set the default here, and then only overwrite it if we
! recognize one of the other valid queries.
query_location = real(loc%which_vert, r8) ! this is really an int
if (.not. present(attr)) return
select case(attr)
case ('lat','LAT')
query_location = loc%lat
case ('lon','LON')
query_location = loc%lon
case ('vloc','VLOC')
query_location = loc%vloc
case ('which_vert','WHICH_VERT')
! already set
case default
call error_handler(E_ERR, 'query_location:', &
'Only "lon","lat","vloc","which_vert" are legal attributes to request from location', &
source)
end select
end function query_location
!----------------------------------------------------------------------------
subroutine write_location(locfile, loc, fform, charstring)
! Writes a location to a file.
! most recent change: adding the optional charstring option. if present,
! locfile is ignored, and a pretty-print formatting is done into charstring.
! the locations are converted to lat/lon, and the vert is put into more
! common units (e.g. hPa, km, etc)
integer, intent(in) :: locfile
type(location_type), intent(in) :: loc
character(len = *), intent(in), optional :: fform
character(len = *), intent(out), optional :: charstring
integer :: charlength
logical :: writebuf
character(len=129) :: string1
! 10 format(1x,3(f22.14,1x),i4) ! old
10 format(1X,3(G25.16,1X),I2)
if ( .not. module_initialized ) call initialize_module()
! writing to a file (normal use) or to a character buffer?
writebuf = present(charstring)
! output file; test for ascii or binary, write what's asked, and return
if (.not. writebuf) then
if (ascii_file_format(fform)) then
write(locfile, '(''loc3d'')' )
write(locfile, 10) loc%lon, loc%lat, loc%vloc, loc%which_vert
else
write(locfile) loc%lon, loc%lat, loc%vloc, loc%which_vert
endif
return
endif
! you only get here if you're writing to a buffer and not
! to a file, and you can't have binary format set.
if (.not. ascii_file_format(fform)) then
call error_handler(E_ERR, 'write_location', &
'Cannot use string buffer with binary format', source)
endif
! format the location to be more human-friendly; meaning
! degrees instead of radians, and kilometers for height,
! hectopascals instead of pascals for pressure, etc.
! this must be the sum of the longest of the formats below.
charlength = 72
if (len(charstring) < charlength) then
write(msgstring, *) 'charstring buffer must be at least ', charlength, ' chars long'
call error_handler(E_ERR, 'write_location', msgstring, source)
endif
! format the horizontal into a temp string
write(string1, '(A,F12.8,1X,F12.8,1X,A)') 'Lon/Lat(deg): ', loc%lon*RAD2DEG, &
loc%lat*RAD2DEG, ' Vert:'
! then pretty up the vertical choices, trying to get them to line up in
! case the caller is listing out locations with different vert units.
! concatinate the vertical on the end of the horizontal and put it all
! into the return string.
select case (loc%which_vert)
case (VERTISUNDEF)
write(charstring, '(A,A)') trim(string1), ' Undefined'
case (VERTISSURFACE)
write(charstring, '(A,F13.5,A)') trim(string1), loc%vloc, ' surface (m)'
case (VERTISLEVEL)
write(charstring, '(A,F13.6,A)') trim(string1), loc%vloc, ' level'
case (VERTISPRESSURE)
write(charstring, '(A,F13.7,A)') trim(string1), loc%vloc / 100.0_r8, ' hPa'
case (VERTISHEIGHT)
write(charstring, '(A,F13.7,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
case (VERTISSCALEHEIGHT)
write(charstring, '(A,F13.7,A)') trim(string1), loc%vloc, ' scale ht'
case default
write(msgstring, *) 'unrecognized key for vertical type: ', loc%which_vert
call error_handler(E_ERR, 'write_location', msgstring, source)
end select
end subroutine write_location
!----------------------------------------------------------------------------
function read_location(locfile, fform)
! Reads a location from a file that was written by write_location.
! See write_location for additional discussion.
integer, intent(in) :: locfile
character(len = *), intent(in), optional :: fform
type(location_type) :: read_location
character(len=5) :: header
if ( .not. module_initialized ) call initialize_module()
if (ascii_file_format(fform)) then
read(locfile, '(a5)' ) header
if(header /= 'loc3d') then
write(msgstring,*)'Expected location header "loc3d" in input file, got ', header
call error_handler(E_ERR, 'read_location', msgstring, source)
endif
! Now read the location data value
read(locfile, *)read_location%lon, read_location%lat, &
read_location%vloc, read_location%which_vert
else
read(locfile)read_location%lon, read_location%lat, &
read_location%vloc, read_location%which_vert
endif
end function read_location
!--------------------------------------------------------------------------
subroutine interactive_location(location, set_to_default)
! Allows for interactive input of a location. Also gives option of selecting
! a uniformly distributed random location.
type(location_type), intent(out) :: location
logical, intent(in), optional :: set_to_default
real(r8) :: lon, lat, minlon, maxlon, minlat, maxlat
if ( .not. module_initialized ) call initialize_module()
! If set_to_default is true, then just zero out and return
if(present(set_to_default)) then
if(set_to_default) then
location%lon = 0.0
location%lat = 0.0
location%vloc = 0.0
location%which_vert = 0 ! note that 0 isn't a valid vert type
return
endif
endif
write(*, *)'Vertical coordinate options'
write(*, *)VERTISUNDEF,' --> vertical coordinate undefined'
write(*, *)VERTISSURFACE,' --> surface'
write(*, *)VERTISLEVEL,' --> model level'
write(*, *)VERTISPRESSURE,' --> pressure'
write(*, *)VERTISHEIGHT,' --> height'
write(*, *)VERTISSCALEHEIGHT,' --> scale height'
100 read(*, *) location%which_vert
if(location%which_vert == VERTISLEVEL ) then
write(*, *) 'Vertical coordinate model level'
read(*, *) location%vloc
else if(location%which_vert == VERTISPRESSURE ) then
write(*, *) 'Vertical coordinate Pressure (in hPa)'
read(*, *) location%vloc
location%vloc = 100.0 * location%vloc
else if(location%which_vert == VERTISHEIGHT ) then
write(*, *) 'Vertical coordinate height (in meters)'
read(*, *) location%vloc
else if(location%which_vert == VERTISSURFACE ) then
! most 3d sphere users want height in meters, not pressure.
! original code asked for pressure:
!write(*, *) 'Vertical coordinate surface pressure (in hPa)'
!location%vloc = 100.0 * location%vloc ! only applies to pressure
write(*, *) 'Vertical coordinate height'
read(*, *) location%vloc
else if(location%which_vert == VERTISSCALEHEIGHT ) then
write(*, *) 'Vertical coordinate scale height (-ln(p/ps))'
read(*, *) location%vloc
else if(location%which_vert == VERTISUNDEF ) then
! a valid floating point value, but should be unused.
location%vloc = MISSING_R8
else
write(*, *) 'Wrong choice of which_vert try again between ',VERTISUNDEF, &
' and ',VERTISHEIGHT
go to 100
end if
write(*, *) 'Input longitude: value 0 to 360.0 or a negative number for '
write(*, *) 'Uniformly distributed random location in the horizontal'
read(*, *) lon
do while(lon > 360.0_r8)
write(*, *) 'Input value greater than 360.0 is illegal, please try again'
read(*, *) lon
end do
if(lon < 0.0_r8) then
! Need to make sure random sequence is initialized
if(.not. ran_seq_init) then
call init_random_seq(ran_seq)
ran_seq_init = .TRUE.
endif
write(*, *) 'Input minimum longitude (0 to 360.0)'
read(*, *) minlon
minlon = minlon * DEG2RAD
write(*, *) 'Input maximum longitude (0 to 360.0)'
read(*, *) maxlon
maxlon = maxlon * DEG2RAD
! Longitude is random from minlon to maxlon
location%lon = random_uniform(ran_seq) * (maxlon-minlon) + minlon
write(*, *) 'Input minimum latitude (-90.0 to 90.0)'
read(*, *) minlat
minlat = sin(minlat * DEG2RAD)
write(*, *) 'Input maximum latitude (-90.0 to 90.0)'
read(*, *) maxlat
maxlat = sin(maxlat * DEG2RAD)
! Latitude must be area weighted
location%lat = asin(random_uniform(ran_seq) * (maxlat-minlat) + minlat)
write(*, *) 'random location is ', location%lon / DEG2RAD, &
location%lat / DEG2RAD
else
write(*, *) 'Input latitude: value -90.0 to 90.0'
read(*, *) lat
do while(lat < -90.0_r8 .or. lat > 90.0_r8)
write(*, *) 'Input value < -90.0 or > 90.0 is illegal, please try again'
read(*, *) lat
end do
location%lat = lat*DEG2RAD
location%lon = lon*DEG2RAD
end if
end subroutine interactive_location
!----------------------------------------------------------------------------
function is_location_in_region(loc, minl, maxl)
! Returns true if the given location is inside the rectangular
! region defined by minl as the lower left, maxl the upper right.
! test is inclusive; values on the edges are considered inside.
! Periodic in longitude (box can cross the 2PI -> 0 line)
logical :: is_location_in_region
type(location_type), intent(in) :: loc, minl, maxl
if ( .not. module_initialized ) call initialize_module()
! maybe could use VERTISUNDEF in the minl and maxl args to indicate
! we want to test only in horizontal? and if not, vtypes must match?
!if ( (minl%which_vert /= maxl%which_vert) .or. &
! ((minl%which_vert /= loc%which_vert).and.(minl%which_vert /= VERTISUNDEF))) then
! write(msgstring,*)'which_vert (',loc%which_vert,') must be same in all args'
! call error_handler(E_ERR, 'is_location_in_region', msgstring, source)
!endif
! assume failure and return as soon as we are confirmed right.
! set to success only at the bottom after all tests have passed.
is_location_in_region = .false.