-
Notifications
You must be signed in to change notification settings - Fork 4
/
dataassim_user.f
1838 lines (1669 loc) · 69.2 KB
/
dataassim_user.f
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
c dataassim_user.f
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine DATAASSIM_USER(nx,ny,icorr_factor_index,
& corr_factor,max_iter,deltax,deltay,xmn,ymn,nobs_dates,
& print_inc,iday_init,imonth_init,iyear_init,depth_assim,
& ihrestart_flag)
c Perform the required correction (precipitation and melt) factor
c calculations. To run the data assimilation routines requires
c additional model inputs to be added here (inputs in addition
c to those in the snowmodel.par file). Then you must recompile
c the code before running it.
c This program works when you have data at one or many points
c (many individual grid cells), for one or more times. And for
c data (e.g., average swe) over areas of grid cells; there can be
c many of these areas, at many different times.
c The input (swe) data file is assumed to be in /snowmodel/data/.
c See below for the details of how this file is to be configured.
c The data assimilation produces a couple of extra files the user
c can look at to study what the model did with the input data
c file (fname_sweobs) that was provided. First, there is a text
c file written to /snowmodel/data/ that provides a summary of the
c calculations that were done (see unit=77 in the code below for
c what is written to this file). Second, there is a copy of the
c precipitation correction surfaces that were generated. This
c is also placed in /snowmodel/data/, the file is a GrADS file
c (called corr_factor.gdat), and there is also a corr_factor.ctl
c there that can be modified to fit the current simulation. The
c data layers in corr_factor.gdat are structured as follows:
c The number of layers equals the total number of observation
c times that were assimilated, plus the number of years in the
c assimilation. The order is: the correction surface (cf) for
c the time between the start of the simulation and observation
c time 1 in year 1, the cf for the time between obs time 1 in
c year 1 and obs time 2 in year 1, etc., then a cf==1.0 for the
c time between the last obs time and the end of year 1 (or the
c end of the simulation for a 1-year run). Then this order
c repeats for the rest of the years of the simulation. In the
c GrADS control file (corr_factor.ctl) these layers are assumed
c to correspond to different times in the data file (although
c the actual time increment defined in the .ctl file is not
c really relevant: for example, t=1 corresponds to the cf for
c obs 1, t=2 is for obs 2, t=3 is for the time between the last
c obs time and the end of year 1, t=4 is for the obs 1 in year
c 2, etc.).
implicit none
include 'snowmodel.inc'
real deltax,deltay,beta,areas_flag,print_inc
real sprec_ratio(max_obs_dates),smelt_ratio(max_obs_dates)
real corr_factor(nx_max,ny_max,max_obs_dates)
real areas_mask(nx_max,ny_max)
double precision xmn,ymn
double precision xstn(nx_max*ny_max),ystn(nx_max*ny_max)
integer icorr_factor_index(max_time_steps)
integer iobs_rec(max_obs_dates)
integer nobs_dates,nx,ny,max_iter,local_assim_flag,iday_init,
& imonth_init,iyear_init,nyear,nyears,nobs_total,nobs_total_cfi,
& ihrestart_flag
integer depth_assim
character*80 fname_swed,fname_sspr,fname_ssmt,fname_sden
character*80 fname_sweobs
character*80 fname_sweobs_barnes_mask
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c BEGIN USER EDIT SECTION.
c BEGIN USER EDIT SECTION.
c Define how many years there are in this simulation.
nyears = 1
c A single data file describes the observation information that
c will be used in the data assimilation. This file contains the
c following information in the following format. The id can be
c any number, x and y are easting and northing in m, and swe is
c in m).
c total_number_of_observation_dates_for_this_year
c iyr imo idy (for this observation date)
c number_of_stations_for_this_observation_date
c id x y swe
c id x y swe
c id x y swe
c iyr imo idy (for this observation date)
c number_of_stations_for_this_observation_date
c id x y swe
c id x y swe
c For example:
c 2
c 2014 3 15
c 3
c 101 3456.7 23677.4 0.42
c 102 3556.3 25079.3 0.52
c 103 3106.2 29089.3 0.59
c 2014 4 1
c 2
c 101 3456.7 23677.4 0.48
c 103 3106.2 29089.3 0.62
c Then this repeats for each year of the assimilation (the input
c file looks like the above for a single-year run, and for
c multi-year runs the data for each following year is just
c stacked on top of the previous year (like the example below
c for a two-year assimilation run). The code requires a
c "total_number_of_observation_dates_for_this_year"
c line for each year. If you have a year with no data, this
c can be set to be 0 (zero). If this is 0, the code sets the
c correction factor to equal 1.0 for that simulation year, and
c no adjustments are made to the precipitation for that year.
c 2
c 2014 3 15
c 3
c 101 3456.7 23677.4 0.42
c 102 3556.3 25079.3 0.52
c 103 3106.2 29089.3 0.59
c 2014 4 1
c 2
c 101 3456.7 23677.4 0.48
c 103 3106.2 29089.3 0.62
c 1
c 2015 3 25
c 2
c 101 3456.7 23677.4 0.23
c 102 3556.3 25079.3 0.32
c Provide the name of the data file that contains the observed swe
c information (as described above).
c fname_sweobs = 'data/Full_pkSWE.txt'
fname_sweobs = 'data/SDVEl_X276562.5_pkSWEbl.txt'
c Define the file names of the swe depth (swed), annual summed snow
c precipitation (sspr), and annual summed snowmelt (ssmt) outputs
c from the first iteration of the data assimilation run. In this
c implementation of the data assimilation code, I have assumed
c that the output files are those created by outputs_user.f,
c where there is an individual file for each variable.
fname_swed = 'outputs/wo_assim/swed.gdat'
fname_sspr = 'outputs/wo_assim/sspr.gdat'
fname_ssmt = 'outputs/wo_assim/ssmt.gdat'
fname_sden = 'outputs/wo_assim/sden.gdat'
c THE PARAMETERS BELOW ARE RARELY CHANGED, UNLESS YOU ARE DOING AN
c AREAS ASSIMILATION (INSTEAD OF ASSIMILATING POINT DATA).
c Beta controls the interpolation distance weights. Beta = 1.0
c will give you a very smooth field, and correction factor
c distributions that may not produce swe's that exactly match
c the observations. Beta << 1.0 will give you correction factor
c fields that go right through the data. If you just have one
c data point/area, beta is not used.
beta = 1.0
c beta = 0.1
c beta = 0.5
c Define whether this simulation will be processing areas (data
c within groups of grid cells: areas_flag = 1.0), or points
c (single grid cells: areas_flag = 0.0). Note that if you have
c a combination of areas and points, you have to use the areas
c option and treat each point like a single-grid-cell (small)
c area.
areas_flag = 0.0
c areas_flag = 1.0
c If this is an areas simulation, open and read in the areas mask
c data. Note that here I assume that the area mask is a nx by ny
c file with undef values everywhere except at the area 'stations'.
c And that each 'station' area is given a 1.0, 2.0, etc. that
c corresponds to the order of the station listing in the 'station'
c data input file (the first 'station' listed has mask value = 1.0,
c the second listed has mask value = 2.0, etc.
if (areas_flag.eq.1.0) then
c open(63,file='sweobs/nea.obsmask.100m.2009.gdat',
c & form='unformatted',access='direct',recl=4*nx*ny)
c read(63,rec=1) ((areas_mask(i,j),i=1,nx),j=1,ny)
c If you have two masks for two different observation dates, then
c do something like the following.
c open(63,file='data/zack_obs_mask.gdat',
c & form='unformatted',access='direct',recl=4*nx*ny)
c read(63,rec=1) ((areas_mask(i,j,1),i=1,nx),j=1,ny)
c read(63,rec=2) ((areas_mask(i,j,2),i=1,nx),j=1,ny)
endif
c Define whether this simulation is going to restrict the
c assimilation influence to some local area surrounding each
c data point that is assimilated. This was implemented for
c ANWR simulations where we only had observations in a corner
c of the simulation domain and I didn't want the corrections
c to extend too far outside that local observation area. So,
c this is an example of what can be done, and it is not written
c for general application. If you want to do something similar,
c the associated subroutine can be edited for your specific
c simulation of interest. For yes, local_assim_flag = 1, for
c no, local_assim_flag = 0.
local_assim_flag = 0
c Identify the file that contains the local data assimilation
c mask. This is only used if local_assim_flag = 1.
fname_sweobs_barnes_mask =
& '../swe_obs/2014/barnes/obs.gridded.gdat'
c END USER EDIT SECTION.
c END USER EDIT SECTION.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c LOOP THROUGH THE YEARS IN THE SIMULATION.
do nyear=1,nyears
c Run the data assimilation routines.
call data_assimilation(nx,ny,deltax,deltay,beta,
& areas_flag,sprec_ratio,smelt_ratio,corr_factor,
& areas_mask,xmn,ymn,xstn,ystn,nobs_dates,iobs_rec,
& local_assim_flag,fname_swed,fname_sspr,fname_ssmt,
& fname_sweobs,fname_sweobs_barnes_mask,iday_init,
& imonth_init,iyear_init,nyear,nobs_total,depth_assim,
& fname_sden)
c Build an array indicating the appropriate correction factor to
c use at any given time during the simulation. What this does
c is define an index array that contains the record number that
c gets used at every model time step during the second model run
c loop. This record number corresponds to the record (krec) of
c the corr_factor(i,j,krec) array that was generated and saved
c in the subroutine above.
call corr_factor_index(nobs_dates,icorr_factor_index,
& iobs_rec,max_iter,sprec_ratio,smelt_ratio,print_inc,
& iday_init,imonth_init,iyear_init,nyear,nobs_total_cfi,
& nyears,ihrestart_flag)
enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine data_assimilation(nx,ny,deltax,deltay,beta,
& areas_flag,sprec_ratio,smelt_ratio,corr_factor,
& areas_mask,xmn,ymn,xstn,ystn,nobs_dates,iobs_rec,
& local_assim_flag,fname_swed,fname_sspr,fname_ssmt,
& fname_sweobs,fname_sweobs_barnes_mask,iday_init,
& imonth_init,iyear_init,nyear,nobs_total,depth_assim,
& fname_sden)
implicit none
include 'snowmodel.inc'
real deltax,deltay,undef,dn,beta,areas_flag,swe_count,
& cf_min,ro_avg
real sprec_ratio(max_obs_dates),smelt_ratio(max_obs_dates)
real corr_factor(nx_max,ny_max,max_obs_dates)
real swe_tmp(nx_max,ny_max),sum_sprec_tmp1(nx_max,ny_max),
& sum_sprec_tmp2(nx_max,ny_max),grid(nx_max,ny_max),
& areas_mask(nx_max,ny_max),sum_smelt_tmp1(nx_max,ny_max),
& sum_smelt_tmp2(nx_max,ny_max),ro_tmp(nx_max,ny_max)
real corr_factor_tmp(nx_max*ny_max),swe_obs(nx_max*ny_max),
& swe_model(nx_max*ny_max),sumsprec_model(nx_max*ny_max),
& delta_old(nx_max*ny_max),obsid(nx_max*ny_max),
& sumsmelt_model(nx_max*ny_max),obsid_old(nx_max*ny_max),
& delta_old_tmp(nx_max*ny_max),depth_obs(nx_max*ny_max),
& ro_model(nx_max*ny_max)
c real corr_offset(nx_max*ny_max)
double precision xmn,ymn
double precision xstn(nx_max*ny_max),ystn(nx_max*ny_max)
integer iobs_rec(max_obs_dates)
integer ii(nx_max*ny_max),jj(nx_max*ny_max)
integer iobs_num,irec1,irec2,nobs_dates,nx,ny,i,j,ifill,
& iobsint,k,nstns,nstns_old,kk,local_assim_flag,iiyr,iimo,
& iidy,iobs_rec_tmp,iday_init,imonth_init,iyear_init,nyear,
& krec,nobs_total
integer depth_assim,count
character*80 fname_swed,fname_sspr,fname_ssmt,fname_sden
character*80 fname_sweobs
character*80 fname_sweobs_barnes_mask
c Perform some initialization steps.
if (nyear.eq.1) then
c Define some of the constants and parameters used in the data
c assimilation. ifill should be = 1; in that case undef is not
c really used (so it does not have to be the same as defined
c in the .par file).
undef = -9999.0
ifill = 1
iobsint = 0
c Open the swe input file.
open (unit=61,file=fname_sweobs)
c Open a file to write some basic correction factor information.
c This just saves information that the user might want to look
c at.
open (unit=77,file='data/corr_factor.txt')
c Open an output file for the correction factor array.
open(62,file='data/corr_factor.gdat',
& form='unformatted',access='direct',recl=4*nx*ny)
c Initialize the number-of-observations counter.
nobs_total = 0
endif
c Read the number of observation dates for this year.
read(61,*) nobs_dates
c If you have observations for this year, generate the correction
c factors. For the case of no observations for this year, set
c the correction factor equal to 1.0.
if (nobs_dates.gt.0) then
c Loop through the observation dates.
do iobs_num=1,nobs_dates
c Increment the number-of-observations counter.
nobs_total = nobs_total + 1
c Read the date corresponding to this observation.
read(61,*) iiyr,iimo,iidy
c Convert this date to the corresponding record number in the
c original SnowModel output files. Note that this has assumed
c that daily data files were written out since the start of
c the simulation.
call get_obs_record(iday_init,imonth_init,iyear_init,
& iidy,iimo,iiyr,iobs_rec_tmp)
iobs_rec(iobs_num) = iobs_rec_tmp
c For this observation date, read in the data describing the
c location and swe values for each observation. For areas
c simulations, xstn, and ystn correspond to the center of the
c area domain and they are not really used.
read(61,*) nstns
do k=1,nstns
if (depth_assim.eq.1) then
read(61,*) obsid(k),xstn(k),ystn(k),depth_obs(k)
else
read(61,*) obsid(k),xstn(k),ystn(k),swe_obs(k)
endif
enddo
c Convert the x and y locations to (ii,jj) locations.
do k=1,nstns
ii(k) = 1 + nint((xstn(k) - xmn) / deltax)
jj(k) = 1 + nint((ystn(k) - ymn) / deltay)
enddo
c If you do a data assimilation run from start to finish, it is
c not required to close and reopen these files. But if you are
c doing a history restart then these files are no longer open
c so you must do this.
close (238)
close (239)
close (240)
close (237)
c Open the required inputs from the initial assimilation loop.
c Open swe depth (swe_depth).
c /outputs/wo_assim/swed.gdat is unit 238 in outputs_user.f
c Open sum snow precip (sum_sprec).
c /outputs/wo_assim/sspr.gdat is unit 239 in outputs_user.f
c Open sum snow melt (sum_smelt).
c /outputs/wo_assim/ssmt.gdat is unit 240 in outputs_user.f
open (238,file=fname_swed,
& form='unformatted',access='direct',recl=4*1*nx*ny)
open (239,file=fname_sspr,
& form='unformatted',access='direct',recl=4*1*nx*ny)
open (240,file=fname_ssmt,
& form='unformatted',access='direct',recl=4*1*nx*ny)
print *,fname_ssmt
print *,fname_sden
open (237,file=fname_sden,
& form='unformatted',access='direct',recl=4*1*nx*ny)
c Read the model output for the first observation time.
if (iobs_num.eq.1) then
irec1 = iobs_rec(iobs_num)
read(238,rec=irec1) ((swe_tmp(i,j),i=1,nx),j=1,ny)
read(239,rec=irec1) ((sum_sprec_tmp1(i,j),i=1,nx),j=1,ny)
read(240,rec=irec1) ((sum_smelt_tmp1(i,j),i=1,nx),j=1,ny)
read(237,rec=irec1) ((ro_tmp(i,j),i=1,nx),j=1,ny)
count = 0
ro_avg = 0
do i=1,nx
do j=1,ny
if (ro_tmp(i,j).ge.0) then
ro_avg = ro_avg + ro_tmp(i,j)
count = count + 1
endif
enddo
enddo
if (count.gt.0) ro_avg = ro_avg/count
c print *,swe_tmp(1,1), sum_sprec_tmp1(1,1),
c & sum_smelt_tmp1(1,1)
c stop
c For points, just pull the data at the appropriate grid cell.
c For areas, average the data over the masked out area for each
c 'station'.
do k=1,nstns
if (areas_flag.eq.0.0) then
swe_model(k) = swe_tmp(ii(k),jj(k))
ro_model(k) = ro_tmp(ii(k),jj(k))
sumsprec_model(k) = sum_sprec_tmp1(ii(k),jj(k))
sumsmelt_model(k) = sum_smelt_tmp1(ii(k),jj(k))
elseif (areas_flag.eq.1.0) then
swe_model(k) = 0.0
ro_model(k) = 0.0
sumsprec_model(k) = 0.0
sumsmelt_model(k) = 0.0
swe_count = 0.0
do j=1,ny
do i=1,nx
if (areas_mask(i,j).eq.obsid(k)) then
c The following is used if the mask changes with observation time.
c if (areas_mask(i,j,iobs_num).eq.obsid(k)) then
swe_count = swe_count + 1.0
swe_model(k) = swe_model(k) + swe_tmp(i,j)
ro_model(k) = ro_model(k) + ro_tmp(i,j)
sumsprec_model(k) = sumsprec_model(k) +
& sum_sprec_tmp1(i,j)
sumsmelt_model(k) = sumsmelt_model(k) +
& sum_smelt_tmp1(i,j)
endif
enddo
enddo
swe_model(k) = swe_model(k) / swe_count
ro_model(k) = ro_model(k) / swe_count
sumsprec_model(k) = sumsprec_model(k) / swe_count
sumsmelt_model(k) = sumsmelt_model(k) / swe_count
endif
enddo
endif
c Read the model output for any additional observation times (irec1
c = current obs time, irec2 = previous obs time).
if (iobs_num.gt.1) then
irec1 = iobs_rec(iobs_num)
irec2 = iobs_rec(iobs_num-1)
read(238,rec=irec1) ((swe_tmp(i,j),i=1,nx),j=1,ny)
read(239,rec=irec1) ((sum_sprec_tmp1(i,j),i=1,nx),j=1,ny)
read(239,rec=irec2) ((sum_sprec_tmp2(i,j),i=1,nx),j=1,ny)
read(240,rec=irec1) ((sum_smelt_tmp1(i,j),i=1,nx),j=1,ny)
read(240,rec=irec2) ((sum_smelt_tmp2(i,j),i=1,nx),j=1,ny)
read(237,rec=irec1) ((ro_tmp(i,j),i=1,nx),j=1,ny)
count = 0
ro_avg = 0
do i=1,nx
do j=1,ny
if (ro_tmp(i,j).ge.0) then
ro_avg = ro_avg + ro_tmp(i,j)
count = count + 1
endif
enddo
enddo
if (count.gt.0) ro_avg = ro_avg/count
c For points, just pull the data at the appropriate grid cell.
c For areas, average the data over the masked out area for each
c 'station'.
do k=1,nstns
if (areas_flag.eq.0.0) then
swe_model(k) = swe_tmp(ii(k),jj(k))
ro_model(k) = ro_tmp(ii(k),jj(k))
sumsprec_model(k) = sum_sprec_tmp1(ii(k),jj(k)) -
& sum_sprec_tmp2(ii(k),jj(k))
sumsmelt_model(k) = sum_smelt_tmp1(ii(k),jj(k)) -
& sum_smelt_tmp2(ii(k),jj(k))
elseif (areas_flag.eq.1.0) then
swe_model(k) = 0.0
ro_model(k) = 0.0
sumsprec_model(k) = 0.0
sumsmelt_model(k) = 0.0
swe_count = 0.0
do j=1,ny
do i=1,nx
if (areas_mask(i,j).eq.obsid(k)) then
c The following is used if the mask changes with observation time.
c if (areas_mask(i,j,iobs_num).eq.obsid(k)) then
swe_count = swe_count + 1.0
swe_model(k) = swe_model(k) + swe_tmp(i,j)
ro_model(k) = ro_model(k) + ro_tmp(i,j)
sumsprec_model(k) = sumsprec_model(k) +
& sum_sprec_tmp1(i,j) - sum_sprec_tmp2(i,j)
sumsmelt_model(k) = sumsmelt_model(k) +
& sum_smelt_tmp1(i,j) - sum_smelt_tmp2(i,j)
endif
enddo
enddo
swe_model(k) = swe_model(k) / swe_count
ro_model(k) = ro_model(k) / swe_count
sumsprec_model(k) = sumsprec_model(k) / swe_count
sumsmelt_model(k) = sumsmelt_model(k) / swe_count
endif
enddo
endif
c To avoid a divide by zero later on, make sure sumsprec_model and
c sumsmelt_model are not both zero.
do k=1,nstns
sumsprec_model(k) = sumsprec_model(k) + 1.0e-6
enddo
c Determine whether we will adjust the precipitation or melt. To
c do this, calculate the relative contributions of precipitation
c and melt inputs for this correction period. This can be
c different for each observation interval. Calculate the average
c over all of the stations/areas in the domain.
sprec_ratio(iobs_num) = 0.0
smelt_ratio(iobs_num) = 0.0
do k=1,nstns
sprec_ratio(iobs_num) = sprec_ratio(iobs_num) +
& sumsprec_model(k) / (sumsprec_model(k)+sumsmelt_model(k))
smelt_ratio(iobs_num) = smelt_ratio(iobs_num) +
& sumsmelt_model(k) / (sumsprec_model(k)+sumsmelt_model(k))
enddo
sprec_ratio(iobs_num) = sprec_ratio(iobs_num) / real(nstns)
smelt_ratio(iobs_num) = smelt_ratio(iobs_num) / real(nstns)
c Initialize the delta swe variable.
if (iobs_num.eq.1) then
do k=1,nstns
delta_old(k) = 0.0
enddo
else
do k=1,nstns
delta_old(k) = 0.0
enddo
do k=1,nstns
do kk=1,nstns_old
if(obsid(k).eq.obsid_old(kk))
& delta_old(k) = delta_old_tmp(kk)
enddo
enddo
c write (77,*)
c do k=1,nstns
c write (77,*) 'k, delta_old(k)',k,100.*delta_old(k)
c enddo
c write (77,*)
endif
c Calculate the correction factor to be used in the next model
c iteration. Let the correction factor equal 1.0 during
c periods where we have no swe observations. Also, note that the
c reason for the delta_old variable is to account for the fact
c that that delta will be fixed with the previous date correction
c time period. This is one of the things that allows the
c correction to be done in two model iterations.
c If sumsprec_model or sumsmelt_model are too small to be used in
c the assimilation (like less than 1 mm), set corr_factor_tmp = 1.0
c so no adjustments are performed for this observation interval.
cf_min = 0.1
do k=1,nstns
if (depth_assim.eq.1) then
if (ro_model(k).le.0) then
swe_obs(k) = ro_avg*depth_obs(k)/1000.0
else
swe_obs(k) = ro_model(k)*depth_obs(k)/1000.0
endif
endif
if (sprec_ratio(iobs_num).ge.smelt_ratio(iobs_num)) then
if (sumsprec_model(k).lt.1.0e-3) then
corr_factor_tmp(k) = 1.0
else
corr_factor_tmp(k) = 1.0 +
& (swe_obs(k) - swe_model(k) - delta_old(k)) /
& sumsprec_model(k)
corr_factor_tmp(k) = max(cf_min,corr_factor_tmp(k))
endif
else
if (sumsmelt_model(k).lt.1.0e-3) then
corr_factor_tmp(k) = 1.0
else
corr_factor_tmp(k) = 1.0 +
& (swe_model(k) - swe_obs(k) + delta_old(k)) /
& sumsmelt_model(k)
corr_factor_tmp(k) = max(cf_min,corr_factor_tmp(k))
endif
endif
c Save some information about the model calculations.
c write (77,*) '---'
c write (77,*) k,swe_obs(k)
c write (77,*) k,swe_model(k)
c write (77,*) k,delta_old(k)
c write (77,*) k,swe_obs(k)-swe_model(k)-delta_old(k)
c write (77,*) k,sumsprec_model(k)
c write (77,*) k,sumsmelt_model(k)
c write (77,*) k,corr_factor_tmp(k)
c write (77,*) '---'
c Save some data from this observation time for use at the next
c observation time.
nstns_old = nstns
obsid_old(k) = obsid(k)
delta_old_tmp(k) = swe_obs(k) - swe_model(k)
enddo
c Now that I have the correction factors calculated at each
c observation point, interpolate those over the simulation domain.
c Use the barnes oi scheme to create the distribution. If there is
c only a single station, distribute those data uniformly over
c the domain.
if (nstns.ge.2) then
call get_dn(nx,ny,deltax,deltay,nstns,dn,iobsint)
c Modify the size of dn.
dn = beta * dn
call barnes_oi(nx,ny,deltax,deltay,xmn,ymn,
& nstns,xstn,ystn,corr_factor_tmp,dn,grid,undef,ifill)
elseif (nstns.eq.1) then
call single_stn(nx,ny,nstns,corr_factor_tmp,grid)
endif
c The following calculations are done if you want to implement the
c special case where you limit the data assimilation corrections
c to a certain area of your simulation domain. Edits to this
c subroutine will certainly be required to make it work for your
c specific application. This subroutine generates correction
c factors with 1.0's in places outside the local obs influences.
if (local_assim_flag.eq.1) then
call mk_local_cfs(nx,ny,undef,xmn,ymn,deltax,deltay,
& fname_sweobs,fname_sweobs_barnes_mask,nobs_dates,
& corr_factor_tmp,beta,iobsint,ifill)
endif
c Define the correction surface record that corresponds to this
c year and observation.
krec = nobs_total + (nyear - 1)
if (krec.gt.max_obs_dates) then
print *, 'max_obs_dates must be increased in snowmodel.inc'
print *, 'krec = ',krec,' max_obs_dates = ',max_obs_dates
stop
endif
c Use the gridded output file to build the corr_factor array.
do j=1,ny
do i=1,nx
corr_factor(i,j,krec) = grid(i,j)
corr_factor(i,j,krec) =
& max(cf_min,corr_factor(i,j,krec))
enddo
enddo
c Note that the interpolation scheme may have produced correction
c factors that do not produce exact matches with the
c observations (like happens with the case of having a single
c data point). If you are interested, calculate the difference
c between the exact value and the actual calculated value, and
c then write it out as done below.
c do k=1,nstns
c if (sprec_ratio(iobs_num).ge.smelt_ratio(iobs_num)) then
c corr_offset(k) = sumsprec_model(k) *
c & (corr_factor(ii(k),jj(k),iobs_num) - corr_factor_tmp(k))
c else
c corr_offset(k) = sumsmelt_model(k) *
c & (corr_factor(ii(k),jj(k),iobs_num) - corr_factor_tmp(k))
c endif
c enddo
c Write some information to the text file.
write (77,*) '***************************************'
write (77,*) ' sprec_ratio =',sprec_ratio(iobs_num),
& ' smelt_ratio =',smelt_ratio(iobs_num)
write (77,*)
do k=1,nstns
write (77,*) k,' swe diff =',
& 100.0*abs(swe_obs(k)-swe_model(k)),' SWE OBS =',
& 100.0*swe_obs(k)
write (77,*) 'sumsprec =',sumsprec_model(k)*100.,
& ' SWE MODEL =',swe_model(k)*100.
write (77,*) 'iobs_num =',iobs_num,
& ' CORR_FACTOR =',corr_factor_tmp(k)
cc write (77,*) 'corr_offset =',100.*corr_offset(k),
cc & ' ij',ii(k),jj(k)
cc write (77,*) ' delta_old =',100.*delta_old(k),
cc & ' corr fact used =',corr_factor(ii(k),jj(k),iobs_num)
c write (77,*)
c write (77,*) k,' sumsprec_model(k) =',sumsprec_model(k)
c write (77,*) k,' sumsmelt_model(k) =',sumsmelt_model(k)
c write (77,*)
enddo
write (77,*) '***************************************'
c Write the output data to a grads file.
write(62,rec=krec) ((corr_factor(i,j,krec),i=1,nx),j=1,ny)
enddo
c Fill corr_factor with 1.0 for the period following the last obs
c date in the current year. This is also required for the history
c restart to work correctly. Without the history restart this was
c already done as part of the model initialization.
if (krec+1.gt.max_obs_dates) then
print *, 'max_obs_dates must be increased in snowmodel.inc'
print *, 'krec+1 = ',krec+1,' max_obs_dates = ',max_obs_dates
stop
endif
do j=1,ny
do i=1,nx
corr_factor(i,j,krec+1) = 1.0
enddo
enddo
write(62,rec=krec+1) ((corr_factor(i,j,krec+1),i=1,nx),j=1,ny)
c The met, topo, and veg files must be closed for the next model
c iteration.
close (20)
close (37)
close (38)
close (238)
close (239)
close (240)
close (237)
else
c For the case of no observations for this year, set the correction
c factor equal to 1.0.
krec = nobs_total + nyear
if (krec.gt.max_obs_dates) then
print *, 'max_obs_dates must be increased in snowmodel.inc'
print *, 'krec = ',krec,' max_obs_dates = ',max_obs_dates
stop
endif
do j=1,ny
do i=1,nx
corr_factor(i,j,krec) = 1.0
enddo
enddo
write(62,rec=krec) ((corr_factor(i,j,krec),i=1,nx),j=1,ny)
endif
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine corr_factor_index(nobs_dates,icorr_factor_index,
& iobs_rec,max_iter,sprec_ratio,smelt_ratio,print_inc,
& iday_init,imonth_init,iyear_init,nyear,nobs_total_cfi,
& nyears,ihrestart_flag)
implicit none
include 'snowmodel.inc'
integer icorr_factor_index(max_time_steps)
integer kk,istart,iend,nobs_dates,iter,max_iter,krec,nyear,
& nobs_total_cfi,ioptn,julian_start,iiyr,julian_end,iday_init,
& imonth_init,iyear_init,nyears,ihrestart_flag
integer iobs_rec(max_obs_dates)
real sprec_ratio(max_obs_dates),smelt_ratio(max_obs_dates)
real print_inc
if (ihrestart_flag.eq.0) print_inc = 24.0
c Initialize the number-of-observations counter.
if (nyear.eq.1) nobs_total_cfi = 0
c Build an array indicating the appropriate correction factor to
c use at each time step during the simulation.
if (nobs_dates.gt.0) then
c Loop through the observation dates.
do kk=1,nobs_dates+1
c Increment the number-of-observations counter.
if (kk.le.nobs_dates) nobs_total_cfi = nobs_total_cfi + 1
c FIRST, FROM THE SIMULATION START UNTIL THE FIRST OBSERVATION, FOR
c EACH YEAR.
if (kk.eq.1) then
c Here istart equals the first model time step of each year.
ioptn = 3
call calndr (ioptn,iday_init,imonth_init,iyear_init,
& julian_start)
c Find the Julian day for this data record.
iiyr = iyear_init + (nyear - 1)
call calndr (ioptn,iday_init,imonth_init,iiyr,julian_end)
c Calculate istart.
istart = 1 + (julian_end - julian_start) * nint(print_inc)
c Here iend equals the model time step corresponding to the end of
c the first observation day. Take advantage of the data-file
c record that was calculated before. Convert from daily data
c output records to model time steps using print_inc.
iend = iobs_rec(kk) * nint(print_inc)
c Define the corr_factor data array record.
krec = nobs_total_cfi + (nyear - 1)
c Fill the index for each model time step.
do iter=istart,iend
if (sprec_ratio(kk).ge.smelt_ratio(kk)) then
icorr_factor_index(iter) = krec
else
icorr_factor_index(iter) = -krec
endif
enddo
c SECOND, BETWEEN THE LAST OBSERVATION AND THE END OF THE SIMULATION.
elseif (kk.eq.nobs_dates+1) then
istart = iobs_rec(kk-1) * nint(print_inc) + 1
c Here iend equals the last time step in the year of interest.
c Find the Julian day for this data record.
iiyr = iyear_init + nyear
call calndr (ioptn,iday_init,imonth_init,iiyr,julian_end)
c Calculate iend for this year.
iend = (julian_end - julian_start) * nint(print_inc)
iend = min(iend,max_iter)
c Define the corr_factor data array record.
krec = nobs_total_cfi + nyear
c Fill the index for each model time step.
do iter=istart,iend
icorr_factor_index(iter) = krec
enddo
c THIRD, ANY PERIODS BETWEEN OBSERVATIONS.
else
istart = iobs_rec(kk-1) * nint(print_inc) + 1
iend = iobs_rec(kk) * nint(print_inc)
c Define the corr_factor data array record.
krec = nobs_total_cfi + (nyear - 1)
c Fill the index for each model time step.
do iter=istart,iend
if (sprec_ratio(kk).ge.smelt_ratio(kk)) then
icorr_factor_index(iter) = krec
else
icorr_factor_index(iter) = -krec
endif
enddo
endif
enddo
else
c Create an array indes for the case of no observations for this
c year. Here istart equals the first model time step of each year.
ioptn = 3
call calndr (ioptn,iday_init,imonth_init,iyear_init,
& julian_start)
c Find the Julian day for this data record.
iiyr = iyear_init + (nyear - 1)
call calndr (ioptn,iday_init,imonth_init,iiyr,julian_end)
c Calculate istart.
istart = 1 + (julian_end - julian_start) * nint(print_inc)
c Calculate iend for this year.
iiyr = iyear_init + nyear
call calndr (ioptn,iday_init,imonth_init,iiyr,julian_end)
iend = (julian_end - julian_start) * nint(print_inc)
iend = min(iend,max_iter)
c Define the corr_factor data array record.
krec = nobs_total_cfi + nyear
c Fill the index for each model time step.
do iter=istart,iend
icorr_factor_index(iter) = krec
enddo
endif
c SAVE A VERSION OF THE INDEX THAT THE USER CAN EASILY LOOK AT.
if (nyear.eq.nyears) then
print *
do iter=1,max_iter
write (77,*) iter,icorr_factor_index(iter)
enddo
print *
endif
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine mk_local_cfs(nx,ny,undef,xmn,ymn,deltax,deltay,
& fname_sweobs,fname_sweobs_barnes_mask,nobs_dates,
& corr_factor_tmp,beta,iobsint,ifill)
implicit none
include 'snowmodel.inc'
integer nstns,k,nx,ny,i,j,ntstations,icount,nobs_dates,ifill,
& iobsint,nstns2
real dummy,stnid,undef,deltax,deltay,beta,dn
real xmask(nx_max,ny_max),grid(nx_max,ny_max)
real corr_factor_tmp(nx_max*ny_max),cf_tmp2(nx_max*ny_max),
& obsid(nx_max*ny_max)
real var(nstns_max)
double precision xmn,ymn
double precision yg(nx_max,ny_max),xg(nx_max,ny_max)
double precision xstn(nstns_max),ystn(nstns_max)
character*80 fname_sweobs
character*80 fname_sweobs_barnes_mask
print *
print *,'You are doing local assimilation, this requires'
print *,'an observation mask to have been generated before'
print *,'the model run and saved in the file called:'
print *,'fname_sweobs_barnes_mask'
print *
print *,'This was also not tested after some big changes to'
print *,'the data assimilation code. So I strongly suggest'
print *,'you make sure this is doing what you want when you'
print *,'first start using it.'
print *
stop
if (nobs_dates.gt.1) then
print *,'THIS HAS NOT BEEN MADE TO WORK WITH MORE THAN'
print *,'ONE OBS TIME.'
stop
endif
print *
c Save the correction factors for the local-influence assimilation
c scheme.
open (unit=78,file='data/cf.txt')
do k=1,nstns
write (78,*) corr_factor_tmp(k)
enddo
close (78)
c These are the obs and correction factors calculated from the
c first loop in SnowModel.
open (721,file=fname_sweobs)
open (731,file='data/cf.txt')
read (721,*) nstns
do k=1,nstns
read (721,*) stnid,xstn(k),ystn(k),dummy
read (731,*) var(k)
c if (var(k).gt.1.5) then
c print *, 'cf > 1.5 found; setting to 1.5',k,var(k)
c var(k) = 1.5
c endif
c if (var(k).lt.0.5) then
c print *, 'cf < 0.5 found; setting to 0.5',k,var(k)
c var(k) = 0.5
c endif
c var(k) = min(1.5,var(k))
c var(k) = max(0.5,var(k))
enddo
close (721)
close (731)
c Create a collection of 'stations' with correction factors of
c 1.0 in areas outside of our traverse regions.
open(741,file=fname_sweobs_barnes_mask,