-
Notifications
You must be signed in to change notification settings - Fork 145
/
model_mod.f90
7554 lines (5801 loc) · 276 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
module model_mod
! MPAS Atmosphere 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.
! This revision of the model_mod supports both a global MPAS grid and
! a regional grid. For the regional grid only observations which
! are completely inside the interior will be assimilated, meaning obs
! which need interpolation information from the boundary cells
! (in any of the 7 boundary layers) will be rejected. However, during the
! assimilation phase all locations in the local grid will be impacted,
! even locations in the boundary layers if there are obs close to the
! boundaries. A post-processing step will smooth the GFS external
! values with the values updated by the assimilation in the boundary layers.
! Note that to reject obs during interpolation requires the model_interpolate()
! routine to check and return an error, but during the vertical conversion and
! get_close routines all state point operations must succeed, even those
! in the boundary layers. Pay close attention to which internal routines are used
! by each to make sure the intended actions are what happens.
use types_mod, only : r4, r8, i8, 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_type, set_location, get_location, &
write_location, vertical_localization_on, &
VERTISUNDEF, VERTISSURFACE, VERTISLEVEL, &
VERTISPRESSURE, VERTISHEIGHT, VERTISSCALEHEIGHT, &
loc_get_close_obs => get_close_obs, &
loc_get_close_state => get_close_state, &
is_vertical, set_vertical_localization_coord
use netcdf_utilities_mod, only : nc_add_global_attribute, nc_synchronize_file, &
nc_add_global_creation_time, nc_check, &
nc_begin_define_mode, nc_end_define_mode, &
nc_open_file_readonly, nc_close_file, &
nc_add_attribute_to_variable, nc_define_dimension, &
nc_define_unlimited_dimension, nc_define_character_variable, &
nc_define_real_variable, nc_get_variable, nc_put_variable, &
nc_get_dimension_size, nc_variable_exists, nc_dimension_exists, &
nc_define_integer_variable
use location_io_mod, only : nc_write_location_atts, nc_write_location
use default_model_mod, only : nc_write_model_vars, adv_1step, &
init_time => fail_init_time, &
init_conditions => fail_init_conditions
use xyz_location_mod, only : xyz_location_type, xyz_set_location, xyz_get_location, &
xyz_get_close_type, xyz_get_close_init, xyz_get_close_destroy, &
xyz_find_nearest, xyz_use_great_circle_dist
use utilities_mod, only : register_module, error_handler, get_unit, &
E_ERR, E_WARN, E_MSG, E_ALLMSG, logfileunit, &
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, scalar
use obs_kind_mod, only : get_index_for_quantity, &
get_name_for_quantity, &
get_quantity_for_type_of_obs, &
QTY_SURFACE_ELEVATION, &
QTY_SURFACE_PRESSURE, &
QTY_10M_U_WIND_COMPONENT, &
QTY_10M_V_WIND_COMPONENT, &
QTY_2M_TEMPERATURE, &
QTY_2M_SPECIFIC_HUMIDITY, &
QTY_VERTICAL_VELOCITY, &
QTY_POTENTIAL_TEMPERATURE, &
QTY_EDGE_NORMAL_SPEED, &
QTY_TEMPERATURE, &
QTY_U_WIND_COMPONENT, &
QTY_V_WIND_COMPONENT, &
QTY_PRESSURE, &
QTY_DENSITY, &
QTY_VAPOR_MIXING_RATIO, &
QTY_CLOUDWATER_MIXING_RATIO, &
QTY_RAINWATER_MIXING_RATIO, &
QTY_ICE_MIXING_RATIO, &
QTY_SNOW_MIXING_RATIO, &
QTY_GRAUPEL_MIXING_RATIO, &
QTY_SPECIFIC_HUMIDITY, &
QTY_GEOPOTENTIAL_HEIGHT, &
QTY_PRECIPITABLE_WATER, &
QTY_SKIN_TEMPERATURE, & ! for rttov
QTY_SURFACE_TYPE, & ! for rttov
QTY_CLOUD_FRACTION ! for rttov
use mpi_utilities_mod, only: my_task_id, all_reduce_min_max
use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars
use distributed_state_mod
! 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
use state_structure_mod, only : add_domain, get_model_variable_indices, &
state_structure_info, get_index_start, get_index_end, &
get_num_variables, get_domain_size, get_varid_from_varname, &
get_variable_name, get_num_dims, get_dim_lengths, &
get_dart_vector_index, get_num_varids_from_kind
implicit none
private
! these routines must be public and you cannot change
! the arguments - they will be called *from* the DART code.
! routines in this list have code in this module
public :: get_model_size, &
get_num_vars, &
get_state_meta_data, &
model_interpolate, &
shortest_time_between_assimilations, &
static_init_model, &
end_model, &
nc_write_model_atts, &
pert_model_copies, &
get_close_obs, &
get_close_state, &
convert_vertical_obs, &
convert_vertical_state, &
read_model_time, &
write_model_time
! code for these routines are in other modules
public :: init_time, &
init_conditions, &
adv_1step, &
nc_write_model_vars
! generally useful routines for various support purposes.
! the interfaces here can be changed as appropriate.
public :: get_init_template_filename, &
analysis_file_to_statevector, &
statevector_to_analysis_file, &
statevector_to_boundary_file, &
get_analysis_time, &
get_grid_dims, &
get_xland, &
get_surftype, &
get_cell_center_coords, &
get_bdy_mask, &
print_variable_ranges, &
find_closest_cell_center, &
find_triangle, &
read_2d_from_nc_file, &
find_height_bounds, &
cell_ok_to_interpolate, &
is_global_grid, &
uv_cell_to_edges
! try adjusting what static_init_model does, before it is called.
! the main update_bc program, for example, can call these before
! calling static_init_model() and that will modify what it does.
public :: set_lbc_variables, &
force_u_into_state
! set_lbc_variables sets the lbc_variables string array
! force_u_into_state sets a logical add_u_to_state_list that forces u to be in state
! version controlled file description for error handling, do not edit
character(len=*), parameter :: source = 'models/mpas_atm/model_mod.f90'
character(len=*), parameter :: revision = ''
character(len=*), parameter :: revdate = ''
! module global storage; maintains values between calls, accessible by
! any subroutine
character(len=512) :: string1, string2, string3
character(len=256) :: locstring
logical, save :: module_initialized = .false.
! length of an mpas (also wrf) time string: YYYY-MM-DD_hh:mm:ss
integer, parameter :: TIMELEN = 19
! Real (physical) constants as defined exactly in MPAS.
! redefined here for consistency with the model.
real(r8), parameter :: rgas = 287.0_r8 ! Constant: Gas constant for dry air [J kg-1 K-1]
real(r8), parameter :: rv = 461.6_r8 ! Constant: Gas constant for water vapor [J kg-1 K-1]
real(r8), parameter :: cp = 7.0_r8*rgas/2.0_r8 ! = 1004.5 Constant: Specific heat of dry air at constant pressure [J kg-1 K-1]
real(r8), parameter :: cv = cp - rgas ! = 717.5 Constant: Specific heat of dry air at constant volume [J kg-1 K-1]
real(r8), parameter :: p0 = 100000.0_r8
real(r8), parameter :: rcv = rgas/(cp-rgas)
real(r8), parameter :: rvord = rv/rgas
! earth radius; needed to convert lat/lon to x,y,z cartesian coords.
! for the highest accuracy this should match what the model uses.
real(r8), parameter :: radius = 6371229.0_r8 ! meters
! roundoff error tolerance
! set in static_init_model() to 1e-5 or 1e-12 depending on compiled precision
real(r8) :: roundoff
! 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(:)
logical :: search_initialized = .false.
! compile-time control over whether grid information is written to the
! diagnostic files or not. if it is, the files are self-contained (e.g. for
! ease of plotting), but are also much larger than they would be otherwise.
! change this private variable to control whether the grid information is
! written or not.
logical :: add_static_data_to_diags = .false.
! variables which are in the module namelist
character(len=256) :: init_template_filename = 'mpas_init.nc'
character(len=256) :: bdy_template_filename = ''
integer :: vert_localization_coord = VERTISHEIGHT
integer :: assimilation_period_days = 0
integer :: assimilation_period_seconds = 21600
real(r8) :: model_perturbation_amplitude = 0.0001 ! tiny amounts
logical :: log_p_vert_interp = .true. ! if true, interpolate vertical pressure in log space
character(len=32) :: calendar = 'Gregorian'
real(r8) :: highest_obs_pressure_mb = 100.0_r8 ! do not assimilate obs higher than this level.
real(r8) :: sfc_elev_max_diff = -1.0_r8 ! do not assimilate if |model - station| height is larger than this [m].
integer :: xyzdebug = 0
integer :: debug = 0 ! turn up for more and more debug messages
! this is not in the namelist or supported generally.
! (setting this to true avoids the surface elevation max diff
! test for elevation and surface pressure.)
logical :: always_assim_surf_altimeters = .false.
integer :: anl_domid = -1 ! For state_structure_mod access
integer :: lbc_domid = -1
! 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.
! if set by: call force_u_into_state() *BEFORE* calling static_init_model(),
! the code will add U (edge normal winds) to the mpas state vector even if it
! isn't in the namelist.
logical :: add_u_to_state_list = .false.
! if > 0, amount of distance in fractional model levels that a vertical
! point can be above or below the top or bottom of the grid and still be
! evaluated without error. if extrapolate is true, extrapolate from the
! first or last model level. if extrapolate is false, use level 1 or N.
real(r8) :: outside_grid_level_tolerance = -1.0_r8
logical :: extrapolate = .false.
! if the calling code updates an existing file it simply writes the state variable
! data. if it needs to create a file from scratch it calls nc_write_model_atts()
! to optionally add grid info or any other non-state variables or attributes.
! setting this to .true. adds the mpas grid info to the file. .false. does
! not and results in smaller diag/output files.
logical :: write_grid_to_diag_files = .false.
! in converting to scale height for the vertical, set this to .false. to
! use simply the log of the pressure. to normalize by the surface pressure
! (backwards compatible with previous code), set this to .true.
logical :: no_normalization_of_scale_heights = .true.
! for regional MPAS
real(r8) :: dxmax ! max distance between two adjacent cell centers in the mesh (in meters)
! when updating boundary files for regional mpas, note whether the boundary
! file has the reconstructed winds (lbc_ur, lbc_vr) or not. (this is set by
! looking at the bdy template file and seeing if those variables are there.)
! if not, the other two options are ignored.
! if they are in the lbc file, then the other logicals control whether to use
! them instead of updating the U edge winds directly, and whether to use the
! reconstructed increments or not.
! the latter two options could be added to the namelist if someone wanted to
! explore the options for how the edge winds are updated in the boundary file.
! for now they're not - they're hardcoded true.
! note that these are for the boundary file update only - there are separate
! options for how to update the U winds in the main assimilation state.
logical :: lbc_file_has_reconstructed_winds = .false.
namelist /model_nml/ &
init_template_filename, &
vert_localization_coord, &
assimilation_period_days, &
assimilation_period_seconds, &
model_perturbation_amplitude, &
log_p_vert_interp, &
calendar, &
debug, &
xyzdebug, &
use_u_for_wind, &
use_rbf_option, &
update_u_from_reconstruct, &
use_increments_for_u_update, &
highest_obs_pressure_mb, &
outside_grid_level_tolerance, &
extrapolate, &
sfc_elev_max_diff, &
write_grid_to_diag_files, &
no_normalization_of_scale_heights
! 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 )
! this is special and not in the namelist. boundary files have a fixed
! set of variables with fixed names.
character(len=NF90_MAX_NAME) :: lbc_variables(max_state_variables) = ''
namelist /mpas_vars_nml/ mpas_state_variables, mpas_state_bounds
integer :: nfields
!>@todo FIXME - REMOVE AS MUCH OF THIS AS POSSIBLE.
!> some of this information is in the state structure now.
!> the duplicate progvar stuff should be removed and the
!> state routines used instead. this duplicates work and
!> makes us keep up code in 2 different places.
! original code:
! Everything needed to describe a variable
type progvartype
private
character(len=NF90_MAX_NAME) :: varname
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
integer :: nSoilLevels = -1 ! Number of soil layers
! 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 :: dcEdge(:) ! distance between two adjacent cell centers (in meters)
real(r8), allocatable :: xland(:) ! land-ocean mask (1=land including sea-ice ; 2=ocean)
real(r8), allocatable :: seaice(:) ! sea-ice flag (0=no seaice; =1 otherwise) - for rttov
real(r8), allocatable :: skintemp(:)! ground or water surface temperature - for rttov
real(r8), allocatable :: zGridFace(:,:) ! geometric height at cell faces (nVertLevelsP1,nCells)
real(r8), allocatable :: zGridCenter(:,:) ! geometric height at cell centers (nVertLevels, nCells)
real(r8), allocatable :: zGridEdge(:,:) ! geometric height at edge centers (nVertLevels, nEdges)
!real(r8), allocatable :: zEdgeFace(:,:) ! geometric height at edges faces (nVertLevelsP1,nEdges)
!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 :: bdyMaskCell(:)
integer, allocatable :: bdyMaskEdge(:)
integer, allocatable :: maxLevelCell(:)
integer, allocatable :: surftype(:) ! ! surface type (land=0, water=1, seaice = 2) - for rttov
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
interface write_model_time
module procedure write_model_time_file
module procedure write_model_time_restart
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()
!>@todo FIXME - can we add an optional arg here for update_bc use?
! 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, lbc_nfields
logical :: both
real(r8) :: variable_bounds(max_state_variables, 2)
integer :: variable_qtys(max_state_variables)
character(len=*), parameter :: routine = 'static_init_model'
if ( module_initialized ) return ! only need to do this once.
! Since this routine calls other routines that could call this routine
! we'll say we've been initialized pretty dang early.
module_initialized = .true.
! Print module information to log file and stdout.
call register_module(source, revision, revdate)
! 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 window is ',dd,' days ',ss,' seconds'
call error_handler(E_MSG,routine,string1,source,revision,revdate)
!---------------------------------------------------------------
! 1) get grid dimensions
! 2) allocate space for the grids
! 3) read them from the analysis file
ncid = nc_open_file_readonly(init_template_filename, routine)
! move this up - some of the routines below depend on having
! the variable table parsed already, and U added if needed.
call verify_state_variables( mpas_state_variables, ncid, init_template_filename, &
nfields, variable_table)
! get sizes
call read_grid_dims(ncid)
allocate(latCell(nCells), lonCell(nCells))
allocate(zGridFace(nVertLevelsP1, nCells))
allocate(zGridCenter(nVertLevels, nCells))
allocate(cellsOnVertex(vertexDegree, nVertices))
allocate(dcEdge(nEdges))
allocate(nEdgesOnCell(nCells))
allocate(xland(nCells))
allocate(seaice(nCells)) ! for rttov
allocate(skintemp(nCells)) ! for rttov
allocate(surftype(nCells)) ! for rttov
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(variable_table, lbc_variables)
if(data_on_edges) then
allocate(zGridEdge(nVertLevels, nEdges))
allocate(xEdge(nEdges), yEdge(nEdges), zEdge(nEdges))
allocate(latEdge(nEdges), lonEdge(nEdges))
endif
! is this a global or regional grid?
! if regional, allocate and read in the boundary info here.
call set_global_grid(ncid)
! fill in the grid values
call get_grid(ncid)
! vertical faces are in the input file. compute vertical center locations here.
do kloc=1, nCells
do iloc=1, nVertLevels
zGridCenter(iloc,kloc) = (zGridFace(iloc,kloc) + zGridFace(iloc+1,kloc))*0.5_r8
enddo
enddo
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 height 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.
!
! THIS CODE SHOULD BE REMOVED - it is done by the add_domain code.
!
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(init_template_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(init_template_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))
! 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
case ('nSoilL') ! nSoilLevels
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 > 11 .and. do_output()) call dump_progvar(ivar)
enddo
call nc_close_file(ncid, routine)
! FIXME: moved below to after add_domain() calls
!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,i9))') &
nCells, nEdges, nVertices, nVertLevels
write( * ,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i9))') &
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 is NOT a global grid. Lateral boundaries exist '
write( * , *)'static_init_model: grid is NOT a global grid. Lateral boundaries exist '
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
! set the domain(s) in the state structure here
variable_bounds(1:nfields, 1) = progvar(1:nfields)%range(1)
variable_bounds(1:nfields, 2) = progvar(1:nfields)%range(2)
variable_qtys(1:nfields) = progvar(1:nfields)%dart_kind
anl_domid = add_domain(init_template_filename, nfields, &
var_names = variable_table (1:nfields,1), &
kind_list = variable_qtys(1:nfields), &
clamp_vals = variable_bounds(1:nfields,:) )
model_size = get_domain_size(anl_domid)
if ( debug > 4 .and. do_output()) print*,'model_size(anl_domid)=',model_size ! HA
lbc_nfields = 0
! if we have a lateral boundary file, add it to the domain
! so we have access to the corresponding lbc_xxx fields.
!>@todo FIXME: if we want to do increments, we could also add a
! third domain which is the original forecast fields before
! the assimilation (so we can compute increments)
if (.not. global_grid .and. lbc_variables(1) /= '') then
! regional: count number of lbc fields to read in
COUNTUP: do i=1, max_state_variables
if (lbc_variables(i) /= '') then
lbc_nfields = lbc_nfields + 1
else
exit COUNTUP
endif
enddo COUNTUP
if( debug > 4 .and. do_output()) print*, 'Regional: number of lbc fields to read in = ', lbc_nfields
lbc_domid = add_domain(bdy_template_filename, lbc_nfields, &
var_names = lbc_variables)
! FIXME clamp_vals = variable_bounds(1:nfields,:) )
if( debug > 4 .and. do_output()) print*, 'model_size, lbc_domid =',model_size,lbc_domid
model_size = model_size + get_domain_size(lbc_domid)
else
lbc_domid = -1
endif
if ( debug > 4 .and. do_output()) then
call state_structure_info(anl_domid)
print*, 'model_size =',model_size
if(lbc_domid >= 0) print*, 'get_domain_size(lbc_domid):',get_domain_size(lbc_domid)
endif
!if ( debug > 4 .and. do_output() .and. lbc_domid >= 0) call state_structure_info(lbc_domid)
! 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:
!print *, get_num_varids_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE)
!print *, get_num_varids_from_kind(anl_domid, QTY_DENSITY)
!print *, get_num_varids_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO)
if ((get_num_varids_from_kind(anl_domid, QTY_POTENTIAL_TEMPERATURE) < 0) .or. &
(get_num_varids_from_kind(anl_domid, QTY_DENSITY) < 0) .or. &
(get_num_varids_from_kind(anl_domid, QTY_VAPOR_MIXING_RATIO) < 0)) then
write(string1, *) 'State vector is missing one or more of the following fields:'
write(string2, *) 'Potential Temperature (theta), Density (rho), Vapor Mixing Ratio (qv).'
write(string3, *) 'Cannot convert between height/pressure nor compute sensible temps.'
call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, &
text2=string2, text3=string3)
endif
if (extrapolate) then
call error_handler(E_MSG,'static_init_model',&
'extrapolate not supported yet; will use level 1 or N values', &
source, revision, revdate)
endif
! tell the location module how we want to localize in the vertical
call set_vertical_localization_coord(vert_localization_coord)
! set an appropriate value for roundoff tests based
! on this code being compiled single or double precision.
! set to 1e-5 (for single) or 1e-12 (for double precision).
if (r8 == digits12) then
roundoff = 1.0e-12_r8
else
roundoff = 1.0e-5_r8
endif
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(i8), intent(in) :: index_in
type(location_type), intent(out) :: location
integer, optional, intent(out) :: var_type
! Local variables
integer :: nzp, iloc, vloc, nf, ndim
real(r8) :: height
type(location_type) :: new_location
if ( .not. module_initialized ) call static_init_model
! get the local indicies and type from dart index. kloc is a dummy variable for this subroutine
call find_mpas_indices(index_in, iloc, vloc, ndim, nf)
nzp = progvar(nf)%numvertical
! 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
height = zGridEdge(vloc,iloc)
else
call error_handler(E_ERR, 'get_state_meta_data', 'no support for edges at face heights', &