-
Notifications
You must be signed in to change notification settings - Fork 2
/
ssib3.bak
5553 lines (5263 loc) · 272 KB
/
ssib3.bak
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=======================================================================
C THIS IS A SIMPLIFIED SIB MODEL
C=======================================================================
CM THIS CODE OF YONGKANG XUE'S HAS BEEN MODIFIED TO RUN AS A SUBROUTINE
CM INSIDE AN ENKF RADIOMETRIC DATA ASSIMILATION FOR ESTIMATION OF SWE.
CM THE VEGIN1, VEGIN2 AND CNTROLS SUBROUTINES NO LONGER READ DATA FROM
CM FILES, BUT EXTRACT IT FROM THESE ARRAYS; SOME VALUES THAT ARE
CM STATIC HAVE BEEN HARDCODED INTO THESE SUBROUTINES. THEREFORE, THE
CM run.ctl FILE WHICH USED TO CONTAIN THE NAMES OF ALL OF THESE FILES
CM IS NOW OBSELETE. FURTHERMORE, THE ALBEDO ADJUSTMENT BASED ON THE
CM xadj VARIABLE IS ALSO OBSELETE AND NOT IMPLEMENTED HERE. THE CODE
CM TO COMPUTE rnoffs AND rnoffb WAS ALSO REMOVED, AS IT WAS NOT BEING
CM USED IN THE PRESENT APPLICATION.
CM CONSTANT VEGETATION DATA IS CONTAINED IN THE VEGIN1ARG VARIABLE,
CM VEGETATION DATA WHICH CHANGES EACH MONTH IN THE VEGIN2ARG VARIABLE.
CM N_U - NUMBER OF FORCING DATA AT EACH TIMESTEP
CM N_STEPS - NUMBER OF STEPS AT WHICH TO RUN THE MODEL
CM UARG - FORCING DATA
CM ALPHA - PARAMETERS DICTATING THE GROWTH OF THE SNOW GRAINS
CM T0 - INITIAL TIME DATA
CM X0 - INITIAL AUXILIARY DATA NEEDED TO INITIALIZE THE MODEL
CM Y0 - INITIAL SNOW STATE VARIABLES
CM Z - THE LATITUDE AND LONGITUDE OF THIS RUN
CM YOUT - THE OUTPUT ARRAY CONTAINING THE SNOW DATA
subroutine ssib3(IVEG_TYPE,n_u,n_steps,uarg,alpha,t0,
& x0,y0,z_in,yout,PIXEL,n_y,n_x,xout,zout,tend,replicate,rank,
& meas,n_a,f_veg,vcov_dat)
COMMON /XRIB/RIB,temprib
character*7 run
INCLUDE 'comsib.in'
include 'snow4.in'
cm these lines to define sizes for arrays now used to initialize
cm ssib and call it as a subroutine
integer n_u,n_steps,PIXEL,n_y,n_x,replicate,rank,meas
real uarg(n_u,n_steps)
real alpha(n_a),t0(5),x0(12),y0(n_y),z_in(2),tend(5)
real vcov_dat
cm the dimension of f_veg is consistent with mvnrnd
real yout(n_y,n_steps),xout(n_x,n_steps),zout(n_steps),f_veg(1,1)
real f_radt(1,1)
cm print *, 'ssib3: f_radt=',f_radt(1,1)
CM THE ctlpa AND nroot VARIABLE WILL BE HARDCODED IN (read from file before)
c control stomatal resistance;
c final stomatal resistance=ctlpa * stomatal resistance
ctlpa=1.0
c control root depth, nroot=1: root depth has no control: nroot not
c =1: root depth is controled by rootp, which are read in vegin2.
nroot=1
cm print out all inputs...
c print *, iveg_type,n_u,n_steps,alpha,t0,x0,y0,z_in,PIXEL,
c & n_y,n_x,tend,replicate,rank,meas,n_a,f_veg,vcov_dat
CM SOME OF THESE FORMAT STATEMENTS WERE OBSELTE AND SO DELETED
2 format(a12)
3 format(a40)
CM INITIALIZE SOME VARIABLES
cm put in a check that n_x is 13, not 12, since i used 12 for so long...
if(n_x.ne.13)then
print *, 'error! n_x=',n_x,'; n_x should be 13...'
end if
CALL CNTROLS(MDAY,alpha,t0,x0,y0,z_in,pixel,n_a)
CALL CONSTS
CUMERROR = 0
ERRORCUM = 0
xqlast1 = 0.
xqlast2 = 0.
xqlast3 = 0.
xlastwet =0.
xlastswe = 0.
xlastcap1 = 0.
c four lines changed by mike to allow time to be read in and used
c iyr=1979
iyr=YEAR
mthst = MONTH
ndyst = MDAY
nhrst = TIME
mthend = 240
ndyend = 31
nhrend = 24
nobs = 0
XPG1 = 0.
XPG2 = 0.
c these are for my own time keeping which replaces time in forcing file
nyy=iyr
nmm=mthst
ndd=ndyst
nhh=nhrst
C ** ASSIGN VEGETATION PARAMETERS FROM INPUT ARRAYS
call vegin
CALL VEGIN1(IVEG_TYPE)
CALL VEGIN2(IVEG_TYPE,nmm,f_veg,vcov_dat)
c
cjyj--control parameter of snow model
MDLSNO=0
SD_CR=0.05
cm this part is changed by mike to allow snow density to be read in
cm snden=3.75
isnow=1
cm this park changed by mike
SDEP=CAPAC(2)*snden
call getinput
swe=CAPAC(2)
cm this part is changed by mike, to allow dz(nd) to be read in
snowdepth=dz(1)+dz(2)+dz(3)
if(mdlsno.eq.0.and.snowdepth.gt.sd_cr) then
isnow=0
call layern(tgs,0,nmm,ndd,nhh)
endif
CS 10/13/98
C ** CREATE THE COEFFICIENTS FOR SURFACE ALBEDO
XHC = 0.
XHG = 0.
XCI = 0.
XCT = 0.
XGI = 0.
XGT = 0.
XGS = 0.
NNX = 0
isnowhr=0
iday=0
cm ********************************************
cm START OF TIME LOOP
cm ********************************************
DO ICTRL=1,N_STEPS
cm extract forcing data from uarg array
c print *, ictrl
icrash=1
swdown=uarg(1,ICTRL)
cm obsnow and rainf: precip variables should have units mm/s
rainf=uarg(2,ICTRL)
obsnow=uarg(3,ICTRL)
dirdown=uarg(4,ICTRL)
psur=uarg(5,ICTRL)
tm=uarg(6,ICTRL)
qair=uarg(7,ICTRL)
windn=uarg(8,ICTRL)
winde=uarg(9,ICTRL)
afac=1.0
c if(ictrl.eq.icrash) print *, 'U=',uarg(1:9,ictrl)
cm check on whether it should be rain or snow...
if(tm.gt.273.16.and.obsnow.gt.0.0)then
rainf=obsnow
obsnow=0.0
end if
if(obsnow.gt.0.0)then
isnowhr=isnowhr+1
end if
cm revised by R.K(2014.7.7)
if(tm.le.273.16.and.rainf.gt.0.0)then
obsnow=rainf
rainf=0.0
end if
c if(obsnow.gt.0.0)then
c isnowhr=isnowhr+1
c end if
cm conversions of forcing data
UM=(windn**2+winde**2)**0.5
TPREC=(rainf+obsnow)*3600
cm here surface pressure is converted from Pa to mbar
psur = psur / 100.0
EM = (psur*qair)/0.6220
cm update of time variables
idkmax = 31
lstmth = nmm
do 3456, ikahan=0,19
if ((nmm.eq.(4+12*ikahan)).or.
& (nmm.eq.(9+12*ikahan)).or.
& (nmm.eq.(11+12*ikahan)).or.
& (nmm.eq.(6+12*ikahan))) then
idkmax = 30
endif
3456 continue
do 3457, ikahan=0,19
if (nmm.eq.(2+12*ikahan)) then
idkmax=28
if ((nyy.eq.1980).or.(nyy.eq.1984).or.(nyy.eq.1988).
& or.(nyy.eq.1992).or.(nyy.eq.1996).or.
& (nyy.eq.2000).or.(nyy.eq.2004)) then
idkmax = 29
endif
endif
3457 continue
if (nhh.eq.24) then
nhh = 0
ndd = ndd + 1
if (ndd.eq.(idkmax+1)) then
ndd = 1
nmm = nmm + 1
do 3458, ikahan=1,20
if (nmm.eq.(1+12*ikahan)) then
nmm = 1
nyy = nyy + 1
endif
3458 continue
endif
endif
nhh=nhh+1
qsoil=0.
wfsoil=0.
solsoil=0.
snroff=0.
if (swdown.lt.0) swdown = 0.0
if (tprec.le.0.) tprec = 0.0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (lstmth.ne.nmm) then
if ((nmm.eq.13).or.(nmm.eq.25).or.(nmm.eq.37).or.
& (nmm.eq.49).or.(nmm.eq.61).or.(nmm.eq.73).or.
& (nmm.eq.85).or.(nmm.eq.97).or.(nmm.eq.109).or.
& (nmm.eq.121).or.(nmm.eq.133).or.(nmm.eq.145).or.
& (nmm.eq.157).or.(nmm.eq.169).or.(nmm.eq.181).or.
& (nmm.eq.193).or.(nmm.eq.205).or.(nmm.eq.217).or.
& (nmm.eq.229)) rewind(2)
CALL VEGIN2(IVEG_TYPE,nmm,f_veg,vcov_dat)
endif
c if (ictrl.eq.icrash) print *, 'ictrl=',ictrl,'after timekeeping'
cm observation counter updated
nobs = nobs + 1
IF (MDLSNO.eq.0) THEN
iptype=2
cm Mike is changing this to the SNTHERM cutoff temperature for snow
If (TM.ge.273.16) iptype=1
If (TPREC.le.0.) iptype=0
END IF
c
C ** START THE MAIN PROGRAM
c
c if (ictrl.eq.icrash) print *, 'ictrl=',ictrl,'before main prog.'
RADN(3,2) = dirdown
UM = AMAX1(UM,0.25)
SWDOWN = AMAX1(SWDOWN,0.1)
IHOUR = NHH
C
c ** calculate solar zenith angle
CALL RADC2(SUNANG)
SUNANG = AMAX1(0.01,SUNANG)
c
C ** CALCULATE THE CLOUD COVER USING AN EMPIRICAL EQUATION
CLOUD = (1160.*SUNANG - SWDOWN) / (963. * SUNANG)
CLOUD = AMAX1(CLOUD,0.)
CLOUD = AMIN1(CLOUD,1.)
CLOUD = AMAX1(0.58,CLOUD)
C ** SEPERATE THE SHORT WAVE RADIATION INTO DIFFERENT SPECTRUM
DIFRAT = 0.0604 / ( SUNANG-0.0223 ) + 0.0683
IF ( DIFRAT .LT. 0. ) DIFRAT = 0.
IF ( DIFRAT .GT. 1. ) DIFRAT = 1.
DIFRAT = DIFRAT + ( 1. - DIFRAT ) * CLOUD
VNRAT = ( 580. - CLOUD*464. ) / ( ( 580. - CLOUD*499. )
& + ( 580. - CLOUD*464. ) )
C
FRAC(1,1) = (1.-DIFRAT)*VNRAT
FRAC(1,2) = DIFRAT*VNRAT
FRAC(2,1) = (1.-DIFRAT)*(1.-VNRAT)
FRAC(2,2) = DIFRAT*(1.-VNRAT)
RADN(1,1) = FRAC(1,1)*SWDOWN
RADN(1,2) = FRAC(1,2)*SWDOWN
RADN(2,1) = FRAC(2,1)*SWDOWN
RADN(2,2) = FRAC(2,2)*SWDOWN
Cm initialize snow vars -mike
IF (MDLSNO.eq.0.and.ISNOW.EQ.0) THEN
tsoil=TGS
TGS=tssno(n)
CAPAC(2)=swe
SDEP=snowdepth
ssss=sdep
else
ssss=capac(2)*snden
END IF
c if (ictrl.eq.icrash) print *, 'before radab call'
cm MAIN RADIATION CALLS
CALL RADAB (ISNOW,MDLSNO,afac,ictrl,PIXEL)
cm this is where i perturb the amount of radiation absorbed by the
cm the ground surface...
cm print *, ictrl,radt(2)
c RADT(2)=RADT(2)*exp(f_radt(1,1))
cm print *, 'radt(2)... after perturbation:',radt(2)
c if (ictrl.eq.icrash) print *, 'after radab call'
CALL ROOT1
CALL STOMA1
RSTUN = RST(1)
C
c PPC = TPREC
c PPL = TPREC-PPC
PPL = TPREC
PPC = TPREC-PPL
C
C ** WATER BALANCE CHES.
c
TOTWB = WWW(1) * POROS * ZDEPTH(1)
& + WWW(2) * POROS * ZDEPTH(2)
& + WWW(3) * POROS * ZDEPTH(3)
& + CAPAC(1) + CAPAC(2)
thelastmoist = www(1)*poros*zdepth(1) + www(2)*poros*zdepth(2)
& + www(3)*poros*zdepth(3)
thelastcapac1 = CAPAC(1)
thelastcapac2 = CAPAC(2)
cm if(meas.eq.2) then
cm print *, 'www=',www,'capac=',capac,'totwb=',totwb
cm end if
C
c ** interception and runoff calculations
CS Chanfe INTERC to INTERCS (****,***) ON 10/13/98
c if (ictrl.eq.icrash) print *, 'before intercs calls'
CALL INTERCS (ISNOW,p0,CSOIL,dzsoil,CHISL)
c if (ictrl.eq.icrash) print *, 'after intercs calls; ',
c & 'isnow=',isnow
IF (MDLSNO.eq.0.and.ISNOW.EQ.0) THEN
prcp=p0
tkair=TM
c if (ictrl.eq.icrash) print *, 'before getmet calls'
CALL getmet (iptype,UM,nmm,nhh,ndd,ictrl,pixel)
c ** aerodynamic resistance and flux calculations
solar=0.
DO 1100 IVEG = 2, 2
DO 1100 IWAVE = 1, 2
DO 1100 IRAD = 1, 2
solar=solar+RADFAC(IVEG,IWAVE,IRAD)*RADN(IWAVE,IRAD)
1100 CONTINUE
CALL snow1st (dtt,TM,solsoil,ISNOW,nmm,ndd,nhh,ictrl,pixel,imike)
c if (ictrl.eq.icrash) print *, 'before temrs2 calls'
CALL TEMRS2 (MDLSNO,ISNOW,CHISL,tsoil,solsoil,meas,
& CSOIL,dzsoil,wfsoil,ictrl,pixel,y0,n_y,replicate,rank)
c if (ictrl.eq.icrash) print *, 'after temrs2 calls'
call old
ELSE
cm if there is snowfall, find snowfall density (get_met)and
cm incorporate snowfall in the rest of the snowpack (newsnow)
c if (ictrl.eq.icrash) print *, 'top of if-structure, isnow=1',
c & 'tprec=',tprec,'iptype=',iptype,'gsize=',gsize
if(TPREC.gt.0.and.iptype.ne.1)then
prcp=TPREC/1000.
tkair=TM
c if (ictrl.eq.icrash) print *, 'before getmet'
call getmet(iptype,UM,nmm,ndd,nhh,ictrl,pixel)
c if (ictrl.eq.icrash) print *, 'before newsnow'
cm The 'getmet' subroutine sets 'iptype' to 0 if the depth of
cm newsnowfall is less than some critical value.
if(iptype.ne.0) call newsnow(ISNOW,nmm,ndd,nhh,ictrl,pixel)
c if (ictrl.eq.icrash) print *, 'after getmet, newsnow calls'
endif
if(gsize>0.) call graingrowth(ISNOW)
c if (ictrl.eq.icrash) print *, 'before temrs1 calls'
CALL TEMRS1 (MDLSNO,ISNOW,rank,replicate,pixel,ictrl,icrash)
c if (ictrl.eq.icrash) print *, 'after temrs1 calls'
END IF
cm UPDATE STATE VARIABLES AFTER N-R TEMPERATURE SOLUTION
CALL UPDAT1 (MDLSNO,ISNOW,wfsoil,swe,snroff)
c if (ictrl.eq.icrash) print *, 'after updat1 calls'
cm RELAYER SNOWPACK DEPENDING ON CHANGE IN DEPTH
IF (MDLSNO.eq.0.and.ISNOW.EQ.0) THEN
CAPAC(2)=swe
cm mike is adding this on jan 31 05 to solve an issue
cm that comes up when the snow melts completely away from a
cm pack with snowdepth > 0.05 m
snowdepth=dzo(1)+dzo(2)+dzo(3)
If (snowdepth.lt.SD_CR) Then
ISNOW=1
call LAYER1 (CSOIL,TGS,dzsoil,h,w,snowdepth,
& swe,stemp,nd,gsize,gdia)
Else
ISNOW=0
call modnodenew
End if
ELSE IF(MDLSNO.eq.0.and.ISNOW.GT.0) THEN
If (capac(2)*snden.gt.SD_CR) Then
swe=CAPAC(2)
snowdepth=capac(2)*snden
ISNOW=0
CALL LAYERN (TGS,1,nmm,nhh,ndd)
Else
ISNOW=1
cm in order to output the snow, update snowdepth variable
snowdepth = capac(2)*snden
swe = capac(2)
cm end of mike's changes
End if
END IF
c if (ictrl.eq.icrash) print *, 'after layer calls'
call set0
c if (ictrl.eq.icrash) print *, 'before water balance checks'
ROFF=ROFF+snroff
c this next line was added by dr xue
c rnoffs(nobs) = rnoffs(nobs) + snroff
C
cm WATER AND ENERGY BALANCE CHECKS
ENDWB = WWW(1) * POROS * ZDEPTH(1)
& + WWW(2) * POROS * ZDEPTH(2)
& + WWW(3) * POROS * ZDEPTH(3)
& + CAPAC(1) + CAPAC(2) - PPL/1000. + ETMASS/1000. + ROFF
ERROR = TOTWB - ENDWB
cm if(meas.eq.2) then
cm print *, 'www=',www,'capac=',capac,'ppl=',ppl,'etmass=',etmass,
cm & 'roff=',roff,'endwb=',endwb,'poros=',poros,'zdepth=',zdepth
cm end if
CUMERROR = CUMERROR + ERROR
if ((nmm.eq.240).and.(ndd.eq.31).and.(nhh.eq.24))
& write(6,*) 'CUMULATIVE ERROR = ',CUMERROR
if (abs(error).gt.0.0001) then
print*, ndd,nhh,'STOPPED Water out of balance by ',error,
& 'replicate=',replicate,'rank=',rank,'pixel=',pixel,'ictrl=',
& ictrl
STOP
endif
c if (ictrl.eq.49) print *, 'before energy balance checks, dtt=',
c & dtt
CBAL = RADT(1) - CHF - (ECT+HC+ECI)/DTT
GBAL = RADT(2) - SHF - (EGT+EGI+HG+EGS)/DTT
ZLHS = RADT(1) + RADT(2) - CHF - SHF
ZRHS = HFLUX + (ECT + ECI + EGT + EGI + EGS)/DTT
RORRE = ZLHS-ZRHS
ERRORCUM = ERRORCUM + RORRE
c if (ictrl.eq.49) print *, 'after energy balance calcs'
IF(ABS (ZLHS - ZRHS) .GT. 0.010) THEN
print *, 'warning: energy balance error greater than 0.01'
IF(ABS (ZLHS - ZRHS) .GT. 4.000) THEN
XPG1 = XPG1 + ZLHS
XPG2 = XPG2 + ZRHS
ectw=ect/dtt
eciw=eci/dtt
egtw=egt/dtt
egiw=egi/dtt
egsw=egs/dtt
print *, nmm,ndd,nhh,swdown,radn(3,2),zlwup,
+ HFLUX,CHF,SHF,ectw,eciw,egtw,egiw,egsw,temprib,
+ zlhs,zrhs
print *, 'ssib3: energy balance error. snowdepth=',snowdepth,
& 'bwo=',bwo,'tssno=',tssno,'tgs=',tgs,'gdia=',gdia,'capac(2)=',
& capac(2),'flo=',flo,'pixel=',pixel,'replicate=',replicate,
& 'rank=',rank,'y0=',y0,'obsnow=',obsnow,'egi_old=',x0(9),
& 'Tair=',uarg(6,ictrl),'ictrl=',ictrl,'Tair(i-1)=',
& uarg(6,ictrl-1),'x0=',x0,'zlat=',zlat,'zlong=',zlong,
& 'LWDOWN=',uarg(4,ictrl),'radt(1)=',radt(1),'radt(2)=',radt(2)
STOP
END IF
END IF
c if (ictrl.eq.49) print *, 'after balance checks'
c ** prepare for post-processing
c
c
wav = sqrt(www(1)*www(2))
tenest = phsat * wav ** ( - bee)
tenes1 = phsat * (www(1) ** ( - bee) )
tenes2 = phsat * (www(2) ** ( - bee) )
c
c ** store the output data for post-processing
c
c if (ictrl.eq.49) print *, 'before x* calcs, hlat=',
c & hlat
xhc = hc / dtt + xhc
xhg = hg / dtt + xhg
xci = eci /hlat + xci
xct = ect /hlat + xct
xgi = egi /hlat + xgi
xgt = egt /hlat + xgt
xgs = egs /hlat + xgs
cm CALCULATE OVERALL LATENT HEAT
siblh = (((ECT + ECI + EGT + EGI + EGS)/HLAT
& *(3150.19 - 2.378 * tm)) /3.6)
cs for france data albm=0.7, so Sun add folowing statement albm=0.4
cs on 03/29/99
c albm=0.7
cs sun add above statement on 03/29/99
xmustar=0.7
qmsensh=0.7
albm=0.7
cs sun add above statement on 03/29/99
c if (ictrl.eq.49) print *, 'before output1 call...'
call output1(sunang,siblh,isnow,iday,n_steps,yout,n_y,xout,n_x,
& zout,tend,pixel)
c if (ictrl.eq.49) print *, 'after output1 call...'
END DO
return
end
cm **************************************************************
cm END OF TIME LOOP
cm **************************************************************
C=======================================================================
C
SUBROUTINE CNTROLS(MDAY,alpha,t0,x0,y0,z_in,ipixel,n_a)
C 1 AUGUST 1988
C=======================================================================
C
C INITIALISATION AND SWITCHES.
C
C-----------------------------------------------------------------------
include 'comsib.in'
include 'snow4.in'
dimension alpha(n_a),t0(5),x0(12),y0(14),z_in(2)
C
cm READ(4,*)
cm READ(4,*) DTT, ITRUNK, ILW
cm READ(4,*) ZLAT, ZLONG, TIME,MONTH,MDAY,DAY,YEAR, NITER
cm READ(4,*) ssisnow, snden
cm READ(4,*) TC, TGS, TD, TA, TM, HT, QA
cm READ(4,*) WWW(1),WWW(2),WWW(3),CAPAC(1),CAPAC(2)
cm READ(4,*) DZ(1), DZ(2), DZ(3)
cm READ(4,*) tssn(1),tssn(2),tssn(3)
cm READ(4,*) bw(1), bw(2), bw(3)
cm READ(4,*) fl(1), fl(2), fl(3)
cm READ(4,*) gsize,gdia(1),gdia(2),gdia(3)
cm READ(4,*) sg1,sg2,egi
CM INSTEAD OF READING THESE DATA FROM FILE, THEY ARE PASSED TO SSIB
CM IN THE ALPHA,T0,X0,Y0,Z ARRAYS OR HARDCODED IN
CM FIRST EXTRACT ABOVE ARRAYS TO COMMON VARS
do i=1,3
sg1(i)=alpha(i)
end do
sg2=alpha(4)
sg3=alpha(5)
sg4=alpha(6)
time=t0(1)
month=int(t0(2))
mday=int(t0(3))
day=t0(4)
year=t0(5)
TD=x0(1)
TC=x0(2)
TA=x0(3)
TM=x0(4)
WWW(1)=x0(5)
WWW(2)=x0(6)
WWW(3)=x0(7)
CAPAC(1)=x0(8)
egi=x0(9)
gsize=x0(10)
snden=x0(11)
capac(2)=x0(12)
cm note! there is no need to set roff to x0(13)...
snowdepth=y0(1)
bw(1)=y0(2)
bw(2)=y0(3)
bw(3)=y0(4)
tssn(1)=y0(5)
tssn(2)=y0(6)
tssn(3)=y0(7)
cm initialize fl as 0, then set only if bw>0
fl(1)=0.
fl(2)=0.
fl(3)=0.
if(bw(1).gt.0.) fl(1)=y0(8)/bw(1)*1000
if(bw(2).gt.0.) fl(2)=y0(9)/bw(2)*1000
if(bw(3).gt.0.) fl(3)=y0(10)/bw(3)*1000
gdia(1)=y0(11)
gdia(2)=y0(12)
gdia(3)=y0(13)
tgs=y0(14)
zlat=z_in(1)
zlong=z_in(2)
tssno=0.
flo=0.
bwo=0.
CM DETERMINE DZ(1,2,3) FROM SNOWDEPTH USING LAYERING RULE
IF (snowdepth.le.0.05) THEN
dz=0.0
ELSE IF (snowdepth.gt.0.05.and.snowdepth.le.0.06) THEN
dz(1)=0.02
dz(2)=0.02
dz(3)=snowdepth- dz(1)- dz(2)
ELSE IF ( snowdepth.gt.0.06.and.snowdepth.le.0.08) then
dz(3)=0.02
dz(2)=0.02
dz(1)=snowdepth- dz(3)- dz(2)
ELSE IF ( snowdepth.gt.0.08.and.snowdepth.le.0.62) then
dz(3)=0.02
dz(2)=(snowdepth- dz(3))*0.33333333
dz(1)=(snowdepth- dz(3))*0.66666667
ELSE IF ( snowdepth.gt.0.62) then
dz(3)=0.02
dz(2)=0.20
dz(1)=snowdepth- dz(3)- dz(2)
End IF
cm recompute capac(2) based on new snow variables, UNLESS, we have a
cm snowpack that is greater than 0, but less than the critical depth.
cm in that case, capac(2) is used to ke
If(snowdepth.eq.0.)Then
capac(2)=0.
Else If (snowdepth.gt.0.05) Then
capac(2)=(dz(1)*bw(1)+dz(2)*bw(2)+dz(3)*bw(3))/1000.
Else
cm set up one layer variables
cm this old way was to try to be fancy and preserve the densities
cm determined by the update. but i think it caused problems.
cm the new way (directly below) resets the snow density to 4.3
cm capac(2)=snowdepth*(bw(1)+bw(2)+bw(3))/3./1000.
cm snden=snowdepth/capac(2)
snden=4.3
capac(2)=snowdepth/4.3
End IF
CM ENTER HARDCODED VARIABLES
DTT=3600.
c ITRUNK=20
ITRUNK=40
ILW=1
NITER=12
SSISNOW=0.04
HT=999.
QA=999.
C
RETURN
END
C=======================================================================
C
SUBROUTINE CONSTS
C 1 AUGUST 1988
C=======================================================================
C
C INITIALIZATION OF PHYSICAL CONSTANTS
C
C-----------------------------------------------------------------------
include 'comsib.in'
C
PSUR = 1000.
CPAIR = 1010.
RHOAIR = 1.225
cm is this supposed to be the stefan-boltzman constant?
STEFAN = 5.669 * 10E-9
G = 9.81
VKC = 0.41
PIE = 3.14159265
TIMCON = PIE/86400.
CLAI = 4.2 * 1000. * 0.2
CW = 4.2 * 1000. * 1000.
TF = 273.16
C-----------------------------------------------------------------------
C N.B. : HLAT IS EXPRESSED IN J KG-1
C SNOMEL IS EXPRESSED IN J M-1
C-----------------------------------------------------------------------
HLAT = ( 3150.19 - 2.378 * TM ) * 1000.
SNOMEL = 370518.5 * 1000.
PSY = CPAIR / HLAT * PSUR / .622
C
RETURN
END
C=====================================================================
C
SUBROUTINE CROUT(A,N,M)
C DECEM 1988
C=====================================================================
DIMENSION A(N,M)
DO 11 I=2,N
11 A(1,I)=A(1,I)/A(1,1)
DO 3 K=2,N
DO 2 IK=K,N
DO 2 MI=2,K
2 A(IK,K)=A(IK,K)-A(IK,MI-1)*A(MI-1,K)
J1=K+1
DO 3 J=J1,N
DO 4 MJ=2,K
4 A(K,J)=A(K,J)-A(K,MJ-1)*A(MJ-1,J)
3 A(K,J)=A(K,J)/A(K,K)
I1=N+1
DO 5 I=I1,M
5 A(1,I)=A(1,I)/A(1,1)
DO 7 JJ=I1,M
DO 7 L=2,N
DO 8 KL=2,L
8 A(L,JJ)=A(L,JJ)-A(L,KL-1)*A(KL-1,JJ)
7 A(L,JJ)=A(L,JJ)/A(L,L)
DO 10 JI=I1,M
DO 10 K1=2,N
K2=N-K1+2
DO 10 K3=K2,N
K4=N-K1+1
10 A(K4,JI)=A(K4,JI)-A(K4,K3)*A(K3,JI)
RETURN
END
C======================================================================
C
SUBROUTINE DELRN ( RNCDTC, RNCDTG, RNGDTG, RNGDTC )
C
C======================================================================
C
C PARTIAL DERIVATIVES OF RADIATIVE AND SENSIBLE HEAT FLUXES
C
C----------------------------------------------------------------------
include 'comsib.in'
C
TC3 = TC * TC * TC
TG3 = TGS * TGS * TGS
FAC1 = ( 1. - ALBEDO(1,3,2) ) * ( 1.-THERMK ) * VCOVER(1)
FAC2 = 1. - ALBEDO(2,3,2)
C
RNCDTC = - 2. * 4. * FAC1 * STEFAN * TC3
RNCDTG = 4. * FAC1 * FAC2 * STEFAN * TG3
C
RNGDTG = - 4. * FAC2 * STEFAN * TG3
RNGDTC = 4. * FAC1 * FAC2 * STEFAN * TC3
C
RETURN
END
C======================================================================
C
SUBROUTINE DELHF ( HCDTC, HCDTG, HGDTG, HGDTC )
C
C======================================================================
C
C PARTIAL DERIVATIVES OF SENSIBLE HEAT FLUXES
C
C----------------------------------------------------------------------
include 'comsib.in'
C
RCP = RHOAIR * CPAIR
D1 = 1./RA + 1./RB + 1./RD
TA = ( TGS/RD + TC/RB + TM/RA ) / D1
C
HC = RCP * ( TC - TA ) / RB * DTT
HG = RCP * ( TGS - TA ) / RD * DTT
C----------------------------------------------------------------------
C N.B. FLUXES EXPRESSED IN JOULES M-2
C----------------------------------------------------------------------
C
HCDTC = RCP / RB * ( 1./RA + 1./RD ) / D1
HCDTG = - RCP / ( RB * RD ) / D1
C
HGDTG = RCP / RD * ( 1./RA + 1./RB ) / D1
HGDTC = - RCP / ( RD * RB ) / D1
C
RETURN
END
C======================================================================
C
SUBROUTINE DELEF ( ECDTC, ECDTG, EGDTG, EGDTC, DEADTC, DEADTG,
& EC, EG, WC, WG, FC, FG, HR,MDLSNO,ISNOW)
C
C======================================================================
C
C PARTIAL DERIVATIVES OF LATENT HEAT FLUXES
C
C----------------------------------------------------------------------
include 'comsib.in'
C
C RADD = 44.
RCP = RHOAIR * CPAIR
C----------------------------------------------------------------------
C MODIFICATION FOR SOIL DRYNESS : HR = REL. HUMIDITY IN TOP LAYER
C----------------------------------------------------------------------
C
HRR = HR
IF ( FG .LT. .5 ) HRR = 1.
C
RCC = RST(1)*FC + 2. * RB
COC = (1.-WC)/RCC + WC/(2.*RB)
RG = RST(2)*FG
IF (MDLSNO.eq.0.and.ISNOW.eq.0) THEN
RSURF=RSOIL
ELSE
RSURF = RSOIL*FG
END IF
COG1 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)*HRR
& + VCOVER(2)/(RSURF+RD+44.)*HRR
COG2 = VCOVER(2)*(1.-WG)/(RG+RD)+(1.-VCOVER(2))/(RSURF+RD)
& + VCOVER(2)/(RSURF+RD+44.)
COG1 = COG1 + WG/RD * VCOVER(2)
COG2 = COG2 + WG/RD * VCOVER(2)
C
D2 = 1./RA + COC + COG2
TOP = COC * ETC + COG1 * ETGS + EM/RA
EA = TOP / D2
C
EC = ( ETC - EA ) * COC * RCP/PSY * DTT
C
EG = ( ETGS*COG1 - EA*COG2 ) * RCP/PSY * DTT
DEADTC = GETC * COC / D2
DEADTG = GETGS * COG1 / D2
C
ECDTC = ( GETC - DEADTC ) * COC * RCP / PSY
ECDTG = - DEADTG * COC * RCP / PSY
C
EGDTG = ( GETGS*COG1 - DEADTG*COG2 ) * RCP / PSY
EGDTC = - DEADTC * COG2 * RCP / PSY
C
RETURN
END
C
C====================================================================
C
SUBROUTINE FIT2(X1,Y1,DD,K,MM)
C DECEM 1989
C=====================================================================
PARAMETER ( NN = 30 )
DIMENSION X1(NN), Y1(11,NN), Y2(NN), SS(3,4),DD(11,3)
IOUT4 = 48
DO 50 I = 1, 3
DO 50 J = 1, 4
50 SS(I,J) = 0.
DO 100 I = 1, MM
SS(1,2) = SS(1,2) + X1(I)
XX = X1(I) * X1(I)
SS(1,3) = SS(1,3) + XX
SS(2,3) = SS(2,3) + XX * X1(I)
SS(3,3) = SS(3,3) + XX * XX
SS(1,4) = SS(1,4) + Y1(K,I)
SS(2,4) = SS(2,4) + Y1(K,I) * X1(I)
SS(3,4) = SS(3,4) + Y1(K,I) * XX
100 CONTINUE
SS(1,1) = MM
SS(2,1) = SS(1,2)
SS(2,2) = SS(1,3)
SS(3,1) = SS(2,2)
SS(3,2) = SS(2,3)
CALL CROUT(SS,3,4)
AA = SS(1,4)
BB = SS(2,4)
CC = SS(3,4)
DO 200 I = 1, MM
Y2(I) = AA + BB * X1(I) + CC * X1(I) * X1(I)
200 CONTINUE
DD(K,1) = AA
DD(K,2) = BB
DD(K,3) = CC
XY2 = 0.
DO 220 I = 1, MM
220 XY2 = XY2 + (Y2(I) - Y1(K,I))**2
XY2 = SQRT(XY2/MM)
RETURN
END
C =====================================================================
c
SUBROUTINE INTERCS (ISNOW,p0,CSOIL,DZSOIL,CHISL)
C 1 AUGUST 1988
C=======================================================================
C
C CALCULATION OF (1) INTERCEPTION AND DRAINAGE OF RAINFALL AND SNOW
C (2) SPECIFIC HEAT TERMS FIXED FOR TIME STEP
C
C MODIFICATION 30 DEC 1985 : NON-UNIFORM PRECIPITATION
C ------------ CONVECTIVE PPN. IS DESCRIBED BY AREA-INTENSITY
C RELATIONSHIP :-
C
C F(X) = A*EXP(-B*X)+C
C
C THROUGHFALL, INTERCEPTION AND INFILTRATION
C EXCESS ARE FUNCTIONAL ON THIS RELATIONSHIP
C AND PROPORTION OF LARGE-SCALE PPN.
C----------------------------------------------------------------------
include 'comsib.in'
C
DIMENSION CAPACP(2), SNOWP(2), PCOEFS(2,2)
DATA PCOEFS(1,1)/ 20. /, PCOEFS(1,2)/ .206E-8 /,
& PCOEFS(2,1)/ 0.0001 /, PCOEFS(2,2)/ 0.9999 /, BP /20. /
C
AP = PCOEFS(2,1)
CP = PCOEFS(2,2)
TOTALP = PPC + PPL
IF(TOTALP.LT.1.E-8)GO TO 6000
AP = PPC/TOTALP * PCOEFS(1,1) + PPL/TOTALP * PCOEFS(2,1)
CP = PPC/TOTALP * PCOEFS(1,2) + PPL/TOTALP * PCOEFS(2,2)
6000 CONTINUE
C
ROFF = 0.
THRU = 0.
FPI = 0.
C
C----------------------------------------------------------------------
C THERMAL CONDUCTIVITY OF THE SOIL, TAKING INTO ACCOUNT POROSITY
C----------------------------------------------------------------------
C
THETA=WWW(1)*POROS
CHISL=( 9.8E-4+1.2E-3*THETA )/( 1.1-0.4*THETA )
CHISL=CHISL*4.186E2
C
C----------------------------------------------------------------------
C THERMAL DIFFUSIVITY AND HEAT CAPACITYOF THE SOIL
C----------------------------------------------------------------------
C
DIFSL=5.E-7
C
ROCS =CHISL/DIFSL
D1 =SQRT(DIFSL*86400.0)
CSOIL=ROCS*D1/SQRT(PIE)/2.0
C YX2002 (test2)
dzsoil=D1/SQRT(PIE)/2.0
THALAS=0.
OCEANS=0.
POLAR=0.
CSOIL=CSOIL*(1.0-THALAS)+10.E10*OCEANS+POLAR*3.6*4.2E4
C
C
P0 = TOTALP * 0.001
C
C----------------------------------------------------------------------