-
Notifications
You must be signed in to change notification settings - Fork 145
/
Copy pathmodel_mod.f90
8428 lines (6617 loc) · 342 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
!> @mainpage Remote Memory Access version of DART code.
!>
!> Forward operator, vertical conversion
!>
!> WRF
!>
!> For distributed phb array see distributed_phb_model_mod.f90
!> \todo To do list
!> @author dart@ucar.edu
!>
!> \subpage test
!>
!> \page test
!> WRF model mod
module model_mod
! Assimilation interface for WRF model
!> \defgroup wrf model_mod
!> Model mod
!>
!> Model mod for WRF
!> @{
!-----------------------------------------------------------------------
!
! interface for WRF
!
!-----------------------------------------------------------------------
!---------------- m o d u l e i n f o r m a t i o n ------------------
!-----------------------------------------------------------------------
use types_mod, only : r8, i8, deg2rad, missing_r8, ps0, earth_radius, &
gas_constant, gas_constant_v, gravity, pi, &
digits12
use time_manager_mod, only : time_type, set_time, set_calendar_type, GREGORIAN,&
set_date, get_date
use location_mod, only : location_type, get_location, set_location, &
query_location, VERTISUNDEF, VERTISSURFACE, &
VERTISLEVEL, VERTISPRESSURE, VERTISHEIGHT, &
VERTISSCALEHEIGHT, vertical_localization_on, &
set_vertical_localization_coord, &
get_close_type, get_dist, is_vertical, &
loc_get_close => get_close_obs
use utilities_mod, only : file_exist, open_file, close_file, &
register_module, error_handler, E_ERR, E_WARN, &
E_MSG, nmlfileunit, logfileunit, do_output, &
find_namelist_in_file, check_namelist_read, &
find_textfile_dims, file_to_text, &
do_nml_file, do_nml_term, scalar
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
use mpi_utilities_mod, only : my_task_id, task_count, all_reduce_min_max
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
use obs_kind_mod, only : QTY_U_WIND_COMPONENT, QTY_V_WIND_COMPONENT, &
QTY_SURFACE_PRESSURE, QTY_TEMPERATURE, &
QTY_SPECIFIC_HUMIDITY, QTY_SURFACE_ELEVATION, &
QTY_PRESSURE, QTY_VERTICAL_VELOCITY, &
QTY_DENSITY, QTY_FLASH_RATE_2D, &
QTY_RAINWATER_MIXING_RATIO, QTY_HAIL_MIXING_RATIO, &
QTY_GRAUPEL_MIXING_RATIO, QTY_SNOW_MIXING_RATIO, &
QTY_CLOUDWATER_MIXING_RATIO, QTY_ICE_MIXING_RATIO, &
QTY_CLOUD_FRACTION, QTY_CONDENSATIONAL_HEATING, &
QTY_VAPOR_MIXING_RATIO, QTY_2M_TEMPERATURE, &
QTY_2M_SPECIFIC_HUMIDITY, QTY_10M_U_WIND_COMPONENT, &
QTY_10M_V_WIND_COMPONENT, &
QTY_ICE_NUMBER_CONCENTRATION, QTY_GEOPOTENTIAL_HEIGHT, &
QTY_POTENTIAL_TEMPERATURE, QTY_SOIL_MOISTURE, &
QTY_DROPLET_NUMBER_CONCENTR, QTY_SNOW_NUMBER_CONCENTR, &
QTY_RAIN_NUMBER_CONCENTR, QTY_GRAUPEL_NUMBER_CONCENTR, &
QTY_HAIL_NUMBER_CONCENTR, QTY_HAIL_VOLUME, &
QTY_GRAUPEL_VOLUME, QTY_DIFFERENTIAL_REFLECTIVITY, &
QTY_RADAR_REFLECTIVITY, QTY_POWER_WEIGHTED_FALL_SPEED, &
QTY_SPECIFIC_DIFFERENTIAL_PHASE, &
QTY_VORTEX_LAT, QTY_VORTEX_LON, &
QTY_VORTEX_PMIN, QTY_VORTEX_WMAX, &
QTY_SKIN_TEMPERATURE, QTY_LANDMASK, &
QTY_SURFACE_TYPE, &
get_index_for_quantity, get_num_quantities, &
get_name_for_quantity
use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars
use sort_mod, only : sort
use distributed_state_mod, only : get_state
use default_model_mod, only : adv_1step, nc_write_model_vars, &
init_conditions => fail_init_conditions, &
init_time => fail_init_time
use state_structure_mod, only : add_domain, get_model_variable_indices, &
state_structure_info, &
get_index_start, get_index_end, &
get_dart_vector_index
! FIXME:
! the kinds
! QTY_ICE_NUMBER_CONCENTRATION should be QTY_ICE_NUMBER_CONCENTR
! to be consistent with the other concentration names.
!nc -- module_map_utils split the declarations of PROJ_* into a separate module called
!nc -- misc_definitions_module
use map_utils, only : proj_info, map_init, map_set, latlon_to_ij, &
ij_to_latlon, gridwind_to_truewind
use misc_definitions_module, only : PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS, PROJ_CASSINI, &
PROJ_CYL
use netcdf
use typesizes
implicit none
private
!-----
! DART requires 18 specific public interfaces from model_mod.f90
! routines with code in this module
public :: get_model_size, &
get_state_meta_data, &
shortest_time_between_assimilations, &
static_init_model, &
model_interpolate, &
nc_write_model_atts, &
get_close_obs, &
get_close_state, &
convert_vertical_obs, &
convert_vertical_state, &
pert_model_copies, &
read_model_time, &
write_model_time, &
end_model
! required routines where the code is in another module
public :: adv_1step, &
init_time, &
init_conditions, &
nc_write_model_vars
!-----
! Here is the appropriate place for other users to make additional routines
! contained within model_mod available for public use:
public :: get_number_domains, &
get_wrf_static_data, &
model_pressure_distrib, &
model_height_distrib, &
pres_to_zk, &
height_to_zk, &
get_domain_info, &
get_wrf_state_variables, &
fill_default_state_table, &
read_wrf_dimensions, &
get_number_of_wrf_variables, &
get_variable_bounds, &
set_variable_bound_defaults, &
get_variable_size_from_file, &
get_wrf_date, set_wrf_date, &
vert_convert, &
height_diff_check
! public parameters
public :: max_state_variables, &
num_state_table_columns, &
num_bounds_table_columns
! types
public :: wrf_dom, wrf_static_data_for_dart
! Interfaces for array and single value versions of subroutines/functions
! This is because the forward operator works on the whole ensemble, and the
! vertical conversion only uses the mean copy.
! HK ? interp_4pressure interface needed?
!-----------------------------------------------------------------------
! version controlled file description for error handling, do not edit
character(len=*), parameter :: source = 'wrf/model_mod.f90'
character(len=*), parameter :: revision = ''
character(len=*), parameter :: revdate = ''
! miscellaneous
integer, parameter :: max_state_variables = 100
integer, parameter :: num_state_table_columns = 5
integer, parameter :: num_bounds_table_columns = 4
!-----------------------------------------------------------------------
! Model namelist parameters with default values.
!
! center_search_half_length: half length (in meter) of the searching box to locate
! minimum pressure at a grid point
! center_spline_scale: coarse grid to spline interp. fine grid ratio
!-----------------------------------------------------------------------
logical :: output_state_vector = .false. ! output prognostic variables
logical :: default_state_variables = .true. ! use default state list?
character(len=129) :: wrf_state_variables(num_state_table_columns,max_state_variables) = 'NULL'
character(len=129) :: wrf_state_bounds(num_bounds_table_columns,max_state_variables) = 'NULL'
integer :: num_domains = 1
integer :: calendar_type = GREGORIAN
integer :: assimilation_period_seconds = 21600
! Max height a surface obs can be away from the actual model surface
! and still be accepted (in meters)
real (kind=r8) :: sfc_elev_max_diff = -1.0_r8 ! could be something like 200.0_r8
real (kind=r8) :: center_search_half_length = 500000.0_r8
real(r8) :: circulation_pres_level = 80000.0_r8
real(r8) :: circulation_radius = 108000.0_r8
integer :: center_spline_grid_scale = 10
integer :: vert_localization_coord = VERTISHEIGHT
! Allow observations above the surface but below the lowest sigma level.
logical :: allow_obs_below_vol = .false.
! Do the interpolation of pressure values only after taking the log (.true.)
! vs doing a linear interpolation directly in pressure units (.false.)
logical :: log_vert_interp = .true.
logical :: log_horz_interpM = .false.
logical :: log_horz_interpQ = .false.
!nc -- we are adding these to the model.nml until they appear in the NetCDF files
logical :: polar = .false. ! wrap over the poles
logical :: periodic_x = .false. ! wrap in longitude or x
logical :: periodic_y = .false. ! used for single column model, wrap in y
!JPH -- single column model flag
logical :: scm = .false. ! using the single column model
logical :: allow_perturbed_ics = .false. ! should spin the model up for a while after
! obsolete items; ignored by this code.
! non-backwards-compatible change. should be removed,
! but see note below about namelist.
integer :: num_moist_vars
logical :: surf_obs, soil_data, h_diab
! adv_mod_command moved to dart_to_wrf namelist; ignored here.
character(len = 72) :: adv_mod_command = ''
! num_moist_vars, surf_obs, soil_data, h_diab, and adv_mod_command
! are IGNORED no matter what their settings in the namelist are.
! they are obsolete, but removing them here will cause a fatal error
! until users remove them from their input.nml files as well.
namelist /model_nml/ num_moist_vars, &
num_domains, calendar_type, surf_obs, soil_data, h_diab, &
default_state_variables, wrf_state_variables, &
wrf_state_bounds, sfc_elev_max_diff, &
adv_mod_command, assimilation_period_seconds, &
allow_obs_below_vol, vert_localization_coord, &
center_search_half_length, center_spline_grid_scale, &
circulation_pres_level, circulation_radius, polar, &
periodic_x, periodic_y, scm, allow_perturbed_ics
! if you need to check backwards compatibility, set this to .true.
! otherwise, leave it as false to use the more correct geometric height
logical :: use_geopotential_height = .false.
character(len = 20) :: wrf_nml_file = 'namelist.input'
logical :: have_wrf_nml_file = .false.
integer :: num_obs_kinds = 0
logical, allocatable :: in_state_vector(:)
integer, allocatable :: domain_id(:) ! Global storage to interface with state_structure_mod.
!-----------------------------------------------------------------------
! Private definition of domain map projection use by WRF
!nc -- added in CASSINI and CYL according to module_map_utils convention
!JPH -- change variable name from map_sphere to more specific map_latlon
integer, parameter :: map_latlon = 0, map_lambert = 1, map_polar_stereo = 2, map_mercator = 3
integer, parameter :: map_cassini = 6, map_cyl = 5
! Private definition of model variable types
real (kind=r8), PARAMETER :: kappa = 2.0_r8/7.0_r8 ! gas_constant / cp
real (kind=r8), PARAMETER :: ts0 = 300.0_r8 ! Base potential temperature for all levels.
!---- private data ----
! Got rid of surf_var as global private variable for model_mod and just defined it locally
! within model_interpolate
TYPE wrf_static_data_for_dart
integer :: bt, bts, sn, sns, we, wes, sls
real(r8) :: dx, dy, dt, p_top
integer :: map_proj
real(r8) :: cen_lat,cen_lon
type(proj_info) :: proj
! Boundary conditions -- hopefully one day these will be in the global attributes of the
! input NetCDF file ("periodic_x" and "polar" are namelist items in the &bdy_control
! section of a standard WRF "namelist.input" file), but for now we have included them
! in the "model_nml" group of DART's own "input.nml". Above, their default values are
! both set to .false.
logical :: periodic_x
logical :: periodic_y
logical :: polar
logical :: scm
integer :: domain_size
integer :: localization_coord
real(r8), dimension(:), pointer :: znu, dn, dnw, zs, znw
real(r8), dimension(:,:), pointer :: mub, hgt
real(r8), dimension(:,:), pointer :: latitude, latitude_u, latitude_v
real(r8), dimension(:,:), pointer :: longitude, longitude_u, longitude_v
real(r8), dimension(:,:,:), pointer :: phb
! NEWVAR: Currently you have to add a new type here if you want to use
! NEWVAR: a WRF variable which is not one of these types. This will go
! NEWVAR: away eventually, we hope. Search for NEWVAR for other places
! NEWVAR: the code has to change.
! JPH local variables to hold type indices
integer :: type_u, type_v, type_w, type_t, type_qv, type_qr, type_hdiab, &
type_qndrp, type_qnsnow, type_qnrain, type_qngraupel, type_qnice, &
type_qc, type_qg, type_qi, type_qs, type_gz, type_refl, type_fall_spd, &
type_dref, type_spdp, type_qh, type_qnhail, type_qhvol, type_qgvol, &
type_cldfra
integer :: type_u10, type_v10, type_t2, type_th2, type_q2, &
type_ps, type_mu, type_tsk, type_tslb, type_sh2o, &
type_smois, type_2dflash
integer :: number_of_wrf_variables
integer(i8), dimension(:,:), pointer :: var_index
integer, dimension(:,:), pointer :: var_size
integer, dimension(:), pointer :: var_type
integer, dimension(:), pointer :: var_index_list
logical, dimension(:), pointer :: var_update_list
integer, dimension(:), pointer :: dart_kind
integer, dimension(:,:), pointer :: land
real(r8), dimension(:), pointer :: lower_bound,upper_bound
character(len=10), dimension(:),pointer :: clamp_or_fail
character(len=129),dimension(:),pointer :: description, units, stagger, coordinates
end type wrf_static_data_for_dart
type wrf_dom
type(wrf_static_data_for_dart), pointer :: dom(:)
integer(i8) :: model_size
end type wrf_dom
type(wrf_dom) :: wrf
! JPH move map stuff into common (can move back into S/R later?)
real(r8) :: stdlon,truelat1,truelat2 !,latinc,loninc
! have a single, module global error string (rather than
! replicate it in each subroutine and use up more stack space)
character(len=129) :: errstring, msgstring2, msgstring3
contains
!#######################################################################
subroutine static_init_model()
! Initializes class data for WRF
integer :: ncid
integer :: io, iunit
character (len=1) :: idom
logical, parameter :: debug = .false.
integer :: ind, i, j, k, id
integer(i8) :: dart_index
integer :: my_index
integer :: var_element_list(max_state_variables)
logical :: var_update_list(max_state_variables)
real(r8) :: var_bounds_table(max_state_variables,2)
! holds the variable names for a domain when calling add_domain
character(len=129) :: netcdf_variable_names(max_state_variables)
!----------------------------------------------------------------------
! Register the module
call register_module(source, revision, revdate)
! Begin by reading the namelist input
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)
! Temporary warning until this namelist item is removed.
if (adv_mod_command /= '') then
msgstring2 = "Set the model advance command in the &dart_to_wrf_nml namelist"
call error_handler(E_MSG, 'static_init_model:', &
"WARNING: adv_mod_command ignored in &model_mod namelist", &
text2=msgstring2)
endif
allocate(wrf%dom(num_domains))
allocate(domain_id(num_domains))
! get default state variable table if asked
if ( default_state_variables ) then
wrf_state_variables = 'NULL'
call fill_default_state_table(wrf_state_variables)
msgstring2 = 'Set "default_state_variables" to .false. in the namelist'
msgstring3 = 'to use the "wrf_state_variables" list instead.'
call error_handler(E_MSG, 'static_init_model:', &
'Using predefined wrf variable list for dart state vector.', &
text2=msgstring2, text3=msgstring3)
endif
if ( debug ) then
print*,'WRF state vector table'
print*,'default_state_variables = ',default_state_variables
print*,wrf_state_variables
endif
!---------------------------
! set this array so we know exactly which obs kinds are
! allowed to be interpolated (and can give a reasonably
! helpful error message if not).
!---------------------------
num_obs_kinds = get_num_quantities()
allocate(in_state_vector(num_obs_kinds))
call fill_dart_kinds_table(wrf_state_variables, in_state_vector)
! set calendar type
call set_calendar_type(calendar_type)
! Store vertical localization coordinate
! Only 4 are allowed: level(1), pressure(2), height(3), or scale height(4)
! Everything else is assumed height
if (vert_localization_coord == VERTISLEVEL) then
wrf%dom(:)%localization_coord = VERTISLEVEL
elseif (vert_localization_coord == VERTISPRESSURE) then
wrf%dom(:)%localization_coord = VERTISPRESSURE
elseif (vert_localization_coord == VERTISHEIGHT) then
wrf%dom(:)%localization_coord = VERTISHEIGHT
elseif (vert_localization_coord == VERTISSCALEHEIGHT) then
wrf%dom(:)%localization_coord = VERTISSCALEHEIGHT
else
write(msgstring2,*)'vert_localization_coord must be one of ', &
VERTISLEVEL, VERTISPRESSURE, VERTISHEIGHT, VERTISSCALEHEIGHT
write(errstring,*)'vert_localization_coord is ', vert_localization_coord
call error_handler(E_ERR,'static_init_model', errstring, source, revision,revdate, &
text2=msgstring2)
endif
! the agreement amongst the dart/wrf users was that there was no need to
! read the wrf namelist, since the only thing it was extracting was the
! timestep, which is part of the wrf input netcdf file.
! call read_dt_from_wrf_nml()
dart_index = 1
WRFDomains : do id=1,num_domains
write( idom , '(I1)') id
! only print this once, no matter how many parallel tasks are running
if (do_output()) then
write( * ,*) '******************'
write( * ,*) '** DOMAIN # ',idom,' **'
write( * ,*) '******************'
write(logfileunit,*) '******************'
write(logfileunit,*) '** DOMAIN # ',idom,' **'
write(logfileunit,*) '******************'
endif
if(file_exist('wrfinput_d0'//idom)) then
call nc_check( nf90_open('wrfinput_d0'//idom, NF90_NOWRITE, ncid), &
'static_init_model','open wrfinput_d0'//idom )
else
call error_handler(E_ERR,'static_init_model', &
'Please put wrfinput_d0'//idom//' in the work directory.', source, revision,revdate)
endif
if(debug) write(*,*) ' ncid is ',ncid
!-------------------------------------------------------
! read WRF dimensions
!-------------------------------------------------------
call read_wrf_dimensions(ncid,wrf%dom(id)%bt, wrf%dom(id)%bts, &
wrf%dom(id)%sn, wrf%dom(id)%sns, &
wrf%dom(id)%we, wrf%dom(id)%wes, &
wrf%dom(id)%sls)
!-------------------------------------------------------
! read WRF file attributes
!-------------------------------------------------------
call read_wrf_file_attributes(ncid,id)
!-------------------------------------------------------
! assign boundary condition flags
!-------------------------------------------------------
call assign_boundary_conditions(id)
!-------------------------------------------------------
! read static data
!-------------------------------------------------------
call read_wrf_static_data(ncid,id)
!-------------------------------------------------------
! next block set up map
!-------------------------------------------------------
call setup_map_projection(id)
!-------------------------------------------------------
! end block set up map
!-------------------------------------------------------
! get the number of wrf variables wanted in this domain's state
wrf%dom(id)%number_of_wrf_variables = get_number_of_wrf_variables(id,wrf_state_variables,var_element_list, var_update_list)
! allocate and store the table locations of the variables valid on this domain
allocate(wrf%dom(id)%var_index_list(wrf%dom(id)%number_of_wrf_variables))
wrf%dom(id)%var_index_list = var_element_list(1:wrf%dom(id)%number_of_wrf_variables)
! allocation for wrf variable types
allocate(wrf%dom(id)%var_type(wrf%dom(id)%number_of_wrf_variables))
! allocation for update/nocopyback/noupdate
allocate(wrf%dom(id)%var_update_list(wrf%dom(id)%number_of_wrf_variables))
wrf%dom(id)%var_update_list = var_update_list(1:wrf%dom(id)%number_of_wrf_variables)
! allocation for dart kinds
allocate(wrf%dom(id)%dart_kind(wrf%dom(id)%number_of_wrf_variables))
! allocation of var size
allocate(wrf%dom(id)%var_size(3,wrf%dom(id)%number_of_wrf_variables))
! allocation for wrf variable metadata
allocate(wrf%dom(id)%stagger(wrf%dom(id)%number_of_wrf_variables))
allocate(wrf%dom(id)%description(wrf%dom(id)%number_of_wrf_variables))
allocate(wrf%dom(id)%units(wrf%dom(id)%number_of_wrf_variables))
allocate(wrf%dom(id)%coordinates(wrf%dom(id)%number_of_wrf_variables))
! set default bounds checking
allocate(wrf%dom(id)%lower_bound(wrf%dom(id)%number_of_wrf_variables))
allocate(wrf%dom(id)%upper_bound(wrf%dom(id)%number_of_wrf_variables))
allocate(wrf%dom(id)%clamp_or_fail(wrf%dom(id)%number_of_wrf_variables))
call set_variable_bound_defaults(wrf%dom(id)%number_of_wrf_variables, &
wrf%dom(id)%lower_bound, &
wrf%dom(id)%upper_bound, &
wrf%dom(id)%clamp_or_fail)
! build the variable indices
! this accounts for the fact that some variables might not be on all domains
do ind = 1,wrf%dom(id)%number_of_wrf_variables
! actual location in state variable table
my_index = wrf%dom(id)%var_index_list(ind)
wrf%dom(id)%var_type(ind) = ind ! types are just the order for this domain
wrf%dom(id)%dart_kind(ind) = get_index_for_quantity(trim(wrf_state_variables(2,my_index)))
if ( debug ) then
print*,'dart kind identified: ',trim(wrf_state_variables(2,my_index)),' ',wrf%dom(id)%dart_kind(ind)
endif
! get stagger and variable size
call get_variable_size_from_file(ncid,id, &
wrf_state_variables(1,my_index), &
wrf%dom(id)%bt, wrf%dom(id)%bts, &
wrf%dom(id)%sn, wrf%dom(id)%sns, &
wrf%dom(id)%we, wrf%dom(id)%wes, &
wrf%dom(id)%stagger(ind), &
wrf%dom(id)%var_size(:,ind))
! get other variable metadata; units, coordinates and description
call get_variable_metadata_from_file(ncid,id, &
wrf_state_variables(1,my_index), &
wrf%dom(id)%description(ind), &
wrf%dom(id)%coordinates(ind), &
wrf%dom(id)%units(ind) )
if ( debug ) then
print*,'variable size ',trim(wrf_state_variables(1,my_index)),' ',wrf%dom(id)%var_size(:,ind)
endif
! add bounds checking information
call get_variable_bounds(wrf_state_bounds, wrf_state_variables(1,my_index), &
wrf%dom(id)%lower_bound(ind), wrf%dom(id)%upper_bound(ind), &
wrf%dom(id)%clamp_or_fail(ind))
if ( debug ) then
write(*,*) 'Bounds for variable ', &
trim(wrf_state_variables(1,my_index)), &
' are ',wrf%dom(id)%lower_bound(ind), &
wrf%dom(id)%upper_bound(ind), &
wrf%dom(id)%clamp_or_fail(ind)
endif
write(errstring, '(A,I4,2A)') 'state vector array ', ind, ' is ', trim(wrf_state_variables(1,my_index))
call error_handler(E_MSG, 'static_init_model: ', errstring)
enddo
if (do_output()) then
write( * ,*)
write(logfileunit,*)
endif
! close data file, we have all we need
call nc_check(nf90_close(ncid),'static_init_model','close wrfinput_d0'//idom)
! indices into 1D array - hopefully this becomes obsolete
! JPH changed last dimension here from num_model_var_types
!HK allocate(wrf%dom(id)%dart_ind(wrf%dom(id)%wes,wrf%dom(id)%sns,wrf%dom(id)%bts,wrf%dom(id)%number_of_wrf_variables))
! start and stop of each variable in vector
allocate(wrf%dom(id)%var_index(2,wrf%dom(id)%number_of_wrf_variables))
!---------------------------
! end block to be obsolete
!---------------------------
!---------------------------
! at this point we need all information assigned to each variable
! then just loop thru the table
!---------------------------
!HK wrf%dom(id)%dart_ind = 0
! NOTE: this could be put into the loop above if wrf%dom(id)%dart_ind is
! eventually not needed
! Here we use ind instead of type as the 4th dimension. In is linked to the
! specific type via wrf%dom(id)%var_index_list(ind). This saves some
! space from the previous implementation but I am not yet sure of other
! problems that it might cause.
do ind = 1,wrf%dom(id)%number_of_wrf_variables
my_index = wrf%dom(id)%var_index_list(ind)
if ( debug ) then
write(*,*) 'Assigning dart vector indices for var_type ',wrf%dom(id)%var_type(ind)
write(*,*) 'affiliated with WRF variable ',trim(wrf_state_variables(1,my_index)),' of size ',wrf%dom(id)%var_size(:,ind)
endif
wrf%dom(id)%var_index(1,ind) = dart_index
do k=1,wrf%dom(id)%var_size(3,ind)
do j=1,wrf%dom(id)%var_size(2,ind)
do i=1,wrf%dom(id)%var_size(1,ind)
! HK wrf%dom(id)%dart_ind(i,j,k,ind) &
! = dart_index
dart_index = dart_index + 1
enddo
enddo
enddo
wrf%dom(id)%var_index(2,ind) = dart_index - 1
if ( debug ) write(*,*) 'assigned start, stop ',wrf%dom(id)%var_index(:,ind)
enddo ! loop through all viable state variables on this domain
if ( id == 1 ) then
wrf%dom(id)%domain_size = dart_index - 1
else
wrf%dom(id)%domain_size = dart_index - 1
do ind = id-1, 1, -1
wrf%dom(id)%domain_size = wrf%dom(id)%domain_size - wrf%dom(ind)%domain_size
enddo
endif
! NEWVAR: If you add a new wrf array type which is not yet in this list, currently
! NEWVAR: you will have to add it here, and add a type_xx for it, and also add
! NEWVAR: a in_state_vector case in the select statement. search for NEWVAR.
! JPH now that we have the domain ID just go ahead and get type indices once
! NOTE: this is not strictly necessary - can use only stagger info in the future (???)
wrf%dom(id)%type_u = get_type_ind_from_type_string(id,'U')
wrf%dom(id)%type_v = get_type_ind_from_type_string(id,'V')
wrf%dom(id)%type_w = get_type_ind_from_type_string(id,'W')
wrf%dom(id)%type_t = get_type_ind_from_type_string(id,'T')
wrf%dom(id)%type_gz = get_type_ind_from_type_string(id,'PH')
wrf%dom(id)%type_qv = get_type_ind_from_type_string(id,'QVAPOR')
wrf%dom(id)%type_qr = get_type_ind_from_type_string(id,'QRAIN')
wrf%dom(id)%type_qc = get_type_ind_from_type_string(id,'QCLOUD')
wrf%dom(id)%type_qg = get_type_ind_from_type_string(id,'QGRAUP')
wrf%dom(id)%type_qh = get_type_ind_from_type_string(id,'QHAIL')
wrf%dom(id)%type_qi = get_type_ind_from_type_string(id,'QICE')
wrf%dom(id)%type_qs = get_type_ind_from_type_string(id,'QSNOW')
wrf%dom(id)%type_cldfra = get_type_ind_from_type_string(id,'CLDFRA')
wrf%dom(id)%type_qgvol = get_type_ind_from_type_string(id,'QVGRAUPEL')
wrf%dom(id)%type_qhvol = get_type_ind_from_type_string(id,'QVHAIL')
wrf%dom(id)%type_qnice = get_type_ind_from_type_string(id,'QNICE')
wrf%dom(id)%type_qndrp = get_type_ind_from_type_string(id,'QNDRP')
wrf%dom(id)%type_qnsnow = get_type_ind_from_type_string(id,'QNSNOW')
wrf%dom(id)%type_qnrain = get_type_ind_from_type_string(id,'QNRAIN')
wrf%dom(id)%type_qngraupel = get_type_ind_from_type_string(id,'QNGRAUPEL')
wrf%dom(id)%type_qnhail = get_type_ind_from_type_string(id,'QNHAIL')
wrf%dom(id)%type_u10 = get_type_ind_from_type_string(id,'U10')
wrf%dom(id)%type_v10 = get_type_ind_from_type_string(id,'V10')
wrf%dom(id)%type_t2 = get_type_ind_from_type_string(id,'T2')
wrf%dom(id)%type_th2 = get_type_ind_from_type_string(id,'TH2')
wrf%dom(id)%type_q2 = get_type_ind_from_type_string(id,'Q2')
wrf%dom(id)%type_ps = get_type_ind_from_type_string(id,'PSFC')
wrf%dom(id)%type_mu = get_type_ind_from_type_string(id,'MU')
wrf%dom(id)%type_tsk = get_type_ind_from_type_string(id,'TSK')
wrf%dom(id)%type_2dflash = get_type_ind_from_type_string(id,'FLASH_RATE_2D')
wrf%dom(id)%type_tslb = get_type_ind_from_type_string(id,'TSLB')
wrf%dom(id)%type_smois = get_type_ind_from_type_string(id,'SMOIS')
wrf%dom(id)%type_sh2o = get_type_ind_from_type_string(id,'SH2O')
wrf%dom(id)%type_refl = get_type_ind_from_type_string(id,'REFL_10CM')
wrf%dom(id)%type_dref = get_type_ind_from_type_string(id,'DIFF_REFL_10CM')
wrf%dom(id)%type_spdp = get_type_ind_from_type_string(id,'SPEC_DIFF_10CM')
wrf%dom(id)%type_fall_spd = get_type_ind_from_type_string(id,'FALL_SPD_Z_WEIGHTED')
!wrf%dom(id)%type_fall_spd = get_type_ind_from_type_string(id,'VT_DBZ_WT')
wrf%dom(id)%type_hdiab = get_type_ind_from_type_string(id,'H_DIABATIC')
! variable bound table for setting upper and lower bounds of variables
var_bounds_table(1:wrf%dom(id)%number_of_wrf_variables,1) = wrf%dom(id)%lower_bound
var_bounds_table(1:wrf%dom(id)%number_of_wrf_variables,2) = wrf%dom(id)%upper_bound
! List of netcdf variable names in the domain
do i = 1, wrf%dom(id)%number_of_wrf_variables
my_index = wrf%dom(id)%var_index_list(i) ! index in wrf_state_variables
netcdf_variable_names(i) = wrf_state_variables(1, my_index)
enddo
! add domain - not doing anything with domain_id yet so just overwriting it
domain_id(id) = add_domain( 'wrfinput_d0'//idom, &
wrf%dom(id)%number_of_wrf_variables, &
var_names = netcdf_variable_names(1:wrf%dom(id)%number_of_wrf_variables), &
kind_list = wrf%dom(id)%dart_kind, &
clamp_vals = var_bounds_table(1:wrf%dom(id)%number_of_wrf_variables,:), &
update_list = var_update_list(1:wrf%dom(id)%number_of_wrf_variables) )
if (debug) call state_structure_info(domain_id(id))
enddo WRFDomains
wrf%model_size = dart_index - 1
write(errstring,*) ' wrf model size is ',wrf%model_size
call error_handler(E_MSG, 'static_init_model', errstring)
! tell the location module how we want to localize in the vertical
call set_vertical_localization_coord(vert_localization_coord)
end subroutine static_init_model
!#######################################################################
!> Convert from the state structure id to the wrf domain number.
!> These are the same if there is only WRF involved in the assimilation
!> The state structure id may be different if WRF is coupled with another
!> model.
function get_wrf_domain(state_id)
integer, intent(in) :: state_id
integer :: get_wrf_domain
integer :: i
do i = 1, num_domains
if (domain_id(i) == state_id) then
get_wrf_domain = i
return
endif
enddo
call error_handler(E_ERR, 'get_wrf_domain', 'not a valid domain')
end function get_wrf_domain
!#######################################################################
function get_model_size()
integer(i8) :: get_model_size
get_model_size = wrf%model_size
end function get_model_size
!#######################################################################
function get_number_domains()
integer :: get_number_domains
get_number_domains = num_domains
end function get_number_domains
!#######################################################################
function get_wrf_static_data(dom_num)
integer, intent(in) :: dom_num
type(wrf_static_data_for_dart) :: get_wrf_static_data
get_wrf_static_data = wrf%dom(dom_num)
return
end function get_wrf_static_data
!#######################################################################
subroutine get_wrf_state_variables(state_var)
character(len=129), intent(out) :: state_var(num_state_table_columns,max_state_variables)
state_var = wrf_state_variables
end subroutine get_wrf_state_variables
!#######################################################################
function shortest_time_between_assimilations()
type(time_type) :: shortest_time_between_assimilations
integer :: model_dt, assim_dt
! We need to coordinate the desired assimilation window to be a
! multiple of the model time step (which has no precision past integer seconds).
model_dt = nint(wrf%dom(1)%dt)
! The integer arithmetic does its magic.
assim_dt = (assimilation_period_seconds / model_dt) * model_dt
shortest_time_between_assimilations = set_time(assim_dt)
end function shortest_time_between_assimilations
!#######################################################################
subroutine get_state_meta_data(index_in, location, var_type_out, id_out)
! Given an integer index into the DART state vector structure, returns the
! associated location. This is not a function because the more general
! form of the call has a second intent(out) optional argument kind.
! Maybe a functional form should be added?
! this version has an optional last argument that is never called by
! any of the dart code, which can return the wrf domain number.
integer(i8), intent(in) :: index_in
type(location_type), intent(out) :: location
integer, optional, intent(out) :: var_type_out, id_out
integer :: var_type, dart_type
integer(i8) :: index
integer :: ip, jp, kp
integer :: nz, ny, nx
logical :: var_found
real(r8) :: lon, lat, lev
character(len=129) :: string1
integer :: i, id, var_id, state_id
logical, parameter :: debug = .false.
! from the dart index get the local variables indices
call get_model_variable_indices(index_in, ip, jp, kp, var_id=var_id, dom_id=state_id)
! convert from state_structure domain number to wrf.
id = get_wrf_domain(state_id)
! at this point, (ip,jp,kp) refer to indices in the variable's own grid
if(debug) write(*,*) ' ip, jp, kp for index ',ip,jp,kp,index
if(debug) write(*,*) ' Var type: ',var_type
var_type = wrf%dom(id)%var_type(var_id)
dart_type = wrf%dom(id)%dart_kind(var_id)
! first obtain lat/lon from (ip,jp)
call get_wrf_horizontal_location( ip, jp, var_type, id, lon, lat )
! return staggered levels for vertical types? check this later.
if( (var_type == wrf%dom(id)%type_w ) .or. (var_type == wrf%dom(id)%type_gz) ) then
lev = real(kp) - 0.5_r8
else
lev = real(kp)
endif
if(debug) write(*,*) 'lon, lat, lev: ',lon, lat, lev
location = set_location(lon, lat, lev, VERTISLEVEL)
! return DART variable kind if requested
if(present(var_type_out)) var_type_out = dart_type
! return domain id if requested
if(present(id_out)) id_out = id
end subroutine get_state_meta_data
!--------------------------------------------------------------------
!> Distributed version of model interpolate
!> obs_kind is called as location type in assim model model
!> Storing the mean copy level location as the observation location
!> to save recomputation of v_p in vert_convert in get_close_obs
!
! Should this code be simplified so there is not so much repetition?
! This is the main forward operator subroutine for WRF.
! Given an ob (its DART location and kind), the corresponding model
! value is computed at nearest i,j,k. Thus, first i,j,k is obtained
! from ob lon,lat,z and then the state value that corresponds to
! the ob kind is interpolated.
!
! No location conversions are carried out in this subroutine. See
! get_close_obs, where ob vertical location information is converted
! to the requested vertical coordinate type.
subroutine model_interpolate(state_handle, ens_size, location, obs_kind, expected_obs, istatus)
! x: Full DART state vector relevant to what's being updated
! in the filter (mean or individual members).
! istatus: Returned 0 if everything is OK, >0 if error occured.
! 1 = missing domain id
! 2 = bad vertical level
! 3 = unsupported obs kind
! 10 = polar observation while restrict_polar namelist true
! 99 = unknown reason (reached end of routine with missing_r8
! as obs_val)
! modified 26 June 2006 to accomodate vortex attributes
! modified 13 December 2006 to accomodate changes for the mpi version
! modified 22 October 2007 to accomodate global WRF (3.0)
! modified 18 November 2008 to allow unknown kinds to return without stopping
! modified 5 February 2010 to add circulation calc for vortex obs
! Helen Kershaw - Aim: to not require the whole state vector
! arguments
type(ensemble_type), intent(in) :: state_handle
integer, intent(in) :: ens_size
type(location_type), intent(in) :: location
integer, intent(in) :: obs_kind
real(r8), intent(out) :: expected_obs(ens_size)
integer, intent(out) :: istatus(ens_size)
! local
logical, parameter :: debug = .false.
logical, parameter :: restrict_polar = .false.
logical, parameter :: use_old_vortex = .true. ! set to .false. to use circ code
real(r8), parameter :: drad = pi / 18.0_r8
real(r8) :: xloc, yloc, xloc_u, yloc_v, xyz_loc(3)
integer :: i, i_u, j, j_v, k2
real(r8) :: dx,dy,dxm,dym,dx_u,dxm_u,dy_v,dym_v
integer :: id
logical :: surf_var
real(r8) :: a1(ens_size)
real(r8) :: zloc(ens_size)
integer :: k(ens_size)
real(r8) :: dz(ens_size), dzm(ens_size)
real(r8) :: utrue(ens_size), vtrue(ens_size)
! from getCorners
integer, dimension(2) :: ll, lr, ul, ur, ll_v, lr_v, ul_v, ur_v
integer :: rc, i1, i2
integer(i8) :: ill, ilr, iul, iur
real(r8) :: fld(2,ens_size)
real(r8), allocatable, dimension(:,:) :: v_h, v_p
! local vars, used in finding sea-level pressure and vortex center
real(r8), allocatable, dimension(:, :) :: t1d, p1d, qv1d, z1d
real(r8), allocatable, dimension(:, :) :: z1d_1, z1d_2
real(r8), allocatable, dimension(:,:) :: pp, pd
real(r8), allocatable, dimension(:,:,:) :: vfld
real(r8), allocatable, dimension(:,:,:) :: uwnd, vwnd, vort
real(r8), allocatable, dimension(:) :: x1d, y1d, xx1d, yy1d
integer :: center_search_half_size, center_track_xmin, center_track_ymin, &
center_track_xmax, center_track_ymax, circ_half_size, &
circ_xmin, circ_xmax, circ_ymin, circ_ymax, xlen, ylen, &
xxlen, yylen, ii1, ii2, cxlen, cylen, imax, jmax
real(r8) :: clat, clon, cxloc, cyloc, vcrit, &
circ_half_length, asum, distgrid, dgi1, dgi2
real(r8) :: magwnd(ens_size), maxwspd(ens_size), circ(ens_size)
! local vars, used in calculating density, pressure, height
real(r8) :: rho1(ens_size) , rho2(ens_size) , rho3(ens_size), rho4(ens_size)
real(r8) :: pres1(ens_size), pres2(ens_size), pres3(ens_size), pres4(ens_size), pres(ens_size)