-
Notifications
You must be signed in to change notification settings - Fork 4
/
micromet_code.f
3468 lines (2793 loc) · 118 KB
/
micromet_code.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 micromet_code.f
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine MICROMET_CODE(nx,ny,xmn,ymn,deltax,deltay,
& iyear_init,imonth_init,iday_init,xhour_init,dt,undef,
& ifill,iobsint,dn,iter,curve_len_scale,slopewt,curvewt,
& topo,curvature,terrain_slope,slope_az,Tair_grid,
& rh_grid,uwind_grid,vwind_grid,Qsi_grid,prec_grid,
& i_tair_flag,i_rh_flag,i_wind_flag,i_solar_flag,
& i_prec_flag,isingle_stn_flag,igrads_metfile,
& windspd_grid,winddir_grid,windspd_flag,winddir_flag,
& sprec,windspd_min,Qli_grid,i_longwave_flag,vegtype,
& forest_LAI,iyear,imonth,iday,xhour,corr_factor,
& icorr_factor_index,lapse_rate_user_flag,
& iprecip_lapse_rate_user_flag,use_shortwave_obs,
& use_longwave_obs,use_sfc_pressure_obs,sfc_pressure,
& run_enbal,run_snowpack,calc_subcanopy_met,vegsnowd_xy,
& gap_frac,cloud_frac_factor,barnes_lg_domain,n_stns_used,
& k_stn,xlat_grid,xlon_grid,UTC_flag,icorr_factor_loop,
& snowmodel_line_flag,xg_line,yg_line,irun_data_assim,
& wind_lapse_rate,prec_grid_sol,pertPrec)
implicit none
include 'snowmodel.inc'
integer nx ! number of x output values
integer ny ! number of y output values
real deltax ! grid increment in x
real deltay ! grid increment in y
double precision xmn ! center x coords of lower left grid cell
double precision ymn ! center y coords of lower left grid cell
integer nstns_orig ! number of input values
double precision xg_line(nx_max,ny_max),yg_line(nx_max,ny_max)
real snowmodel_line_flag
double precision xstn_orig(nstns_max) ! input stn x coords
double precision ystn_orig(nstns_max) ! input stn y coords
real Tair_orig(nstns_max) ! input values
real rh_orig(nstns_max) ! input values
real winddir_orig(nstns_max) ! input values
real windspd_orig(nstns_max) ! input values
real prec_orig(nstns_max) ! input values
real prec_orig_sol(nstns_max)
real elev_orig(nstns_max) ! station elevation
real dn ! average observation spacing
real topo(nx_max,ny_max) ! grid topography
real xlat_grid(nx_max,ny_max) ! lat (dec deg) of cell centers
real xlon_grid(nx_max,ny_max) ! lon (dec deg) of cell centers
real Tair_grid(nx_max,ny_max) ! output values
real rh_grid(nx_max,ny_max) ! output values
real uwind_grid(nx_max,ny_max) ! output, E-W wind component
real vwind_grid(nx_max,ny_max) ! output, N-S wind component
real windspd_grid(nx_max,ny_max)
real winddir_grid(nx_max,ny_max)
real Qsi_grid(nx_max,ny_max) ! output
real Qli_grid(nx_max,ny_max) ! output
real prec_grid(nx_max,ny_max) ! output
real prec_grid_sol(nx_max,ny_max)
real sprec(nx_max,ny_max) ! output
real sfc_pressure(nx_max,ny_max)
integer iyear,imonth,iday ! model year, month, and day
real xhour ! model decimal hour
real dt ! model time step, in seconds
integer iter ! model iteration
integer iyear_init ! model start year
integer imonth_init ! model start month
integer iday_init ! model start day
real xhour_init ! model start hour
integer J_day ! model Julian day, actually day-of-year
real undef ! undefined value
integer ifill ! flag (=1) forces a value in every cell
integer iobsint ! flag (=1) use dn value from .par file
real curvature(nx_max,ny_max) ! topographic curvature
real slope_az(nx_max,ny_max) ! azimuth of topographic slope
real terrain_slope(nx_max,ny_max) ! terrain slope
real vegtype(nx_max,ny_max)
real vegsnowd_xy(nx_max,ny_max)
real topo_ref_grid(nx_max,ny_max) ! reference surface
real curve_len_scale ! length scale for curvature calculation
real slopewt ! wind model slope weight
real curvewt ! wind model curvature weight
integer i_tair_flag,i_rh_flag,i_wind_flag,i_solar_flag,
& i_prec_flag,i_longwave_flag,isingle_stn_flag,igrads_metfile,
& lapse_rate_user_flag,iprecip_lapse_rate_user_flag,n_stns_used,
& icorr_factor_loop,irun_data_assim
real windspd_flag,winddir_flag,windspd_min,calc_subcanopy_met,
& T_lapse_rate,Td_lapse_rate,precip_lapse_rate,
& use_shortwave_obs,use_longwave_obs,use_sfc_pressure_obs,
& run_enbal,run_snowpack,gap_frac,cloud_frac_factor,
& barnes_lg_domain,UTC_flag,wind_lapse_rate,pertPrec
integer nftypes
parameter (nftypes=5)
real forest_LAI(nftypes)
real corr_factor(nx_max,ny_max,max_obs_dates+1)
integer icorr_factor_index(max_time_steps)
integer k_stn(nx_max,ny_max,5)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Calculate what the current simulation date should be.
call get_model_time(iyear_init,imonth_init,iday_init,
& xhour_init,iter,dt,iyear,imonth,iday,xhour,J_day)
if (irun_data_assim.eq.1) then
write (*,150) icorr_factor_loop,iyear,imonth,iday,xhour
150 format('In Assim Loop #',i1,'; WORKING ON MODEL TIME =',
& i5,2i4,f6.1)
else
write (*,151) iyear,imonth,iday,xhour
151 format(' WORKING ON MODEL TIME =',
& i5,2i4,f6.1)
endif
c Read in the observations for this time step, and build an array of
c valid observations to be interpolated.
call get_obs_data(nstns_orig,Tair_orig,rh_orig,xstn_orig,
& ystn_orig,elev_orig,iyear,imonth,iday,xhour,undef,
& windspd_orig,winddir_orig,prec_orig,isingle_stn_flag,
& igrads_metfile,iter,pertPrec,irun_data_assim,
& icorr_factor_loop)
c Make the topographic calculations required by the wind and solar
c radiation models. These calculations are not fixed in time
c because as the snow depth evolves it modifies the "topography".
call topo_data(nx,ny,deltax,deltay,topo,
& curvature,terrain_slope,slope_az,curve_len_scale)
c Calculate the temperature and dew-point lapse rates to be used in
c the interpoations.
call get_lapse_rates(imonth,iday,T_lapse_rate,
& Td_lapse_rate,xlat_grid(1,1),lapse_rate_user_flag,
& precip_lapse_rate,iprecip_lapse_rate_user_flag)
c Calculate the forest lai for each of the five forest types, and
c for this day of the simulation (the lai varies seasonally for
c the case of deciduous trees).
call get_lai(J_day,forest_LAI)
c TEMPERATURE.
if (i_tair_flag.eq.1) then
c print *,' solving for temperature'
call temperature(nx,ny,deltax,deltay,xmn,ymn,
& nstns_orig,xstn_orig,ystn_orig,Tair_orig,dn,Tair_grid,
& undef,ifill,iobsint,iyear,imonth,iday,xhour,elev_orig,
& topo,T_lapse_rate,barnes_lg_domain,n_stns_used,k_stn,
& snowmodel_line_flag,xg_line,yg_line)
endif
c RELATIVE HUMIDITY.
if (i_rh_flag.eq.1) then
c print *,' solving for relative humidity'
call relative_humidity(nx,ny,deltax,deltay,xmn,ymn,
& nstns_orig,xstn_orig,ystn_orig,rh_orig,dn,rh_grid,undef,
& ifill,iobsint,iyear,imonth,iday,xhour,elev_orig,topo,
& Tair_orig,Tair_grid,Td_lapse_rate,barnes_lg_domain,
& n_stns_used,k_stn,snowmodel_line_flag,xg_line,yg_line)
endif
c WIND SPEED AND DIRECTION.
if (i_wind_flag.eq.1) then
c print *,' solving for wind speed and direction'
call wind(nx,ny,deltax,deltay,xmn,ymn,windspd_orig,
& nstns_orig,xstn_orig,ystn_orig,dn,undef,ifill,
& iobsint,iyear,imonth,iday,xhour,elev_orig,
& winddir_orig,uwind_grid,vwind_grid,slopewt,curvewt,
& curvature,slope_az,terrain_slope,windspd_grid,
& winddir_grid,windspd_flag,winddir_flag,windspd_min,
& vegtype,forest_LAI,calc_subcanopy_met,vegsnowd_xy,
& barnes_lg_domain,n_stns_used,k_stn,snowmodel_line_flag,
& xg_line,yg_line,topo_ref_grid,topo,wind_lapse_rate)
c Provide the ability to read in an alternate wind dataset that
c has been previously generated with another program, like
c NUATMOS. The following assumes you are reading in a single
c file with u and v values. The file name is hardcoded here.
elseif (i_wind_flag.eq.-1) then
call read_wind_file(nx,ny,iter,uwind_grid,vwind_grid,
& windspd_grid,winddir_grid,windspd_flag,winddir_flag,
& windspd_min)
endif
c SOLAR RADIATION.
if (i_solar_flag.eq.1) then
c print *,' solving for solar radiation'
call solar(nx,ny,xhour,J_day,topo,rh_grid,Tair_grid,
& xlat_grid,Qsi_grid,slope_az,terrain_slope,dt,vegtype,
& forest_LAI,T_lapse_rate,Td_lapse_rate,
& calc_subcanopy_met,gap_frac,cloud_frac_factor,UTC_flag,
& xlon_grid)
c If requested, modify the model output to account for shortwave
c radiation observations.
if (use_shortwave_obs.eq.1.0) then
if (barnes_lg_domain.eq.1.0) then
print *,'The model is not configured to assimilate'
print *,' solar data with barnes_lg_domain = 1.0.'
stop
endif
call shortwave_data(nx,ny,deltax,deltay,xmn,ymn,
& iyear,imonth,iday,xhour,undef,Qsi_grid,iter)
endif
endif
c INCOMING LONGWAVE RADIATION.
if (i_longwave_flag.eq.1) then
c print *,' solving for incoming longwave radiation'
call longwave(nx,ny,rh_grid,Tair_grid,Qli_grid,topo,
& vegtype,forest_LAI,T_lapse_rate,Td_lapse_rate,
& calc_subcanopy_met,cloud_frac_factor)
c If requested, modify the model output to account for longwave
c radiation observations.
if (use_longwave_obs.eq.1.0) then
if (barnes_lg_domain.eq.1.0) then
print *,'The model is not configured to assimilate'
print *,' longwave data with barnes_lg_domain = 1.0.'
stop
endif
call longwave_data(nx,ny,deltax,deltay,xmn,ymn,
& iyear,imonth,iday,xhour,undef,Qli_grid,iter)
endif
endif
c PRECIPITATION.
if (i_prec_flag.eq.1) then
c print *,' solving for precipitation'
call precipitation(nx,ny,deltax,deltay,xmn,ymn,
& nstns_orig,xstn_orig,ystn_orig,prec_orig,dn,prec_grid,
& undef,ifill,iobsint,iyear,imonth,iday,xhour,elev_orig,
& topo,Tair_grid,sprec,corr_factor,icorr_factor_index,iter,
& precip_lapse_rate,barnes_lg_domain,n_stns_used,k_stn,
& snowmodel_line_flag,xg_line,yg_line,topo_ref_grid)
c J.PFLUG
elseif (i_prec_flag.eq.-1) then
call read_frozen(nstns_orig,undef,prec_orig_sol,
& isingle_stn_flag,iter,pertPrec)
call precipitation_froz(nx,ny,deltax,deltay,xmn,ymn,
& nstns_orig,xstn_orig,ystn_orig,prec_orig,prec_orig_sol,
& dn,prec_grid,undef,ifill,iobsint,iyear,
& imonth,iday,xhour,elev_orig,topo,Tair_grid,sprec,
& corr_factor,icorr_factor_index,iter,precip_lapse_rate,
& barnes_lg_domain,n_stns_used,k_stn,snowmodel_line_flag,
& xg_line,yg_line,topo_ref_grid,prec_grid_sol)
endif
c END J.PFLUG
c SURFACE PRESSURE.
c Surface pressure is used in EnBal and SnowMass. If needed for
c this SnowModel simulation, calculate the distribution here.
if (run_enbal.eq.1.0 .or. run_snowpack.eq.1.0) then
call pressure(nx,ny,topo,sfc_pressure)
c If requested, modify the model output to account for surface
c pressure observations.
if (use_sfc_pressure_obs.eq.1.0) then
if (barnes_lg_domain.eq.1.0) then
print *,'The model is not configured to assimilate'
print *,' pressure data with barnes_lg_domain = 1.0.'
stop
endif
call sfc_pressure_data(nx,ny,deltax,deltay,xmn,ymn,
& iyear,imonth,iday,xhour,undef,sfc_pressure,iter)
endif
endif
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine precipitation(nx,ny,deltax,deltay,xmn,ymn,
& nstns_orig,xstn_orig,ystn_orig,prec_orig,dn,prec_grid,
& undef,ifill,iobsint,iyear,imonth,iday,xhour,elev_orig,
& topo,Tair_grid,sprec,corr_factor,icorr_factor_index,iter,
& precip_lapse_rate,barnes_lg_domain,n_stns_used,k_stn,
& snowmodel_line_flag,xg_line,yg_line,topo_ref_grid)
c Interpolate the observed precipitation values to the grid. Also
c interpolate the station elevations to a reference surface. Use
c a precipitation "lapse rate", or adjustment factor to define
c the precipitation on the actual elevation grid. The reason the
c interpolated station elevations are used as the topographic
c reference surface (instead of something like sea level), is
c because the precipitation adjustment factor is a non-linear
c function of elevation difference.
c The adjustment factor that is used comes from: Thornton, P. E.,
c S. W. Running, and M. A. White, 1997: Generating surfaces of
c daily meteorological variables over large regions of complex
c terrain. J. Hydrology, 190, 214-251.
implicit none
include 'snowmodel.inc'
integer nx ! number of x output values
integer ny ! number of y output values
real deltax ! grid increment in x
real deltay ! grid increment in y
double precision xmn ! center x coords of lower left grid cell
double precision ymn ! center y coords of lower left grid cell
double precision xg_line(nx_max,ny_max),yg_line(nx_max,ny_max)
real snowmodel_line_flag
integer nstns ! number of input values, all good
integer nstns_orig ! number of input values
double precision xstn(nstns_max) ! input stn x coords
double precision ystn(nstns_max) ! input stn y coords
real prec(nstns_max) ! input values
real elev(nstns_max) ! station elevation
real undef ! undefined value
double precision xstn_orig(nstns_max) ! input stn x coords
double precision ystn_orig(nstns_max) ! input stn y coords
real elev_orig(nstns_max) ! station elevation
real prec_orig(nstns_max) ! input values
real dn ! average observation spacing
real topo(nx_max,ny_max) ! grid topography
real prec_grid(nx_max,ny_max) ! output values
real Tair_grid(nx_max,ny_max) ! input values
real sprec(nx_max,ny_max) ! output values
real topo_ref_grid(nx_max,ny_max) ! reference surface
integer ifill ! flag (=1) forces a value in every cell
integer iobsint ! flag (=1) use dn value from .par file
integer iyear,imonth,iday ! model year, month, and day
real xhour ! model decimal hour
real delta_topo,alfa,Tf,precip_lapse_rate_m,precip_lapse_rate,
& barnes_lg_domain
integer i,j,iter,n_stns_used
integer k_stn(nx_max,ny_max,5)
real corr_factor(nx_max,ny_max,max_obs_dates+1)
integer icorr_factor_index(max_time_steps)
c Filter through the original input data, and eliminate any
c missing values.
call get_good_values1(nstns_orig,xstn_orig,ystn_orig,
& elev_orig,undef,nstns,xstn,ystn,elev,prec_orig,prec)
c Use the barnes oi scheme to interpolate the station elevation data
c to the grid, so that it can be used as a topographic reference
c surface.
call interpolate(nx,ny,deltax,deltay,xmn,ymn,
& nstns,xstn,ystn,elev,dn,topo_ref_grid,undef,ifill,iobsint,
& iyear,imonth,iday,xhour,barnes_lg_domain,n_stns_used,
& k_stn,snowmodel_line_flag,xg_line,yg_line)
c Use the barnes oi scheme to interpolate the station data to
c the grid.
call interpolate(nx,ny,deltax,deltay,xmn,ymn,
& nstns,xstn,ystn,prec,dn,prec_grid,undef,ifill,iobsint,
& iyear,imonth,iday,xhour,barnes_lg_domain,n_stns_used,
& k_stn,snowmodel_line_flag,xg_line,yg_line)
c Convert the precipitation "lapse rate" (km-1) to m-1.
precip_lapse_rate_m = precip_lapse_rate / 1000.0
c Convert the gridded station data to the actual gridded elevations.
do j=1,ny
do i=1,nx
delta_topo = topo(i,j) - topo_ref_grid(i,j)
c Don't let the elevation difference be greater than some number
c (like 1800 meters gives a factor of 4.4). If it is too large
c you get huge precipitation adjustment, a divide by zero, or
c even negative adjustments for high elevations).
delta_topo = min(delta_topo,1800.0)
alfa = precip_lapse_rate_m * delta_topo
prec_grid(i,j) = prec_grid(i,j) * (1.0 + alfa)/(1.0 - alfa)
enddo
enddo
c Convert the precipitation values from mm to m swe. Also, make
c sure the interpolation has not created any negetive
c precipitation values.
do j=1,ny
do i=1,nx
prec_grid(i,j) = prec_grid(i,j) / 1000.0
prec_grid(i,j) = max(0.0,prec_grid(i,j))
enddo
enddo
c Use the temperature distribution to define whether this
c precipitation is falling as rain or snow (generate a
c snow-precipitation array following the air temperature
c threshold defined by Auer (1974) = 2.0 C). Note here that,
c if you ever want it, rain_prec = prec_grid - sprec. This
c snow precipitation (sprec) is assumed to be in meters
c snow-water-equivalent per time step.
Tf = 273.16
do j=1,ny
do i=1,nx
if (Tair_grid(i,j).le.-0.5+Tf) then
if (icorr_factor_index(iter).gt.0.0) then
prec_grid(i,j) =
& corr_factor(i,j,icorr_factor_index(iter)) *
& prec_grid(i,j)
else
prec_grid(i,j) = prec_grid(i,j)
endif
sprec(i,j) = prec_grid(i,j)
else
sprec(i,j) = 0.0
endif
enddo
enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine longwave(nx,ny,rh_grid,Tair_grid,Qli_grid,topo,
& vegtype,forest_LAI,T_lapse_rate,Td_lapse_rate,
& calc_subcanopy_met,cloud_frac_factor)
implicit none
include 'snowmodel.inc'
integer nx ! number of x output values
integer ny ! number of y output values
real Tair_grid(nx_max,ny_max)
real rh_grid(nx_max,ny_max)
real Qli_grid(nx_max,ny_max)
real topo(nx_max,ny_max)
real vegtype(nx_max,ny_max)
real T_lapse_rate,Td_lapse_rate,es,e,emiss_cloud,
& A,B,C,Tf,Stef_Boltz,cloud_frac,E1,X1,Y1,Z1,E2,X2,Y2,Z2,
& Xs,Ys,Zs,forest_frac,E3,X3,Y3,Z3,alfa,calc_subcanopy_met,
& cloud_frac_factor
integer nftypes
parameter (nftypes=5)
real forest_LAI(nftypes)
integer i,j,nveg
c Coeffs for saturation vapor pressure over water (Buck 1981).
c Note: temperatures for Buck`s equations are in deg C, and
c vapor pressures are in mb. Do the adjustments so that the
c calculations are done with temperatures in K, and vapor
c pressures in Pa.
A = 6.1121 * 100.0
B = 17.502
C = 240.97
c Over ice.
c A = 6.1115 * 100.0
c B = 22.452
c C = 272.55
c Define the freezing temperature to be used to convert from C to K.
Tf = 273.16
c Define the Stefan Boltzmann constant.
Stef_Boltz = 5.6696e-8
c Constants required for Iziomon et al. (2003).
E1 = 200.0
X1 = 0.35
Y1 = 0.100
Z1 = 0.224
E2 = 1500.0
X2 = 0.43
Y2 = 0.115
Z2 = 0.320
c Assume the X and Y coefficients increase linearly to 3000 m,
c and adjust Z to create a best fit to the CLPX data.
E3 = 3000.0
X3 = 0.51
Y3 = 0.130
Z3 = 1.100
do j=1,ny
do i=1,nx
c Compute the cloud fraction.
call get_cloudfrac(Tair_grid(i,j),rh_grid(i,j),topo(i,j),
& cloud_frac,T_lapse_rate,Td_lapse_rate,cloud_frac_factor)
c Calculate the vapor pressure.
es = A * exp((B * (Tair_grid(i,j) - Tf))/
& (C + (Tair_grid(i,j) - Tf)))
e = es * rh_grid(i,j) / 100.0
c Compute Qli following Iziomon et al. (2003).
if (topo(i,j).lt.E1) then
Xs = X1
Ys = Y1
Zs = Z1
elseif (topo(i,j).gt.E2) then
Xs = X3
Ys = Y3
Zs = Z3
else
Xs = X1 + (topo(i,j) - E1) * (X3 - X1)/(E3 - E1)
Ys = Y1 + (topo(i,j) - E1) * (Y3 - Y1)/(E3 - E1)
Zs = Z1 + (topo(i,j) - E1) * (Z3 - Z1)/(E3 - E1)
endif
alfa = 1.083
emiss_cloud = alfa *
& (1.0 - Xs * exp((- Ys) * e/Tair_grid(i,j))) *
& (1.0 + Zs * cloud_frac**2)
emiss_cloud = min(1.0,emiss_cloud)
Qli_grid(i,j) = emiss_cloud * Stef_Boltz * Tair_grid(i,j)**4
c Modify the incoming longwave radiation for the forest canopy.
if (vegtype(i,j).le.5.0) then
if (calc_subcanopy_met.eq.1.0) then
c Define the forest-canopy parameters.
nveg = nint(vegtype(i,j))
if (forest_LAI(nveg).lt.0.2) then
forest_frac = 0.5 * forest_LAI(nveg)
else
forest_frac =
& min(1.0,0.55 + 0.29 * log(forest_LAI(nveg)))
endif
Qli_grid(i,j) = Qli_grid(i,j) * (1.0 - forest_frac) +
& (Stef_Boltz * Tair_grid(i,j)**4) * forest_frac
endif
endif
enddo
enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine solar(nx,ny,xhour,J_day,topo,rh_grid,Tair_grid,
& xlat_grid,Qsi_grid,slope_az,terrain_slope,dt,vegtype,
& forest_LAI,T_lapse_rate,Td_lapse_rate,
& calc_subcanopy_met,gap_frac,cloud_frac_factor,UTC_flag,
& xlon_grid)
c First take the surface gridded fields of Tair and RH, and
c calculate Td for the topographic surface. Then use Tair and
c Td, and the associated lapse rates, and calculate Tair and Td
c for the 700 mb level. Use these surfaces to calculate RH at
c 700 mb, and convert these values to cloud fraction following
c Walcek, C. J., 1994: Cloud cover and its relationship to
c relative humidity during a spring midlatitude cyclone. Mon.
c Wea. Rev., 122, 1021-1035.
implicit none
include 'snowmodel.inc'
integer nx ! number of x output values
integer ny ! number of y output values
real topo(nx_max,ny_max)
real Tair_grid(nx_max,ny_max)
real rh_grid(nx_max,ny_max)
real Qsi_grid(nx_max,ny_max)
real slope_az(nx_max,ny_max)
real terrain_slope(nx_max,ny_max)
real vegtype(nx_max,ny_max)
real xlat_grid(nx_max,ny_max)
real xlon_grid(nx_max,ny_max)
real xhour ! model decimal hour
real xxhour
integer J_day ! model day of year
integer i,j,ihrs_day,ihour,nveg
real dt,cloud_frac,Qsi_tmp,Qsi_sum,trans_veg,
& T_lapse_rate,Td_lapse_rate,calc_subcanopy_met,gap_frac,
& cloud_frac_factor,UTC_flag
integer nftypes
parameter (nftypes=5)
real forest_LAI(nftypes)
ihrs_day = 24
do j=1,ny
do i=1,nx
c Compute the cloud fraction.
call get_cloudfrac(Tair_grid(i,j),rh_grid(i,j),topo(i,j),
& cloud_frac,T_lapse_rate,Td_lapse_rate,cloud_frac_factor)
c Compute the incoming solar radiation. The solar_rad subroutine
c calculates the instantaneous incoming solar radiation, so
c if the time step is very long, account for this by calculating
c the incoming solar radiation every 3 hours and then taking the
c average.
if (dt.le.10800.0) then
call solar_rad(Qsi_grid(i,j),J_day,xlat_grid(i,j),
& cloud_frac,xhour,slope_az(i,j),terrain_slope(i,j),
& UTC_flag,xlon_grid(i,j))
elseif (dt.eq.86400.0) then
Qsi_sum = 0.0
do ihour=3,ihrs_day,3
xxhour = real(ihour)
call solar_rad(Qsi_tmp,J_day,xlat_grid(i,j),
& cloud_frac,xxhour,slope_az(i,j),terrain_slope(i,j),
& UTC_flag,xlon_grid(i,j))
Qsi_sum = Qsi_sum + Qsi_tmp
enddo
Qsi_grid(i,j) = Qsi_sum / (real(ihrs_day)/3.0)
else
print *,'The model may not do what you want with this dt'
stop
endif
c Modify the incoming solar radiation for the forest canopy.
if (vegtype(i,j).le.5.0) then
if (calc_subcanopy_met.eq.1.0) then
c Define the forest-canopy transmissivity. 0.71 provided a
c best-fit to the observations, when averaged over the two years
c of hourly data.
nveg = nint(vegtype(i,j))
trans_veg = exp((- 0.71) * forest_LAI(nveg))
c Account for any gaps in the forest canopy that will allow
c direct incoming solar radiation to reach the snow surface.
trans_veg = gap_frac * (1.0 - trans_veg) + trans_veg
Qsi_grid(i,j) = trans_veg * Qsi_grid(i,j)
endif
endif
enddo
enddo
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine get_cloudfrac(Tair_grid,rh_grid,topo,
& cloud_frac,T_lapse_rate,Td_lapse_rate,cloud_frac_factor)
implicit none
real Td_lapse_rate,topo_ref,delta_topo,A,B,C,e,es,dx,
& T_lapse_rate,press_ratio,f_max,one_minus_RHe,f_1,
& Td_700,Tair_700,rh_700,cloud_frac,Tair_grid,rh_grid,
& topo,Td_grid,Tf,cloud_frac_factor
c Coeffs for saturation vapor pressure over water (Buck 1981).
c Note: temperatures for Buck`s equations are in deg C, and
c vapor pressures are in mb. Do the adjustments so that the
c calculations are done with temperatures in K, and vapor
c pressures in Pa.
A = 6.1121 * 100.0
B = 17.502
C = 240.97
c Over ice.
c A = 6.1115 * 100.0
c B = 22.452
c C = 272.55
c Define the freezing temperature to be used to convert from C to K.
Tf = 273.16
c Assume that 700 mb is equivalent to 3000 m in a standard
c atmosphere.
topo_ref = 3000.0
c Define the ratio of 700 mb level pressure to the surface pressure
c (~1000 mb).
press_ratio = 0.7
c Assume dx = 80.0 km, for Walcek (1994).
dx = 80.0
c Walcek coefficients.
f_max = 78.0 + 80.0/15.5
one_minus_RHe = 0.196 + (0.76-80.0/2834.0) * (1.0 - press_ratio)
f_1 = f_max * (press_ratio - 0.1) / 0.6 / 100.0
c Convert the gridded topo-surface RH to Td.
es = A * exp((B * (Tair_grid - Tf))/(C + (Tair_grid - Tf)))
e = es * max(10.0,rh_grid) / 100.0
Td_grid = C * log(e/A) / (B - log(e/A)) + Tf
c Convert the topo-surface temperature values to 700 mb values.
c delta_topo = topo - topo_ref
delta_topo = topo_ref - topo
Td_700 = Td_grid + Td_lapse_rate * delta_topo
Tair_700 = Tair_grid + T_lapse_rate * delta_topo
c Convert each Td to a gridded relative humidity (0-1).
e = A * exp((B * (Td_700 - Tf))/(C + (Td_700 - Tf)))
es = A * exp((B * (Tair_700 - Tf))/(C + (Tair_700 - Tf)))
rh_700 = e/es
rh_700 = min(1.0,rh_700)
rh_700 = max(0.0,rh_700)
c Use this RH at 700 mb to define the cloud fraction (0-1).
cloud_frac = f_1 * exp((rh_700 - 1.0)/one_minus_RHe)
c If the user wants to, reduce the calculate cloud fraction by the
c cloud_frac_factor.
cloud_frac = cloud_frac_factor * cloud_frac
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine solar_rad(Qsi,J_day,xlat,
& cloud_frac,xxhour,slope_az,terrain_slope,
& UTC_flag,xlon)
implicit none
integer J_day
real solar_const,days_yr,Trop_Can,solstice,pi,deg2rad,
& Qsi_direct,Qsi_diffuse,cos_i,cos_Z,Qsi,xlat,sin_z,xxhour,
& cloud_frac,slope_az,terrain_slope,sol_dec,hr_angl,
& trans_direct,trans_diffuse,Qsi_trans_dir,Qsi_trans_dif,
& sun_azimuth,slope_az_S0,xxxhour,UTC_flag,xlon
c Required constants.
solar_const = 1370.
days_yr = 365.25
Trop_Can = 0.41
solstice = 173.
pi = 2.0 * acos(0.0)
deg2rad = pi / 180.0
c COMPUTE THE BASIC SOLAR RADIATION PARAMETERS.
c Compute the solar declination angle (radians).
sol_dec = Trop_Can *
& cos(2.*pi * (real(J_day) - solstice)/days_yr)
c For the case of running UTC time and a latitudinal variation
c in solar radiation, adjust the time to correspond to the
c local time at this longitude position.
if (UTC_flag.ne.0.0) then
xxxhour = xxhour + xlon / 15.0
if (xxxhour.ge.24.0) xxxhour = xxxhour - 24.0
if (xxxhour.lt.0.0) xxxhour = xxxhour + 24.0
else
xxxhour = xxhour
endif
c Compute the sun's hour angle (radians).
hr_angl = (xxxhour * 15.0 - 180.0) * deg2rad
c Compute cos_Z. Note that the sin of the solar elevation angle,
c sin_alfa, is equal to the cosine of the solar zenith angle,
c cos_Z.
cos_Z = sin(sol_dec) * sin(xlat * deg2rad) +
& cos(sol_dec) * cos(xlat * deg2rad) * cos(hr_angl)
cos_Z = max(0.0,cos_Z)
c Account for clouds, water vapor, pollution, etc.
trans_direct = (0.6 + 0.2 * cos_Z) * (1.0 - cloud_frac)
trans_diffuse = (0.3 + 0.1 * cos_Z) * cloud_frac
c Compute the solar radiation transmitted through the atmosphere.
Qsi_trans_dir = solar_const * trans_direct
Qsi_trans_dif = solar_const * trans_diffuse
c COMPUTE THE CORRECTIONS TO ALLOW FOR TOPOGRAPHIC SLOPE AND ASPECT.
c The sine of the solar zenith angle.
sin_Z = sqrt(1.0 - cos_Z*cos_Z)
c Azimuth of the sun, with south having zero azimuth for the
c northern hemisphere.
sun_azimuth =
& asin(max(-1.0,min(1.0,cos(sol_dec)*sin(hr_angl)/sin_Z)))
if (xlat.lt.0.0) then
sun_azimuth = - sun_azimuth
endif
c Make the corrections so that the angles below the local horizon
c are still measured from the normal to the slope.
if (xlat.ge.0.0) then
if (hr_angl.lt.0.0) then
if (hr_angl.lt.sun_azimuth) sun_azimuth = - pi - sun_azimuth
elseif (hr_angl.gt.0.0) then
if (hr_angl.gt.sun_azimuth) sun_azimuth = pi - sun_azimuth
endif
else
c if (hr_angl.lt.0.0) then
c if (hr_angl.lt.sun_azimuth) sun_azimuth = - pi - sun_azimuth
c elseif (hr_angl.gt.0.0) then
c if (hr_angl.gt.sun_azimuth) sun_azimuth = pi - sun_azimuth
c endif
endif
c Build, from the variable with north having zero azimuth, a
c slope_azimuth value with south having zero azimuth. Also
c make north have zero azimuth if in the southern hemsisphere.
if (xlat.ge.0.0) then
if (slope_az.ge.180.0) then
slope_az_S0 = slope_az - 180.0
else
slope_az_S0 = slope_az + 180.0
endif
else
slope_az_S0 = slope_az
endif
c Compute the angle between the normal to the slope and the angle
c at which the direct solar radiation impinges on the sloping
c terrain (radians).
cos_i = cos(terrain_slope * deg2rad) * cos_Z +
& sin(terrain_slope * deg2rad) * sin_Z *
& cos(sun_azimuth - slope_az_S0 * deg2rad)
c Adjust the topographic correction due to local slope so that
c the correction is zero if the sun is below the local horizon
c (i.e., the slope is in the shade) or if the sun is below the
c global horizon.
if (cos_i.lt.0.0) cos_i = 0.0
if (cos_Z.le.0.0) cos_i = 0.0
c Adjust the solar radiation for slope, etc.
Qsi_direct = cos_i * Qsi_trans_dir
Qsi_diffuse = cos_Z * Qsi_trans_dif
c Combine the direct and diffuse solar components.
Qsi = Qsi_direct + Qsi_diffuse
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine wind(nx,ny,deltax,deltay,xmn,ymn,windspd_orig,
& nstns_orig,xstn_orig,ystn_orig,dn,undef,ifill,
& iobsint,iyear,imonth,iday,xhour,elev_orig,
& winddir_orig,uwind_grid,vwind_grid,slopewt,curvewt,
& curvature,slope_az,terrain_slope,windspd_grid,
& winddir_grid,windspd_flag,winddir_flag,windspd_min,
& vegtype,forest_LAI,calc_subcanopy_met,vegsnowd_xy,
& barnes_lg_domain,n_stns_used,k_stn,snowmodel_line_flag,
& xg_line,yg_line,topo_ref_grid,topo,wind_lapse_rate)
c This program takes the station wind speed and direction, converts
c them to u and v components, interpolates u and v to a grid,
c converts the gridded values to speed and direction, and then
c runs a simple wind model that adjusts those speeds and
c directions according to topographic slope and curvature
c relationships. The resulting speeds and directions are
c converted to u and v components and passed back to the main
c program to be written to a file. (All of the conversion
c between u-v and speed-dir is done because of the problems
c with interpolating over the 360/0 direction line.)
implicit none
include 'snowmodel.inc'
integer nx ! number of x output values
integer ny ! number of y output values
real deltax ! grid increment in x
real deltay ! grid increment in y
double precision xmn ! center x coords of lower left grid cell
double precision ymn ! center y coords of lower left grid cell
double precision xg_line(nx_max,ny_max),yg_line(nx_max,ny_max)
real snowmodel_line_flag
integer nstns ! number of input values, all good
integer nstns_orig ! number of input values
double precision xstn(nstns_max) ! input stn x coords
double precision ystn(nstns_max) ! input stn y coords
real elev(nstns_max) ! station elevation
real undef ! undefined value
double precision xstn_orig(nstns_max) ! input stn x coords
double precision ystn_orig(nstns_max) ! input stn y coords
real elev_orig(nstns_max) ! station elevation
real windspd_orig(nstns_max) ! input values
real winddir_orig(nstns_max) ! input values
real speed(nstns_max) ! input values
real dir(nstns_max) ! input values
real u(nstns_max) ! u component of wind
real v(nstns_max) ! v component of wind
real dn ! average observation spacing
real uwind_grid(nx_max,ny_max) ! output values
real vwind_grid(nx_max,ny_max) ! output values
real u_grid(nx_max,ny_max) ! temporary u wind component
real v_grid(nx_max,ny_max) ! temporary v wind component
real winddir_grid(nx_max,ny_max) ! temporary wind direction
real windspd_grid(nx_max,ny_max) ! temporary wind speed
real curvature(nx_max,ny_max) ! topographic curvature
real slope_az(nx_max,ny_max) ! azimuth of topographic slope
real terrain_slope(nx_max,ny_max) ! terrain slope
real vegtype(nx_max,ny_max)
real vegsnowd_xy(nx_max,ny_max)
real topo_ref_grid(nx_max,ny_max) ! reference surface
real topo(nx_max,ny_max)
integer ifill ! flag (=1) forces a value in every cell
integer iobsint ! flag (=1) use dn value from .par file
integer iyear,imonth,iday ! model year, month, and day
real xhour ! model decimal hour
real pi,deg2rad,rad2deg,slopewt,curvewt
integer i,j,k,n_stns_used
integer k_stn(nx_max,ny_max,5)
real windspd_flag,winddir_flag,u_sum,v_sum,windspd_min,
& calc_subcanopy_met,barnes_lg_domain,wind_lapse_rate,
& delta_topo,alfa1,alfa2
integer nftypes
parameter (nftypes=5)
real forest_LAI(nftypes)
c Define the required constants.
pi = 2.0 * acos(0.0)
deg2rad = pi / 180.0
rad2deg = 180.0 / pi
c Filter through the original input data, and eliminate any
c missing values (i.e., make sure each wind direction is paired
c up with a wind speed.
call get_good_values2(nstns_orig,xstn_orig,ystn_orig,
& elev_orig,undef,nstns,xstn,ystn,elev,windspd_orig,winddir_orig,
& speed,dir)
c Convert these station data to u and v wind components.
do k=1,nstns
speed(k) = max(windspd_min,speed(k))
u(k) = (- speed(k)) * sin(deg2rad * dir(k))
v(k) = (- speed(k)) * cos(deg2rad * dir(k))
enddo
c Use the barnes oi scheme to interpolate the station data to
c the grid.
c U component.
call interpolate(nx,ny,deltax,deltay,xmn,ymn,
& nstns,xstn,ystn,u,dn,u_grid,undef,ifill,iobsint,
& iyear,imonth,iday,xhour,barnes_lg_domain,n_stns_used,
& k_stn,snowmodel_line_flag,xg_line,yg_line)
c V component.
call interpolate(nx,ny,deltax,deltay,xmn,ymn,
& nstns,xstn,ystn,v,dn,v_grid,undef,ifill,iobsint,
& iyear,imonth,iday,xhour,barnes_lg_domain,n_stns_used,
& k_stn,snowmodel_line_flag,xg_line,yg_line)
c If desired, impose a wind speed increase with elevation. Here
c the wind_lapse_rate = the wind speed increase per 1-km elevation
c gain. The adjustment is patterned after the precipitation-
c elevation adjustment. topo_ref_grid here comes from the
c precipitation adjustment.
if (wind_lapse_rate.ne.0.0) then
alfa1 = (wind_lapse_rate - 1.0) / (1.0 + wind_lapse_rate)
c Convert to m-1.
alfa1 = alfa1 / 1000.0
do j=1,ny
do i=1,nx
delta_topo = topo(i,j) - topo_ref_grid(i,j)