-
Notifications
You must be signed in to change notification settings - Fork 145
/
model_mod.f90
6441 lines (4937 loc) · 235 KB
/
model_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
!
! $Id$
module model_mod
! MPAS ocean model interface to the DART data assimilation system.
! Code in this module is compiled with the DART executables. It isolates
! all information about the MPAS grids, model variables, and other details.
! There are a set of 16 subroutine interfaces that are required by DART;
! these cannot be changed. Additional public routines in this file can
! be used by converters and utilities and those interfaces can be anything
! that is useful to other pieces of code.
! Units on everything are MKS:
!
! u (velocity): meter / second
! h (depth) : meter
! rho (density) : kilograms / meter^3
! temperature : *potential temperature* degrees C
! salinity : PSU
!
! Note: the 'temperature' variable is *potential* temperature.
! Routines in other modules that are used here.
use types_mod, only : r4, r8, digits12, SECPERDAY, MISSING_R8, &
rad2deg, deg2rad, PI, MISSING_I, obstypelength
use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
print_time, print_date, set_calendar_type, &
operator(*), operator(+), operator(-), &
operator(>), operator(<), operator(/), &
operator(/=), operator(<=)
use location_mod, only : location_type, get_dist, query_location, &
get_close_maxdist_init, get_close_type, &
set_location, get_location, horiz_dist_only, &
write_location, &
vert_is_undef, VERTISUNDEF, &
vert_is_surface, VERTISSURFACE, &
vert_is_level, VERTISLEVEL, &
vert_is_pressure, VERTISPRESSURE, &
vert_is_height, VERTISHEIGHT, &
vert_is_scale_height, VERTISSCALEHEIGHT, &
get_close_obs_init, get_close_obs_destroy, &
loc_get_close_obs => get_close_obs
use xyz_location_mod, only : xyz_location_type, xyz_get_close_maxdist_init, &
xyz_get_close_type, xyz_set_location, xyz_get_location, &
xyz_get_close_obs_init, xyz_get_close_obs_destroy, &
xyz_find_nearest
use utilities_mod, only : register_module, error_handler, &
E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
nc_check, do_output, to_upper, nmlfileunit, &
find_namelist_in_file, check_namelist_read, &
open_file, file_exist, find_textfile_dims, &
file_to_text, close_file, do_nml_file, do_nml_term
use obs_kind_mod, only : get_index_for_quantity, &
get_name_for_quantity, &
QTY_VERTICAL_VELOCITY, &
QTY_POTENTIAL_TEMPERATURE, &
QTY_TEMPERATURE, &
QTY_SALINITY, &
QTY_DRY_LAND, &
QTY_EDGE_NORMAL_SPEED, &
QTY_U_CURRENT_COMPONENT, &
QTY_V_CURRENT_COMPONENT, &
QTY_SEA_SURFACE_HEIGHT, &
QTY_SEA_SURFACE_PRESSURE, &
QTY_TRACER_CONCENTRATION
use mpi_utilities_mod, only: my_task_id
use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
! netcdf modules
use typesizes
use netcdf
! RBF (radial basis function) modules, donated by LANL. currently deprecated
! in this version. they did the job but not as well as other techniques and
! at a much greater execution-time code. they were used to interpolate
! values at arbitrary locations, not just at cell centers. with too small
! a set of basis points, the values were discontinuous at cell boundaries;
! with too many the values were too smoothed out. we went back to
! barycentric interpolation in triangles formed by the three cell centers
! that enclosed the given point.
use get_geometry_mod
use get_reconstruct_mod
implicit none
private
! these routines must be public and you cannot change
! the arguments - they will be called *from* the DART code.
public :: get_model_size, &
adv_1step, &
get_state_meta_data, &
model_interpolate, &
get_model_time_step, &
static_init_model, &
end_model, &
init_time, &
init_conditions, &
nc_write_model_atts, &
nc_write_model_vars, &
pert_model_state, &
get_close_maxdist_init, &
get_close_obs_init, &
get_close_obs, &
ens_mean_for_model
! generally useful routines for various support purposes.
! the interfaces here can be changed as appropriate.
public :: get_model_analysis_filename, &
get_grid_definition_filename, &
analysis_file_to_statevector, &
statevector_to_analysis_file, &
get_analysis_time, &
write_model_time, &
get_grid_dims, &
print_variable_ranges
! version controlled file description for error handling, do not edit
character(len=256), parameter :: source = &
"$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate = "$Date$"
! module global storage; maintains values between calls, accessible by
! any subroutine
character(len=256) :: string1, string2, string3
logical, save :: module_initialized = .false.
! Real (physical) constants as defined exactly in MPAS.
! redefined here for consistency with the model.
real(r8), parameter :: rgas = 287.0_r8
real(r8), parameter :: cp = 1003.0_r8
real(r8), parameter :: cv = 716.0_r8
real(r8), parameter :: p0 = 100000.0_r8
real(r8), parameter :: rcv = rgas/(cp-rgas)
! earth radius; needed to convert lat/lon to x,y,z cartesian coords.
! FIXME: the example ocean files has a global attr with 6371220.0
! the actual mpas code may have hardwired values (it did in the atmosphere)
! need to check what is really going on.
real(r8), parameter :: radius = 6371229.0 ! meters
! roundoff error
real(r8), parameter :: roundoff = 1.0e-12_r8
! Storage for a random sequence for perturbing a single initial state
type(random_seq_type) :: random_seq
! Structure for computing distances to cell centers, and assorted arrays
! needed for the get_close code.
type(xyz_get_close_type) :: cc_gc
type(xyz_location_type), allocatable :: cell_locs(:)
! not part of the namelist because in the ocean it's not clear
! that we need to worry about log pressure in the vertical.
! if this makes a difference, set this to .true.
logical :: log_p_vert_interp = .false. ! if true, interpolate vertical pressure in log space
! variables which are in the module namelist
integer :: vert_localization_coord = VERTISHEIGHT
integer :: assimilation_period_days = 0
integer :: assimilation_period_seconds = 60
real(r8) :: model_perturbation_amplitude = 0.0001 ! tiny amounts
logical :: output_state_vector = .false. ! output prognostic variables (if .false.)
integer :: debug = 0 ! turn up for more and more debug messages
integer :: xyzdebug = 0
character(len=32) :: calendar = 'Gregorian'
character(len=256) :: model_analysis_filename = 'mpas_analysis.nc'
character(len=256) :: grid_definition_filename = 'mpas_analysis.nc'
! if .false. use U/V reconstructed winds tri interp at centers for wind forward ops
! if .true. use edge normal winds (u) with RBF functs for wind forward ops
logical :: use_u_for_wind = .false.
! if using rbf, options 1,2,3 control how many points go into the rbf.
! larger numbers use more points from further away
integer :: use_rbf_option = 2
! if .false. edge normal winds (u) should be in state vector and are written out directly
! if .true. edge normal winds (u) are updated based on U/V reconstructed winds
logical :: update_u_from_reconstruct = .true.
! only if update_u_from_reconstruct is true,
! if .false. use the cell center u,v reconstructed values to update edge winds
! if .true., read in the original u,v winds, compute the increments after the
! assimilation, and use only the increments to update the edge winds
logical :: use_increments_for_u_update = .true.
namelist /model_nml/ &
model_analysis_filename, &
grid_definition_filename, &
vert_localization_coord, &
assimilation_period_days, &
assimilation_period_seconds, &
model_perturbation_amplitude, &
calendar, &
debug, &
xyzdebug, &
use_u_for_wind, &
use_rbf_option, &
update_u_from_reconstruct, &
use_increments_for_u_update
! DART state vector contents are specified in the input.nml:&mpas_vars_nml namelist.
integer, parameter :: max_state_variables = 80
integer, parameter :: num_state_table_columns = 2
integer, parameter :: num_bounds_table_columns = 4
character(len=NF90_MAX_NAME) :: mpas_state_variables(max_state_variables * num_state_table_columns ) = ' '
character(len=NF90_MAX_NAME) :: mpas_state_bounds(num_bounds_table_columns, max_state_variables ) = ' '
character(len=NF90_MAX_NAME) :: variable_table(max_state_variables, num_state_table_columns )
namelist /mpas_vars_nml/ mpas_state_variables, mpas_state_bounds
! FIXME: this shouldn't be a global. the progvar array
! should be allocated at run time and nfields should be part
! of a larger derived type that includes nfields and an array
! of progvartypes.
integer :: nfields
! Everything needed to describe a variable
type progvartype
private
character(len=NF90_MAX_NAME) :: varname
character(len=NF90_MAX_NAME) :: long_name
character(len=NF90_MAX_NAME) :: units
character(len=NF90_MAX_NAME), dimension(NF90_MAX_VAR_DIMS) :: dimname
integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens
integer :: xtype ! netCDF variable type (NF90_double, etc.)
integer :: numdims ! number of dims - excluding TIME
integer :: numvertical ! number of vertical levels in variable
integer :: numcells ! number of horizontal locations (cell centers)
integer :: numedges ! number of horizontal locations (edges for velocity components)
logical :: ZonHalf ! vertical coordinate has dimension nVertLevels
integer :: varsize ! prod(dimlens(1:numdims))
integer :: index1 ! location in dart state vector of first occurrence
integer :: indexN ! location in dart state vector of last occurrence
integer :: dart_kind
character(len=obstypelength) :: kind_string
logical :: clamping ! does variable need to be range-restricted before
real(r8) :: range(2) ! being stuffed back into MPAS analysis file.
logical :: out_of_range_fail ! is out of range fatal if range-checking?
end type progvartype
type(progvartype), dimension(max_state_variables) :: progvar
! Grid parameters - the values will be read from an mpas analysis file.
integer :: nCells = -1 ! Total number of cells making up the grid
integer :: nVertices = -1 ! Unique points in grid that are corners of cells
integer :: nEdges = -1 ! Straight lines between vertices making up cells
integer :: maxEdges = -1 ! Largest number of edges a cell can have
integer :: nVertLevels = -1 ! Vertical levels; count of vert cell centers
integer :: nVertLevelsP1 = -1 ! Vert levels plus 1; count of vert cell faces
integer :: vertexDegree = -1 ! Max number of cells/edges that touch any vertex
! scalar grid positions
! FIXME: we read in a lot of metadata about the grids. if space becomes an
! issue we could consider reading in only the x,y,z arrays for all the items
! plus the radius, and then compute the lat/lon for locations needed by
! get_state_meta_data() on demand. most of the computations we need to do
! are actually easier in xyz coords (no issue with poles).
! FIXME: it may be desirable to read in xCell(:), yCell(:), zCell(:)
! to keep from having to compute them on demand, especially since we
! have converted the radian lat/lon of the cell centers into degrees.
! we have to convert back, then take a few sin and cos to get xyz.
! time/space/accuracy tradeoff here.
real(r8), allocatable :: xVertex(:), yVertex(:), zVertex(:)
real(r8), allocatable :: xEdge(:), yEdge(:), zEdge(:)
real(r8), allocatable :: lonEdge(:) ! edge longitudes (degrees, original radians in file)
real(r8), allocatable :: latEdge(:) ! edge longitudes (degrees, original radians in file)
real(r8), allocatable :: lonCell(:) ! cell center longitudes (degrees, original radians in file)
real(r8), allocatable :: latCell(:) ! cell center latitudes (degrees, original radians in file)
real(r8), allocatable :: zGridFace(:,:) ! geometric depth at cell faces (nVertLevelsP1,nCells)
real(r8), allocatable :: zGridCenter(:,:) ! geometric depth at cell centers (nVertLevels, nCells)
real(r8), allocatable :: zGridEdge(:,:) ! geometric depth at edge centers (nVertLevels, nEdges)
real(r8), allocatable :: zMid(:,:) ! depths at midpoints - may be able to be computed instead
! of requiring it in the input file (save file space). FIXME
real(r8), allocatable :: hZLevel(:) ! layer thicknesses - maybe - FIXME
!real(r8), allocatable :: zEdgeCenter(:,:) ! geometric height at edges faces (nVertLevels ,nEdges)
integer, allocatable :: cellsOnVertex(:,:) ! list of cell centers defining a triangle
integer, allocatable :: verticesOnCell(:,:)
integer, allocatable :: edgesOnCell(:,:) ! list of edges that bound each cell
integer, allocatable :: cellsOnEdge(:,:) ! list of cells that bound each edge
integer, allocatable :: nedgesOnCell(:) ! list of edges that bound each cell
real(r8), allocatable :: edgeNormalVectors(:,:)
! Boundary information might be needed ... regional configuration?
! Read if available.
integer, allocatable :: boundaryCell(:,:) ! logical, cells that are on boundaries
integer, allocatable :: maxLevelEdgeTop(:) !
integer, allocatable :: boundaryEdge(:,:) ! logical, edges that are boundaries
integer, allocatable :: boundaryVertex(:,:) ! logical, vertices that are on boundaries
integer, allocatable :: maxLevelCell(:) ! list of maximum (deepest) level for each cell
real(r8), allocatable :: ens_mean(:) ! needed to convert vertical distances consistently
integer :: model_size ! the state vector length
type(time_type) :: model_timestep ! smallest time to adv model
! useful flags in making decisions when searching for points, etc
logical :: global_grid = .true. ! true = the grid covers the sphere with no holes
logical :: all_levels_exist_everywhere = .true. ! true = cells defined at all levels
logical :: has_edge_u = .false. ! true = has original normal u on edges
logical :: has_uvreconstruct = .false. ! true = has reconstructed at centers
! Do we have any state vector items located on the cell edges?
! If not, avoid reading in or using the edge arrays to save space.
! FIXME: this should be set after looking at the fields listed in the
! namelist which are to be read into the state vector - if any of them
! are located on the edges then this flag should be changed to .true.
! however, the way the code is structured these arrays are allocated
! before the details of field list is examined. since right now the
! only possible field array that is on the edges is the 'u' edge normal
! winds, search specifically for that in the state field list and set
! this based on that. if any other data might be on edges, this routine
! will need to be updated: is_edgedata_in_state_vector()
logical :: data_on_edges = .false.
! currently unused; for a regional model it is going to be necessary to know
! if the grid is continuous around longitudes (wraps in east-west) or not,
! and if it covers either of the poles.
character(len= 64) :: ew_boundary_type, ns_boundary_type
! common names that call specific subroutines based on the arg types
INTERFACE vector_to_prog_var
MODULE PROCEDURE vector_to_1d_prog_var
MODULE PROCEDURE vector_to_2d_prog_var
MODULE PROCEDURE vector_to_3d_prog_var
END INTERFACE
INTERFACE prog_var_to_vector
MODULE PROCEDURE prog_var_1d_to_vector
MODULE PROCEDURE prog_var_2d_to_vector
MODULE PROCEDURE prog_var_3d_to_vector
END INTERFACE
INTERFACE get_analysis_time
MODULE PROCEDURE get_analysis_time_ncid
MODULE PROCEDURE get_analysis_time_fname
END INTERFACE
INTERFACE get_index_range
MODULE PROCEDURE get_index_range_int
MODULE PROCEDURE get_index_range_string
END INTERFACE
!------------------------------------------------
! The regular grid used for triangle interpolation divides the sphere into
! a set of regularly spaced lon-lat boxes. The number of boxes in
! longitude and latitude are set by num_reg_x and num_reg_y. Making the
! number of regular boxes smaller decreases the computation required for
! doing each interpolation but increases the static storage requirements
! and the initialization computation (which seems to be pretty small).
integer, parameter :: num_reg_x = 90, num_reg_y = 90
! The max_reg_list_num controls the size of temporary storage used for
! initializing the regular grid. Two arrays
! of size num_reg_x*num_reg_y*max_reg_list_num are needed. The initialization
! fails and returns an error if max_reg_list_num is too small. A value of
! ??? is sufficient for ???
integer, parameter :: max_reg_list_num = 100
! The triangle interpolation keeps a list of how many and which triangles
! overlap each regular lon-lat box. The number is stored in
! array triangle_num. The allocatable array
! triangle_list lists the uniquen index
! of each overlapping triangle. The entry in
! triangle_start for a given regular lon-lat box indicates
! where the list of triangles begins in the triangle_list.
integer :: triangle_start(num_reg_x, num_reg_y)
integer :: triangle_num (num_reg_x, num_reg_y) = 0
integer, allocatable :: triangle_list(:)
contains
!==================================================================
! All the public REQUIRED interfaces come first - just by convention.
!==================================================================
!------------------------------------------------------------------
subroutine static_init_model()
! Called to do one time initialization of the model.
!
! All the grid information comes from the initialization of
! the dart_model_mod module.
! Local variables - all the important ones have module scope
integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs
character(len=NF90_MAX_NAME) :: varname,dimname
character(len=obstypelength) :: kind_string
integer :: ncid, VarID, numdims, varsize, dimlen
integer :: iunit, io, ivar, i, index1, indexN, iloc, kloc
integer :: ss, dd, z1, m1
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
integer :: cel1, cel2
logical :: both
if ( module_initialized ) return ! only need to do this once.
! Print module information to log file and stdout.
call register_module(source, revision, revdate)
! Since this routine calls other routines that could call this routine
! we'll say we've been initialized pretty dang early.
module_initialized = .true.
! Read the DART namelist for this model
call find_namelist_in_file('input.nml', 'model_nml', iunit)
read(iunit, nml = model_nml, iostat = io)
call check_namelist_read(iunit, io, 'model_nml')
! Record the namelist values used for the run
if (do_nml_file()) write(nmlfileunit, nml=model_nml)
if (do_nml_term()) write( * , nml=model_nml)
! Read the MPAS variable list to populate DART state vector
! Intentionally do not try to dump them to the nml unit because
! they include large character arrays which output pages of white space.
! The routine that reads and parses this namelist will output what
! values it found into the log.
call find_namelist_in_file('input.nml', 'mpas_vars_nml', iunit)
read(iunit, nml = mpas_vars_nml, iostat = io)
call check_namelist_read(iunit, io, 'mpas_vars_nml')
!---------------------------------------------------------------
! Set the time step ... causes mpas namelists to be read.
! Ensures model_timestep is multiple of 'dynamics_timestep'
call set_calendar_type( calendar ) ! comes from model_mod_nml
model_timestep = set_model_time_step()
call get_time(model_timestep,ss,dd) ! set_time() assures the seconds [0,86400)
write(string1,*)'assimilation period is ',dd,' days ',ss,' seconds'
call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
!---------------------------------------------------------------
! 1) get grid dimensions
! 2) allocate space for the grids
! 3) read them from the analysis file
! read_grid_dims() fills in the following module global variables:
! nCells, nVertices, nEdges, maxEdges, nVertLevels, nVertLevelsP1, vertexDegree
call read_grid_dims()
allocate(latCell(nCells), lonCell(nCells))
allocate(zGridFace(nVertLevelsP1, nCells))
allocate(zGridCenter(nVertLevels, nCells))
allocate(cellsOnVertex(vertexDegree, nVertices))
allocate(nEdgesOnCell(nCells))
allocate(edgesOnCell(maxEdges, nCells))
allocate(cellsOnEdge(2, nEdges))
allocate(verticesOnCell(maxEdges, nCells))
allocate(edgeNormalVectors(3, nEdges))
allocate(xVertex(nVertices), yVertex(nVertices), zVertex(nVertices))
! see if U is in the state vector list. if not, don't read in or
! use any of the Edge arrays to save space.
data_on_edges = is_edgedata_in_state_vector(mpas_state_variables)
if(data_on_edges) then
allocate(zGridEdge(nVertLevels, nEdges))
allocate(xEdge(nEdges), yEdge(nEdges), zEdge(nEdges))
allocate(latEdge(nEdges), lonEdge(nEdges))
endif
! this reads in latCell, lonCell, hZLevel, zMid, cellsOnVertex
call get_grid()
! determine which edges are boundaries
! (this requires 'maxLevelEdgeTop' to exist in the restart file)
boundaryEdge(:,1:nEdges) = 1
do iloc=1, nEdges
if(maxLevelEdgeTop(iloc) > 0) then
boundaryEdge(1:maxLevelEdgeTop(iloc),iloc) = 0
endif
end do
! determine which cells are on boundaries
boundaryCell(:,1:nCells) = 0
do kloc=1, nEdges
do iloc=1, nVertLevels
if(boundaryEdge(iloc,kloc) .eq. 1) then
if(cellsOnEdge(1,kloc) > 0) then
boundaryCell(iloc,cellsOnEdge(1,kloc)) = 1
endif
if(cellsOnEdge(2,kloc) > 0) then
boundaryCell(iloc,cellsOnEdge(2,kloc)) = 1
endif
endif
end do
end do
if(data_on_edges) then
! FIXME: This code is supposed to check whether an edge has 2 neighbours or 1 neighbour and then
! compute the depth accordingly. HOWEVER, the array cellsOnEdge does not change with
! depth, but it should as an edge may have 2 neighbour cells at the top but not at depth.
do kloc=1, nEdges
do iloc=1, nVertLevels
cel1 = cellsOnEdge(1,kloc)
cel2 = cellsOnEdge(2,kloc)
if (cel1>0 .and. cel2>0) then
zGridEdge(iloc,kloc) = (zGridCenter(iloc,cel1) + zGridCenter(iloc,cel2))*0.5_r8
else if (cel1>0) then
zGridEdge(iloc,kloc) = zGridCenter(iloc,cel1)
else if (cel2>0) then
zGridEdge(iloc,kloc) = zGridCenter(iloc,cel2)
else !this is bad...
write(string1,*)'Edge ',kloc,' at vertlevel ',iloc,' has no neighbouring cells!'
call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
endif
enddo
enddo
endif
!---------------------------------------------------------------
! Compile the list of model variables to use in the creation
! of the DART state vector. Required to determine model_size.
!
! Verify all variables are in the model analysis file
!
! Compute the offsets into the state vector for the start of each
! different variable type. Requires reading shapes from the model
! analysis file. As long as TIME is the LAST dimension, we're OK.
!
! Record the extent of the data type in the state vector.
call nc_check( nf90_open(trim(model_analysis_filename), NF90_NOWRITE, ncid), &
'static_init_model', 'open '//trim(model_analysis_filename))
call verify_state_variables( mpas_state_variables, ncid, model_analysis_filename, &
nfields, variable_table )
TimeDimID = FindTimeDimension( ncid )
if (TimeDimID < 0 ) then
write(string1,*)'unable to find a dimension named Time.'
call error_handler(E_MSG,'static_init_model', string1, source, revision, revdate)
endif
call nc_check(nf90_Inquire(ncid,nDimensions,nVariables,nAttributes,unlimitedDimID), &
'static_init_model', 'inquire '//trim(model_analysis_filename))
if ( (TimeDimID > 0) .and. (unlimitedDimID > 0) .and. (TimeDimID /= unlimitedDimID)) then
write(string1,*)'IF Time is not the unlimited dimension, I am lost.'
call error_handler(E_MSG,'static_init_model', string1, source, revision, revdate)
endif
index1 = 1;
indexN = 0;
do ivar = 1, nfields
varname = trim(variable_table(ivar,1))
kind_string = trim(variable_table(ivar,2))
progvar(ivar)%varname = varname
progvar(ivar)%kind_string = kind_string
progvar(ivar)%dart_kind = get_index_for_quantity( progvar(ivar)%kind_string )
progvar(ivar)%numdims = 0
progvar(ivar)%numvertical = 1
progvar(ivar)%dimlens = MISSING_I
progvar(ivar)%numcells = MISSING_I
progvar(ivar)%numedges = MISSING_I
string2 = trim(model_analysis_filename)//' '//trim(varname)
call nc_check(nf90_inq_varid(ncid, trim(varname), VarID), &
'static_init_model', 'inq_varid '//trim(string2))
call nc_check(nf90_inquire_variable(ncid, VarID, xtype=progvar(ivar)%xtype, &
dimids=dimIDs, ndims=numdims), 'static_init_model', 'inquire '//trim(string2))
! If the long_name and/or units attributes are set, get them.
! They are not REQUIRED to exist but are nice to use if they are present.
if( nf90_inquire_attribute( ncid, VarID, 'long_name') == NF90_NOERR ) then
call nc_check( nf90_get_att(ncid, VarID, 'long_name' , progvar(ivar)%long_name), &
'static_init_model', 'get_att long_name '//trim(string2))
else
progvar(ivar)%long_name = varname
endif
if( nf90_inquire_attribute( ncid, VarID, 'units') == NF90_NOERR ) then
call nc_check( nf90_get_att(ncid, VarID, 'units' , progvar(ivar)%units), &
'static_init_model', 'get_att units '//trim(string2))
else
progvar(ivar)%units = 'unknown'
endif
! Since we are not concerned with the TIME dimension, we need to skip it.
! When the variables are read, only a single timestep is ingested into
! the DART state vector.
varsize = 1
dimlen = 1
DimensionLoop : do i = 1,numdims
if (dimIDs(i) == TimeDimID) cycle DimensionLoop
write(string1,'(''inquire dimension'',i2,A)') i,trim(string2)
call nc_check(nf90_inquire_dimension(ncid, dimIDs(i), len=dimlen, name=dimname), &
'static_init_model', string1)
progvar(ivar)%numdims = progvar(ivar)%numdims + 1
progvar(ivar)%dimlens(i) = dimlen
progvar(ivar)%dimname(i) = trim(dimname)
varsize = varsize * dimlen
select case ( dimname(1:6) )
case ('nCells')
progvar(ivar)%numcells = dimlen
case ('nEdges')
progvar(ivar)%numedges = dimlen
case ('nVertL') ! nVertLevels, nVertLevelsP1, nVertLevelsP2
progvar(ivar)%numvertical = dimlen
end select
enddo DimensionLoop
! this call sets: clamping, bounds, and out_of_range_fail in the progvar entry
call get_variable_bounds(mpas_state_bounds, ivar)
if (progvar(ivar)%numvertical == nVertLevels) then
progvar(ivar)%ZonHalf = .TRUE.
else
progvar(ivar)%ZonHalf = .FALSE.
endif
if (varname == 'u') has_edge_u = .true.
if (varname == 'uReconstructZonal' .or. &
varname == 'uReconstructMeridional') has_uvreconstruct = .true.
progvar(ivar)%varsize = varsize
progvar(ivar)%index1 = index1
progvar(ivar)%indexN = index1 + varsize - 1
index1 = index1 + varsize ! sets up for next variable
if ( debug > 4 ) call dump_progvar(ivar)
enddo
call nc_check( nf90_close(ncid), &
'static_init_model', 'close '//trim(model_analysis_filename))
model_size = progvar(nfields)%indexN
if ( debug > 0 .and. do_output()) then
write(logfileunit,*)
write( * ,*)
write(logfileunit,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i6))') &
nCells, nEdges, nVertices, nVertLevels
write( * ,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i6))') &
nCells, nEdges, nVertices, nVertLevels
write(logfileunit, *)'static_init_model: model_size = ', model_size
write( * , *)'static_init_model: model_size = ', model_size
if ( global_grid ) then
write(logfileunit, *)'static_init_model: grid is a global grid '
write( * , *)'static_init_model: grid is a global grid '
else
write(logfileunit, *)'static_init_model: grid has boundaries '
write( * , *)'static_init_model: grid has boundaries '
endif
if ( all_levels_exist_everywhere ) then
write(logfileunit, *)'static_init_model: all cells have same number of vertical levels '
write( * , *)'static_init_model: all cells have same number of vertical levels '
else
write(logfileunit, *)'static_init_model: cells have varying number of vertical levels '
write( * , *)'static_init_model: cells have varying number of vertical levels '
endif
endif
! do some sanity checking here:
! if you have at least one of these wind components in the state vector,
! you have to have them both. the subroutine will error out if only one
! is found and not both.
if (has_uvreconstruct) then
call winds_present(z1,m1,both)
endif
! if you set the namelist to use the reconstructed cell center winds,
! they have to be in the state vector.
if (update_u_from_reconstruct .and. .not. has_uvreconstruct) then
write(string1,*) 'update_u_from_reconstruct cannot be True'
write(string2,*) 'because state vector does not contain U/V reconstructed winds'
call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, &
text2=string2)
endif
! if you set the namelist to update the edge normal winds based on the
! updated cell center winds, and you also have the edge normal winds in
! the state vector, warn that the edge normal winds will be ignored
! when going back to the mpas_analysis.nc file. not an error, but the
! updates to the edge normal vectors won't affect the results.
if (update_u_from_reconstruct .and. has_edge_u) then
write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on U/V reconstructed winds'
write(string2,*) 'and not from the updated edge normal values in the state vector'
write(string3,*) 'because update_u_from_reconstruct is True'
call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, &
text2=string2, text3=string3)
endif
! there are two choices when updating the edge normal winds based on the
! winds at the cell centers. one is a direct interpolation of the values;
! the other is to read in the original cell centers, compute the increments
! changed by the interpolation, and then add or substract only the increments
! from the original edge normal wind values.
if (update_u_from_reconstruct) then
if (use_increments_for_u_update) then
write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on the difference'
write(string2,*) 'between the original U/V reconstructed winds and the updated values'
write(string3,*) 'because use_increment_for_u_update is True'
else
write(string1,*) 'edge normal winds (u) in MPAS file will be updated by averaging the'
write(string2,*) 'updated U/V reconstructed winds at the corresponding cell centers'
write(string3,*) 'because use_increment_for_u_update is False'
endif
call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, &
text2=string2, text3=string3)
endif
! basically we cannot do much without having at least these
! three fields in the state vector. refuse to go further
! if these are not present:
!
! TJH if ((get_progvar_index_from_kind(QTY_POTENTIAL_TEMPERATURE) < 0) .or. &
! TJH (get_progvar_index_from_kind(QTY_DENSITY) < 0) .or. &
! TJH (get_progvar_index_from_kind(QTY_VAPOR_MIXING_RATIO) < 0)) then
! TJH write(string1, *) 'State vector is missing one or more of the following fields:'
! TJH write(string2, *) 'Potential Temperature (theta), Density (rho), Vapor Mixing Ratio (qv).'
! TJH write(string3, *) 'Cannot convert between height/pressure nor compute sensible temps.'
! TJH call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, &
! TJH text2=string2, text3=string3)
! TJH endif
string1 = 'WARNING: fix block of required variables - detritus from atmosphere'
call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate)
allocate( ens_mean(model_size) )
end subroutine static_init_model
!------------------------------------------------------------------
subroutine get_state_meta_data(index_in, location, var_type)
! given an index into the state vector, return its location and
! if given, the var kind. despite the name, var_type is a generic
! kind, like those in obs_kind/obs_kind_mod.f90, starting with QTY_
! passed variables
integer, intent(in) :: index_in
type(location_type), intent(out) :: location
integer, optional, intent(out) :: var_type
! Local variables
integer :: nxp, nzp, iloc, vloc, nf, n
integer :: myindx
integer :: istatus
real(r8) :: depth
type(location_type) :: new_location
if ( .not. module_initialized ) call static_init_model
myindx = -1
nf = -1
! Determine the right variable
FindIndex : do n = 1,nfields
if( (progvar(n)%index1 <= index_in) .and. (index_in <= progvar(n)%indexN) ) THEN
nf = n
myindx = index_in - progvar(n)%index1 + 1
exit FindIndex
endif
enddo FindIndex
if( myindx == -1 ) then
write(string1,*) 'Problem, cannot find base_offset, index_in is: ', index_in
call error_handler(E_ERR,'get_state_meta_data',string1,source,revision,revdate)
endif
! Now that we know the variable, find the cell or edge
if ( progvar(nf)%numcells /= MISSING_I) then
nxp = progvar(nf)%numcells
elseif ( progvar(nf)%numedges /= MISSING_I) then
nxp = progvar(nf)%numedges
else
write(string1,*) 'ERROR, ',trim(progvar(nf)%varname),' is not defined on edges or cells'
call error_handler(E_ERR,'get_state_meta_data',string1,source,revision,revdate)
endif
nzp = progvar(nf)%numvertical
iloc = 1 + (myindx-1) / nzp ! cell index
vloc = myindx - (iloc-1)*nzp ! vertical level index
! the zGrid array contains the location of the cell top and bottom faces, so it has one
! more value than the number of cells in each column. for locations of cell centers
! you have to take the midpoint of the top and bottom face of the cell.
if (progvar(nf)%numedges /= MISSING_I) then
if (.not. data_on_edges) then
call error_handler(E_ERR, 'get_state_meta_data', &
'Internal error: numedges present but data_on_edges false', &
source, revision, revdate, text2='variable '//trim(progvar(nf)%varname))
endif
if ( progvar(nf)%ZonHalf ) then
depth = zGridEdge(vloc,iloc)
else
call error_handler(E_ERR, 'get_state_meta_data', 'no support for edges at face depths', &
source, revision, revdate)
endif
else
if ( progvar(nf)%ZonHalf ) then
depth = zGridCenter(vloc,iloc)
else if (nzp <= 1) then
depth = zGridFace(1,iloc)
else
depth = zGridFace(vloc,iloc)
endif
endif
if (progvar(nf)%numedges /= MISSING_I) then
if (.not. data_on_edges) then
call error_handler(E_ERR, 'get_state_meta_data', &
'Internal error: numedges present but data_on_edges false', &
source, revision, revdate, text2='variable '//trim(progvar(nf)%varname))
endif
if (nzp <= 1) then
location = set_location(lonEdge(iloc),latEdge(iloc), depth, VERTISSURFACE)
else
location = set_location(lonEdge(iloc),latEdge(iloc), depth, VERTISHEIGHT)
endif
else ! must be on cell centers
if (nzp <= 1) then
location = set_location(lonCell(iloc),latCell(iloc), depth, VERTISSURFACE)
else
location = set_location(lonCell(iloc),latCell(iloc), depth, VERTISHEIGHT)
endif
endif
! Let us return the vert location with the requested vertical localization coordinate
! hoping that the code can run faster when same points are close to many obs
! since vert_convert does not need to be called repeatedly in get_close_obs any more.
! FIXME: we should test this to see if it's a win (in computation time). for obs
! which are not dense relative to the grid, this might be slower than doing the
! conversions on demand in the localization code (in get_close_obs()).
if ( .not. horiz_dist_only .and. vert_localization_coord /= VERTISHEIGHT ) then
new_location = location
call vert_convert(ens_mean, new_location, progvar(nf)%dart_kind, vert_localization_coord, istatus)
if(istatus == 0) location = new_location
endif
if (debug > 12) then
write(*,'("INDEX_IN / myindx / IVAR / NX, NZ: ",2(i10,2x),3(i5,2x))') index_in, myindx, nf, nxp, nzp
write(*,'(" ILOC, KLOC: ",2(i5,2x))') iloc, vloc
write(*,'(" LON/LAT/HGT: ",3(f12.3,2x))') lonCell(iloc), latCell(iloc), depth
endif
if (present(var_type)) then
var_type = progvar(nf)%dart_kind
endif
end subroutine get_state_meta_data
!------------------------------------------------------------------
subroutine model_interpolate(x, location, obs_type, interp_val, istatus)
! given a state vector, a location, and a QTY_xxx, return the
! interpolated value at that location, and an error code. 0 is success,
! anything positive is an error. (negative reserved for system use)
!
! ERROR codes:
!
! ISTATUS = 99: general error in case something terrible goes wrong...
! ISTATUS = 88: this kind is not in the state vector
! ISTATUS = 11: Could not find a triangle that contains this lat/lon
! ISTATUS = 12: depth vertical coordinate out of model range.
! ISTATUS = 13: Missing value in interpolation.
! ISTATUS = 16: Don't know how to do vertical velocity for now
! ISTATUS = 17: Unable to compute pressure values
! ISTATUS = 18: altitude illegal
! ISTATUS = 19: could not compute u using RBF code
! ISTATUS = 101: Internal error; reached end of subroutine without
! finding an applicable case.
!
! passed variables
real(r8), intent(in) :: x(:)
type(location_type), intent(in) :: location
integer, intent(in) :: obs_type
real(r8), intent(out) :: interp_val
integer, intent(out) :: istatus
! local storage
integer :: ivar, obs_kind
integer :: tvars(3)
logical :: goodkind
real(r8) :: values(3), lpres
if ( .not. module_initialized ) call static_init_model
interp_val = MISSING_R8
istatus = 99 ! must be positive (and integer)
! rename for sanity - we can't change the argument names
! to this subroutine, but this really is a kind.
obs_kind = obs_type
if (debug > 0) then
call write_location(0,location,charstring=string1)
print *, my_task_id(), 'stt x(1), kind, loc ', x(1), obs_kind, trim(string1)
endif
! if we can interpolate any other kinds that are not directly in the
! state vector, then add those cases here.
ivar = get_progvar_index_from_kind(obs_kind)
if (ivar > 0) then
goodkind = .true.
else
goodkind = .false.
! exceptions if the kind isn't directly
! a field in the state vector:
select case (obs_kind)
!case (QTY_TEMPERATURE) ! potential temperature is in state vector not in-situ temp
! goodkind = .true.
!case (QTY_PRESSURE)
! goodkind = .true.
end select
endif
! this kind is not in the state vector and it isn't one of the exceptions
! that we know how to handle.
if (.not. goodkind) then
if (debug > 4) print *, 'kind rejected', obs_kind
istatus = 88
goto 100
endif
! if you add exceptions above, then you need code like this:
! if (obs_kind == QTY_xxx) then
! add code here for how to interpolate it.
!
! compute_scalar_with_barycentric() can take up to 3 kinds at one location
! and return up to 3 values which can be combined or computed on here.
!
! to use the RBF functions, here is an example call. the third arg
! controls whether to return the meridional (.false.) or zonal (.true.) component.
! call compute_u_with_rbf(x, location, .TRUE., interp_val, istatus)
!
! on error, goto 100 to get the final error checking code
!
! else
! direct interpolation, kind is in the state vector
tvars(1) = ivar
call compute_scalar_with_barycentric(x, location, 1, tvars, values, istatus)
interp_val = values(1)
if (debug > 4) print *, 'called generic compute_w_bary, kind, val, istatus: ', obs_kind, interp_val, istatus
if (istatus /= 0) goto 100