-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathETPHR.for
1355 lines (1138 loc) · 48 KB
/
ETPHR.for
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 ETPHR, Subroutine, N.B. Pickering
C Computes canopy ET (mm/h) and gross photosynthesis (痠ol CO2/m2/s)
C for each hour.
C-----------------------------------------------------------------------
C REVISION HISTORY
C 12/05/90 NBP Written
C 11/05/93 KJB Added layers for SLW and leaf N
C 11/23/93 NBP Removed SLW effect from INVEG, all in PGLFEQ
C 12/11/93 NBP Removed AWES1 from calc. of Rsoil.
C 04/21/94 NBP Check in CANPET prevents small trans. amounts for LAI=0.
C 01/10/00 NBP Modular version
C-----------------------------------------------------------------------
C Called from: ETPHOT
C Calls: CANPET,CANOPG,HSOILT
C=======================================================================
SUBROUTINE ETPHR(
& CANHT, CEC, CEN, CLOUDS, CO2HR, DAYTIM, !Input
& DLAYR2, DULE, FNPGL, FNPGN, FRACSH, FRSHV, !Input
& KDIRBL, LAISH, LAISHV, LAISL, LAISLV, LLE, !Input
& LMXREF, LNREF, LWIDTH, MEEVP, MEPHO, NLAYR, !Input
& NSLOPE, PARSH, PARSUN, QEREF, RABS, RCUTIC, !Input
& REFHT, RHUMHR, RNITP, RWUH, SHCAP, SLAAD, !Input
& SLWREF, SLWSLO, STCOND, SWE, TAIRHR, TA, !Input
& TMIN, TYPPGL, TYPPGN, WINDHR, XHLAI, !Input
& XLMAXT, YLMAXT, !Input
& AGEFAC, EHR, LFMXSH, LFMXSL, PCNLSH, PCNLSL, !Output
& PGHR, SLWSH, SLWSL, T0HR, TCAN, THR, TSHR, !Output
& TSURF) !Output
! ------------------------------------------------------------------
USE ModuleDefs !Definitions of constructed variable types,
! which contain control information, soil
! parameters, hourly weather data.
IMPLICIT NONE
SAVE
CHARACTER MEEVP*1,MEPHO*1,TYPPGN*3,TYPPGL*3
INTEGER ITER,NLAYR
LOGICAL DAYTIM,REPEAT,STRESS
REAL AGEFAC,CANHT,CEC,CEN,CLOUDS,CO2HR,CONDSH,
& CONDSL,CSHPRV,CSHSTR,CSLPRV,CSLSTR,DLAYR2(NL),DULE,EHR,EMAXTC,
& EMAXTR,ERRBND,ERRTC,ERRTR,FNPGN(4),FNPGL(4),FRACSH,FRSHV,HOLD,
& KDIRBL,LAISH,LAISHV,LAISL,LAISLV,LFMXSH,LFMXSL,LLE,LNREF,
& LWIDTH,LMXREF,NSLOPE,PARSH,PARSUN(3),PCNLSH,PCNLSL,PGHR,QEREF,
& RA,RABS(3),RCUTIC,REFHT,RHUMHR,RNITP,RWUH,SHCAP(NL),SLAAD,
& SLWREF,SLWSH,SLWSL,SLWSLO,STCOND(NL),SWE,T0HR,TAIRHR,TA,TMIN,
& TCAN,TCPREV,THR,TPREV,TSHR(NL),TSUM,TSURF(3,1),USTAR,
& WINDHR,XHLAI,XLMAXT(6),YLMAXT(6)
PARAMETER (ERRBND=0.01)
C Initialize.
EMAXTC = ERRBND * 30.0
EMAXTR = ERRBND * 1.0
AGEFAC = 1.0
EHR = 0.0
LFMXSH = 0.0
LFMXSL = 0.0
PCNLSH = 0.0
PCNLSL = 0.0
PGHR = 0.0
RA = 0.0
USTAR = 0.0
SLWSH = 0.0
SLWSL = 0.0
T0HR = 0.0
TCAN = TAIRHR
THR = 0.0
TSURF(1,1) = TAIRHR
TSURF(2,1) = TAIRHR
TSURF(3,1) = TAIRHR
STRESS = .FALSE.
C Daylight hours with canopy.
IF (DAYTIM .AND. XHLAI .GT. 0.0) THEN
REPEAT = .TRUE.
ITER = 1
TSUM = 0.0
C Loop until evapotranspiration and photosynthesis are stable.
IF (MEEVP.EQ.'Z' .AND. MEPHO.EQ.'L') THEN
DO WHILE (REPEAT .AND. ITER .LE. 5)
TCPREV = TCAN
CALL CANOPG(
& CO2HR, FNPGL, FNPGN, LAISH, LAISL, LMXREF, !Input
& LNREF, NSLOPE, PARSH, PARSUN, QEREF, RNITP, !Input
& SLAAD, SLWSLO, TMIN, TSURF, TYPPGL, TYPPGN, !Input
& XLMAXT, YLMAXT, !Input
& AGEFAC, CONDSH, CONDSL, CSHSTR, CSLSTR, !Output
& LFMXSH, LFMXSL, PCNLSH, PCNLSL, PGHR, !Output
& SLWREF, SLWSH, SLWSL, STRESS) !Output
CALL CANPET(
& CANHT, CEC, CEN, CLOUDS, CONDSH, CONDSL, !Input
& DLAYR2, FRACSH, FRSHV, KDIRBL, LAISH, !Input
& LAISHV, LAISL, LAISLV, LWIDTH, RABS, !Input
& RCUTIC, REFHT, RHUMHR, STCOND, TAIRHR, !Input
& WINDHR, !Input
& EHR, RA, TCAN, THR, TSHR, TSURF, USTAR) !Output
TSUM = TSUM + TCAN
IF (ITER .GT. 5) THEN
TCAN = TSUM / ITER
ELSE
TCAN = (TCAN+TCPREV) / 2.0
ENDIF
ERRTC = ABS(TCAN-TCPREV)
REPEAT = ERRTC .GT. EMAXTC
ITER = ITER + 1
ENDDO
T0HR = THR
C Water stress loop. Evapotranspiration limited to by soil water
C supply by adjusting leaf conductances. Photosynthesis recalulated.
IF (THR .GT. RWUH) THEN
CSLPRV = CONDSL
C CONDSL = CONDSL * (THR-RWUH)/THR
CONDSL = CONDSL * RWUH/THR
CSHPRV = CONDSH
C CONDSH = CONDSH * (THR-RWUH)/THR
CONDSH = CONDSH * RWUH/THR
REPEAT = .TRUE.
ITER = 1
TSUM = 0.0
DO WHILE (REPEAT .AND. ITER .LE. 5)
TCPREV = TCAN
TPREV = THR
CALL CANPET(
& CANHT, CEC, CEN, CLOUDS, CONDSH, CONDSL, !Input
& DLAYR2, FRACSH, FRSHV, KDIRBL, LAISH, !Input
& LAISHV, LAISL, LAISLV, LWIDTH, RABS, !Input
& RCUTIC, REFHT, RHUMHR, STCOND, TAIRHR, !Input
& WINDHR, !Input
& EHR, RA, TCAN, THR, TSHR, TSURF, USTAR) !Output
TSUM = TSUM + TCAN
IF (ITER .GT. 5) THEN
TCAN = TSUM / ITER
ELSE
TCAN = (TCAN+TCPREV) / 2.0
ENDIF
ERRTC = ABS(TCAN-TCPREV)
ERRTR = ABS(THR-RWUH)
HOLD = CSLPRV+(CONDSL-CSLPRV)/(THR-TPREV)*(RWUH-TPREV)
CSLPRV = CONDSL
CONDSL = MAX(1.0/RCUTIC,HOLD)
HOLD = CSHPRV+(CONDSH-CSHPRV)/(THR-TPREV)*(RWUH-TPREV)
CSHPRV = CONDSH
CONDSH = MAX(1.0/RCUTIC,HOLD)
REPEAT = ERRTR .GT. EMAXTR .AND. ERRTC .GT. EMAXTC
ITER = ITER + 1
ENDDO
CSLSTR = CONDSL
CSHSTR = CONDSH
STRESS = .TRUE.
CALL CANOPG(
& CO2HR, FNPGL, FNPGN, LAISH, LAISL, LMXREF, !Input
& LNREF, NSLOPE, PARSH, PARSUN, QEREF, RNITP, !Input
& SLAAD, SLWSLO, TMIN, TSURF, TYPPGL, TYPPGN, !Input
& XLMAXT, YLMAXT, !Input
& AGEFAC, CONDSH, CONDSL, CSHSTR, CSLSTR, !Output
& LFMXSH, LFMXSL, PCNLSH, PCNLSL, PGHR, !Output
& SLWREF, SLWSH, SLWSL, STRESS) !Output
STRESS = .FALSE.
ENDIF
ELSE
CALL CANOPG(
& CO2HR, FNPGL, FNPGN, LAISH, LAISL, LMXREF, !Input
& LNREF, NSLOPE, PARSH, PARSUN, QEREF, RNITP, !Input
& SLAAD, SLWSLO, TMIN, TSURF, TYPPGL, TYPPGN, !Input
& XLMAXT, YLMAXT, !Input
& AGEFAC, CONDSH, CONDSL, CSHSTR, CSLSTR, !Output
& LFMXSH, LFMXSL, PCNLSH, PCNLSL, PGHR, !Output
& SLWREF, SLWSH, SLWSL, STRESS) !Output
ENDIF
C Night hours or bare soil.
ELSE
IF (MEEVP .EQ. 'Z') THEN
CONDSH = 0.0
CONDSL = 0.0
REPEAT = .TRUE.
ITER = 1
TSUM = 0.0
DO WHILE (REPEAT .AND. ITER .LE. 5)
TCPREV = TCAN
CALL CANPET(
& CANHT, CEC, CEN, CLOUDS, CONDSH, CONDSL, !Input
& DLAYR2, FRACSH, FRSHV, KDIRBL, LAISH, !Input
& LAISHV, LAISL, LAISLV, LWIDTH, RABS, !Input
& RCUTIC, REFHT, RHUMHR, STCOND, TAIRHR, !Input
& WINDHR, !Input
& EHR, RA, TCAN, THR, TSHR, TSURF, USTAR) !Output
TSUM = TSUM + TCAN
IF (ITER .GT. 5) THEN
TCAN = TSUM / ITER
ELSE
TCAN = (TCAN+TCPREV) / 2.0
ENDIF
ERRTC = ABS(TCAN-TCPREV)
REPEAT = ERRTC .GT. EMAXTC
ITER = ITER + 1
ENDDO
T0HR = THR
ENDIF
ENDIF
PGHR = MAX(PGHR,0.0)
T0HR = MAX(T0HR,0.0)
THR = MAX(THR,0.0)
EHR = MAX(EHR,0.0)
C Update soil moisture in upper layer, soil and canopy-air temperature
C difference.
IF (MEEVP .EQ. 'Z') THEN
CALL HSOILT(
& DLAYR2, NLAYR, SHCAP, STCOND, TA, TSURF(3,1), !Input
& TSHR) !Output
SWE = SWE
DULE = DULE
LLE = LLE
CEN = CEN
c SWE = MAX(SWE-EHR,0.0)
c CEN = (DULE-SWE) / (DULE-LLE) * 100.0
ENDIF
RETURN
END SUBROUTINE ETPHR
C========================================================================
C CANOPG, Subroutine, K.J. Boote, J.W. Jones, G. Hoogenboom
C Computes instantaneous canopy photosynthesis (痠ol CO2/m2/s) and leaf
C CO2 conductance (cm/s) of shaded and sunlit leaves. Uses shaded
C and sunlit leaf areas, light intensity on each, and integrates
C over leaf angle classes.
C-----------------------------------------------------------------------
C REVISION HISTORY
C 05/01/89 Written
C 11/15/90 NBP Reorganized
C 11/05/93 KJB Added layers for SLW and leaf N
C-----------------------------------------------------------------------
C Called from: ETPHR
C Calls: PGLFEQ,PGLEAF
C=======================================================================
SUBROUTINE CANOPG(
& CO2HR, FNPGL, FNPGN, LAISH, LAISL, LMXREF, !Input
& LNREF, NSLOPE, PARSH, PARSUN, QEREF, RNITP, !Input
& SLAAD, SLWSLO, TMIN, TSURF, TYPPGL, TYPPGN, !Input
& XLMAXT, YLMAXT, !Input
& AGEFAC, CONDSH, CONDSL, CSHSTR, CSLSTR, !Output
& LFMXSH, LFMXSL, PCNLSH, PCNLSL, PGHR, !Output
& SLWREF, SLWSH, SLWSL, STRESS) !Output
IMPLICIT NONE
SAVE
CHARACTER TYPPGN*3,TYPPGL*3
INTEGER I
LOGICAL STRESS
REAL AGEFAC,AGMXSH,AGMXSL,CO2HR,CONDSH,CONDSL,CONSUM,CONSUN,
& LAISH,LAISL,LMXREF,LFMXSH,LFMXSL,LNREF,NSLOPE,PARSH,
& PARSUN(3),PARSL,PGHR,PGSUM,PGSUN,PGSH,PGSL,QEREF,QEFFSH,
& QEFFSL,RNITP,SLAAD,TEMPSH,TEMPSL,TSURF(3,1),FNPGN(4),
& FNPGL(4),XLMAXT(6),XHLAI,YLMAXT(6),CSLSTR,CSHSTR,SLWSL,
& SLWSH,PCNLSL,PCNLSH,SLWSLO,SLWREF,TMIN
C Initialize.
TEMPSL = TSURF(1,1)
TEMPSH = TSURF(2,1)
XHLAI = LAISL + LAISH
PARSL = PARSUN(2)
C Calculate leaf photosynthesis parameters with separate layering
C of SLW and leaf N. SLW and leaf N decrease linearly with
C increasing LAI (Sinclair et al., 1993). Assume sunlit (SL)
C leaves are physically above shaded (SH) leaves.
SLWSL = 1./SLAAD + 0.5*SLWSLO*LAISH
SLWSH = 1./SLAAD - 0.5*SLWSLO*LAISL
PCNLSL = RNITP + 0.5*NSLOPE*LAISH
PCNLSH = RNITP - 0.5*NSLOPE*LAISL
C Calulate leaf photosynthesis for SH and SL leaves.
CALL PGLFEQ(
& CO2HR, FNPGL, FNPGN, LMXREF, LNREF, QEREF, !Input
& PCNLSH, SLWSH, SLWREF, TEMPSH, TMIN, TYPPGL, !Input
& TYPPGN, XLMAXT, YLMAXT, !Input
& AGMXSH, LFMXSH, QEFFSH) !Output
CALL PGLFEQ(
& CO2HR, FNPGL, FNPGN, LMXREF, LNREF, QEREF, !Input
& PCNLSL, SLWSL, SLWREF, TEMPSL, TMIN, TYPPGL, !Input
& TYPPGN, XLMAXT, YLMAXT, !Input
& AGMXSL, LFMXSL, QEFFSL) !Output
C Gaussian integration of photosynthesis and leaf CO2 conductance
C over three leaf classes for sunlit leaves.
PGSUM = 0.0
CONSUM = 0.0
DO I=1,3
CALL PGLEAF(
& CO2HR, LFMXSL, PARSUN(I), QEFFSL, TEMPSL, !Input
& CONSUN, PGSUN) !Output
IF (I .EQ. 2) THEN
PGSUM = PGSUM + PGSUN*1.6
CONSUM = CONSUM + CONSUN*1.6
ELSE
PGSUM = PGSUM + PGSUN
CONSUM = CONSUM + CONSUN
ENDIF
ENDDO
PGSL = PGSUM / 3.6
CONDSL = CONSUM / 3.6
C Compute photosynthesis and leaf CO2 conductance for shaded leaves
CALL PGLEAF(
& CO2HR, LFMXSH, PARSH, QEFFSH, TEMPSH, !Input
& CONDSH, PGSH) !Output
C Compute canopy photosynthesis (痠ol CO2/m2/s).
IF (STRESS) THEN
IF (CONDSL .GT. 0.0) THEN
PGSL = PGSL * CSLSTR/CONDSL
ELSE
PGSL = 0.0
ENDIF
IF (CONDSH .GT. 0.0) THEN
PGSH = PGSH * CSHSTR/CONDSH
ELSE
PGSH = 0.0
ENDIF
ENDIF
PGHR = PGSL*LAISL + PGSH*LAISH
AGEFAC = (LAISL*AGMXSL+LAISH*AGMXSH) / XHLAI
RETURN
END SUBROUTINE CANOPG
C=======================================================================
C PGLFEQ, Subroutine, K.J. Boote, N.B. Pickering
C Calculates leaf quantum efficiency and maximum photosynthesis for
C leaf photosynthesis equation as a function of temp, CO2, leaf N,
C and specific leaf area.
C-----------------------------------------------------------------------
C REVISION HISTORY
C 05/01/89 KJB Written
C 12/01/90 NBP Modified
C 11/05/93 KJB Added layers for SLW and leaf N
C-----------------------------------------------------------------------
C Called from: CANOPG
C Calls: CURV,TABEX
C Notes : Standard conditions and suggested values for QEREF and LXREF.
C LXREF is for upper sunlit leaves in the canopy. Lower
C leaves will have a lower rate because of less SLW and N.
C QEREF : 30 oC, 350 無/L CO2, 2% O2, <100 痠ol/m2/S PFD
C 0.0541 痠ol/痠ol (Ehleringer and Bjorkman, 1977)
C (converted from 30 oC, 325 無/L CO2)
C LXREF: 30 oC, 350 無/L CO2, 2% O2, 2000 痠ol/m2/S PFD
C measured at SLWREF and LNREF
C BEAN=28, PEANUT=28, SOYBEAN=28 痠ol/m2/s
C========================================================================
SUBROUTINE PGLFEQ(
& CO2HR, FNPGL, FNPGN, LMXREF, LNREF, QEREF, !Input
& RNITP, SLW, SLWREF, TEMPHR, TMIN, TYPPGL, !Input
& TYPPGN, XLMAXT, YLMAXT, !Input
& AGEMXL, LFMAX, QEFF) !Output
IMPLICIT NONE
SAVE
CHARACTER TYPPGN*3,TYPPGL*3
REAL AGEMXL,AGEQE,CICA,CINT,CO2HR,CO2MAX,CO2QE,
& CURV,FNPGN(4),FNPGL(4),GAMST,LFMAX,LNREF,LMXREF,LXREF,
& O2,QEFF, QEREF,RGAS,RNITP,RT,SLW,SLWMAX,SLWREF,TABEX,TAU,
& TEMPHR,TEMPMX,TK,XLMAXT(6),YLMAXT(6),TMIN,CHILL
PARAMETER (O2=210000.0,RGAS=8.314)
C Initialization. Convert LMXREF from mgCO2/m2/s to 痠ol/m2/s.
TK = TEMPHR + 273.
RT = RGAS * TK
LXREF = LMXREF * 1000.0 / 44.0
C Temperature and CO2 effects on QEFF AND LMXREF are both modeled using
C Farquhar and Caemmerer's (1982) equation (16.60 a,b) for limiting RuBP,
C combined with the temperature effect on the specificity factor (TAU)
C and the compensation point in the absence of dark respiration (GAMST).
!CHP 4/15/03 Prevent overflow
IF (RT .GT. 1000.) THEN
TAU = EXP(-3.949 + 28990.0/RT)
GAMST = 0.5 * O2 / TAU
ELSE
TAU = 1E10
GAMST = 0.0
ENDIF
C EFFECTS ON MAXIMUM LEAF PHOTOSYNTHESIS (LMXREF).
C SLW effect on LMXREF assumed linear (Dornhoff and Shibles, 1970).
SLWMAX = SLW / SLWREF
C Temperature and non-saturating CO2.
C For the computation of LMXREF, Ci/Ca = 0.7 for CO2=350 無/L. The factor
C 7.179 scales CO2MAX to 1.0 at 30 oC and 350 無/L CO2.
C CICA = 0.4+0.6*EXP(-0.002*CO2HR)
CICA = 0.7
CINT = CICA*CO2HR + (1.0-CICA)*GAMST
CINT = MAX(CINT,GAMST)
CO2MAX = 7.179 * (CINT-GAMST) / (4.0*CINT+8.0*GAMST)
C Temperature and saturating CO2.
C Temperature effect on LMXREF at saturating CO2 via lookup table.
C Sawtooth shape i.e. linear increase to peak, then linear decrease.
C Based on analysis of LMXREF using Tenhunen's (1976) data with QEFF
C from Ehleringer and Bjorkman (1977). TEMPMX scaled to 1.0 at 30 oC.
TEMPMX = TABEX(YLMAXT,XLMAXT,TEMPHR,6)/TABEX(YLMAXT,XLMAXT,30.0,6)
C Minimum night temp effect on Pmax next day, quadratic function
C after Mike Bell, peanut. ONLY the first two numbers are used.
C KJB, 9/7/94. OKAY TO USE NIGHT MINIMUM CANOPY TEMPERATURE LATER
CHILL = CURV(TYPPGL,FNPGL(1),FNPGL(2),FNPGL(3),FNPGL(4),TMIN)
C Nitrogen effects on LMXREF. Photosynthesis is affected both by
C plant nutrition and age. Quadratic from (1) to (2). AGEMXL scaled
C to 1.0 at LNREF.
AGEMXL = CURV(TYPPGN,FNPGN(1),FNPGN(2),FNPGN(3),FNPGN(4),RNITP) /
& CURV(TYPPGN,FNPGN(1),FNPGN(2),FNPGN(3),FNPGN(4),LNREF)
C EFFECTS ON QUANTUM EFFICIENCY (QEFF).
C Temperature and non-saturating CO2.
C For the computation of QEFF, Ci/Ca = 1.0. The factor 6.225 scales CO2QE
C to 1.0 at 30 oC and 350 無/L CO2.
CINT = MAX(CO2HR,GAMST)
CO2QE = 6.225 * (CINT-GAMST) / (4.*CINT+8.*GAMST)
C Nitrogen effects on QEFF. Photosynthesis is affected both by
C plant nutrition (small in legumes) and age.
C AGEQE = 0.0094 + (1.0-EXP(-4.4*AGEMAX))
C 7/17/94. KJB INCREASED EXTENT OF N EFFECT ON QEFF. COEFFICIENT
C CHANGED FROM -4.4 TO -3.0, NOW 0.97 AT AGEMAX=0.75, 0.819 AT
C 0.50, AND 0.560 AT 0.25. ALSO, NORMALIZED. THE 0.0094 VALUE
C NOW KEEPS THE VALUE AT 0.01
C CHANGE AGAIN 9/24/95 KJB, CHANGE FROM -2.5 TO -2.0, NOW 0.923
C AT 0.8, 0.734 AT 0.5, 0.388 AT 0.2, 1.0 AT 1, 0.01 AT 0.
C
AGEQE = (0.0094 + (1.0-EXP(-2.0*AGEMXL))) /
& (0.0094 + (1.0-EXP(-2.0*1.0)))
AGEQE = MIN(MAX(AGEQE,0.0),1.0)
C Calculate QEFF and LFMAX at ambient conditions.
QEFF = QEREF * CO2QE * AGEQE
LFMAX = LXREF * SLWMAX * TEMPMX * AGEMXL * CO2MAX * CHILL
RETURN
END SUBROUTINE PGLFEQ
C=======================================================================
C PGLEAF, Subroutine, K.J.Boote, J.W.Jones, G.Hoogenboom
C Calculate instantaneous leaf photosynthesis as a function of PAR
C and leaf characteristics (痠ol/m2/s). Leaf conductance calculated
C as a function of net photosynthesis and Ci/Ca ratio (cm/s)
C-----------------------------------------------------------------------
C REVISION HISTORY
C 05/01/89 GH Written
C 12/10/90 NBP Modified to calculated leaf conductance to H2O
C 11/23/93 NBP Modified for layer input of SLW
C-----------------------------------------------------------------------
C Called from: CANOPG
C Calls:
C=======================================================================
SUBROUTINE PGLEAF(
& CO2HR, LFMAX, PARLF, QEFF, TEMPHR, !Input
& CONDLF, PGLF) !Output
IMPLICIT NONE
SAVE
REAL A,B,C,CICA,CINT,CO2HR,CCO2LF,CONDLF,CVTURE,GAMST,LFMAX,QEFF,
& PARLF,PATM,PGLF,PNLF,RGAS,RT,TAU,TEMPHR
PARAMETER (CVTURE=0.8, PATM=101300.0, RGAS=8.314)
C Initialization.
RT = RGAS * (TEMPHR+273.0)
C Calculate leaf photosynthesis (痠ol CO2/m2/s) using a non-rectangular
C hyperbola (Rabinowitch, 1951; Lommen et al, 1971; Evans and Farquhar,
C Norman and Arkebauer, Gutschick, In: Boote and Loomis, 1991)
A = CVTURE
B = (QEFF*PARLF) + LFMAX
C = QEFF * PARLF * LFMAX
! CHP Added checks for floating underflow 1/16/03
! IF (LFMAX .GT. 0.0) THEN
IF (LFMAX .GT. 0.0 .AND. (QEFF*PARLF/LFMAX) .LT. 20.) THEN
C PGLF = (B - SQRT(B**2-4.*A*C)) / (2.*A)
PGLF = LFMAX * (1.0 - EXP(-QEFF*PARLF/LFMAX))
ELSE
PGLF = MAX(LFMAX, 0.0)
ENDIF
PNLF = PGLF
C Calculate leaf CO2 conductance (mol/m2/s) using assumption of constant
C CI/CA ratio. Compensation point (GAMST) is temperature dependent.
C Leaf respiration neglected so PNLF=PGLF.
!CHP 4/15/03 Prevent overflow
IF (RT .GT. 1000.) THEN
TAU = EXP(-3.9489 + 28990.0/RT)
GAMST = 1.0E6 * 0.5 * 0.21 / TAU
ELSE
TAU = 1E10
GAMST = 0.0
ENDIF
C CICA = 0.4+0.6*EXP(-0.002*CO2HR)
CICA = 0.7
CINT = CICA*CO2HR + (1.0-CICA)*GAMST
CCO2LF = MAX(PNLF/(CO2HR-CINT),0.0)
C Convert units from mol/m2/s CO2 to m/s H20.
CONDLF = 1.6 * CCO2LF * RT / PATM
RETURN
END SUBROUTINE PGLEAF
C=======================================================================
C CANPET, Subroutine, N.B. Pickering
C Computes instantaneous canopy evapotranspiration (mm/h) and
C surface temperatures of soil and shaded and sunlit leaves (C).
C Soil temperatures also computed.
C-----------------------------------------------------------------------
C REVISION HISTORY
C 12/12/90 NBP Written
C 04/21/94 NBP Check in CANPET prevents small trans. amounts for LAI=0.
C-----------------------------------------------------------------------
C Called from: ETPHR
C Calls: ETRES,ETSOLV,RADB,VPSAT
C=======================================================================
SUBROUTINE CANPET(
& CANHT, CEC, CEN, CLOUDS, CONDSH, CONDSL, !Input
& DLAYR2, FRACSH, FRSHV, KDIRBL, LAISH, !Input
& LAISHV, LAISL, LAISLV, LWIDTH, RABS, !Input
& RCUTIC, REFHT, RHUMHR, STCOND, TAIRHR, !Input
& WINDHR, !Input
& EHR, RA, TCAN, THR, TSHR, TSURF, USTAR) !Output
! ------------------------------------------------------------------
USE ModuleDefs !Definitions of constructed variable types,
! which contain control information, soil
! parameters, hourly weather data.
IMPLICIT NONE
SAVE
INTEGER I
REAL CANHT,CONDSH,CONDSL,CEC,CEN,DLAYR1,DLAYR2(NL),
& EAIRHR,ECAN,ESATHR,EHR,ETHR,FRACSH,FRSHV,G,KDIRBL,LAISH,
& LAISL,LH,LHEAT(3,1),LHVAP,LWIDTH,PATM,PSYCON,RA,REFHT,RHUMHR,
& RCUTIC,RL(3,3),RABS(3),RNTOT,RNET(3,1),RS(3,3),SHAIR,STCND1,
& STCOND(NL),TAIRHR,TCAN,THR,TSHR1,TSHR(NL),TSURF(3,1),VHCAIR,
& VPD(3,1),VPSAT,WINDHR,CLOUDS,DAIR,DAIRD,DVAPOR,Q,SH,SHEAT(3,1),
& SHAIRD,TK,MWATER,RGAS,MAIR,LAISHV,LAISLV,RADBK(3),
& USTAR,XHLAI,ZERO
PARAMETER (RGAS=8.314,MWATER=0.01802,MAIR=0.02897,PATM=101300.0,
& SHAIRD=1005.0, ZERO=1.0E-6)
C Initialize.
XHLAI = LAISH + LAISL
STCND1 = STCOND(1)
TSHR1 = TSHR(1)
DLAYR1 = DLAYR2(1) / 100.0
C Calculate air/water properties as a function of temperature.
ESATHR = VPSAT(TAIRHR) ! Pa
EAIRHR = ESATHR * RHUMHR / 100.0 ! Pa
TK = TAIRHR + 273.0 ! K
DVAPOR = MWATER * EAIRHR / (RGAS * TK) ! kg/m3
DAIRD = MAIR * (PATM-EAIRHR) / (RGAS * TK) ! kg/m3
DAIR = DVAPOR + DAIRD ! kg/m3
Q = DVAPOR / DAIR
SHAIR = SHAIRD * (1.0+0.84*Q) ! J/kg/K
VHCAIR = DAIR * SHAIR ! J/m3/K
LHVAP = (2501.0-2.373*TAIRHR) * 1000.0 ! J/kg
PSYCON = SHAIR * PATM / (0.622*LHVAP) ! Pa/K
C Create vpd and resistance matrices.
DO I=1,3
VPD(I,1) = ESATHR - EAIRHR
ENDDO
CALL ETRES(
& CANHT, CEC, CEN, CONDSH, CONDSL, FRACSH, FRSHV, !Input
& KDIRBL, LAISH, LAISL, LWIDTH, RCUTIC, REFHT, !Input
& TAIRHR, TCAN, WINDHR, !Input
& RA, RL, RS, USTAR) !Output
C Calculate NET total by subtracting net (back) longwave radiation.
CALL RADB(
& CLOUDS, EAIRHR, FRSHV, LAISHV, !Input
& LAISLV, TAIRHR, TCAN, !Input
& RADBK) !Output
RNET(1,1) = RABS(1) - RADBK(1)
RNET(2,1) = RABS(2) - RADBK(2)
RNET(3,1) = RABS(3) - RADBK(3)
RNTOT = RNET(1,1) + RNET(2,1) + RNET(3,1)
C Solve 3-zone model for ET and E (mm/h).
CALL ETSOLV(
& DLAYR1, EAIRHR, PSYCON, RL, RNET, RS, STCND1, !Input
& TAIRHR, TSHR1, VHCAIR, VPD, !Input
& ECAN, G, LH, LHEAT, SH, SHEAT, TCAN, TSURF) !Output
ETHR = LH / LHVAP * 3600.0
EHR = LHEAT(3,1) / LHVAP * 3600.0
IF (XHLAI .LE. ZERO) THEN
TSURF(1,1) = 0.0
TSURF(2,1) = 0.0
THR = 0.0
EHR = ETHR
ELSE
THR = ETHR - EHR
ENDIF
RETURN
END SUBROUTINE CANPET
C========================================================================
C ETRES, Subroutine, N.B. Pickering
C Computes resistances to latent and sensible heat and inserts them
C into matrices for ETSOLV.
C-----------------------------------------------------------------------
C REVISION HISTORY
C 12/12/90 NBP Written
C------------------------------------------------------------------------
C Called from: CANPET
C Calls: RESBLR
C========================================================================
SUBROUTINE ETRES(
& CANHT, CEC, CEN, CONDSH, CONDSL, FRACSH, FRSHV, !Input
& KDIRBL, LAISH, LAISL, LWIDTH, RCUTIC, REFHT, !Input
& TAIRHR, TCAN, WINDHR, !Input
& RA, RL, RS, USTAR) !Output
IMPLICIT NONE
SAVE
INTEGER I,J
REAL CANHT,CEC,CEN,CONDSH,CONDSL,FRACSH,FRSHV,KDIRBL,
& LAISH,LAISL,LWIDTH,RCUTIC,RA,RB(3),REFHT,RL(3,3),
& RMAX,RS(3,3),RSSH,RSSL,RSSS,RSURF(3),TAIRHR,TCAN,
& WINDHR,XHLAI,USTAR
PARAMETER (RMAX=1.0E4)
C Initialization.
XHLAI = LAISH + LAISL
C Calculate canopy and soil boundary layer resistances.
CALL RESBLR(
& CANHT, FRACSH, FRSHV, KDIRBL, LAISH, LAISL, !Input
& LWIDTH, REFHT, TAIRHR, TCAN, WINDHR, !Input
& RA, RB, USTAR) !Output
C Calculate leaf surface resistances.
IF (XHLAI .GT. 0.0) THEN
IF (CONDSL .GT. 0.0) THEN
RSSL = 1.0/(CONDSL*LAISL) - RB(1)
RSSL = MAX(RSSL,1.0)
ELSE
RSSL = RCUTIC / LAISL
ENDIF
IF (CONDSH .GT. 0.0) THEN
RSSH = 1.0/(CONDSH*LAISH) - RB(2)
RSSH = MAX(RSSH,1.)
ELSE
RSSH = RCUTIC / LAISL
ENDIF
ELSE
C For XHLAI=0.0 set RSURF = 0.0 RBLYR = RMAX in RESBLR, so
C this means RLEAF = RMAX for both heat and vapor.
RSSL = 0.0
RSSH = 0.0
ENDIF
RSURF(1) = MIN(RSSL,RMAX)
RSURF(2) = MIN(RSSH,RMAX)
C Calculate soil surface resistance. Radiation component (a+b*rna)
C reduced to constant 0.0117 using the average RNA from Jagtap (1976).
C If RNET is included again, may need a RNMIN too.
IF (CEN .LE. CEC) THEN
RSSS = 100.0
ELSE
C RSSS = 100.0 + 154.0*(EXP(0.0117*(CEN-CEC)**1.37)-1.0)
RSSS = 100.0 + 154.0*(EXP(0.0117*(CEN-CEC)**1.275)-1.0)
ENDIF
RSURF(3) = MIN(RSSS,RMAX)
C Create resistance matrices (RS, RL).
DO I=1,3
DO J=1,3
IF (I .EQ. J) THEN
RS(I,J) = RA + RB(I)
RL(I,J) = RA + RB(I) + RSURF(I)
ELSE
RS(I,J) = RA
RL(I,J) = RA
ENDIF
ENDDO
ENDDO
RETURN
END SUBROUTINE ETRES
C=======================================================================
C RESBLR, Subroutine, N.B. Pickering
C Calculates canopy and soil boundary-layer resistances to heat and
C vapor. Based on Choudhury and Monteith (1988).
C-----------------------------------------------------------------------
C REVISION HISTORY
C 02/09/93 NBP Written
C 04/24/94 NBP Added check to prevent -ve wind speed at top of canopy.
C-----------------------------------------------------------------------
C Called from: ETRES
C Calls:
C=======================================================================
SUBROUTINE RESBLR(
& CANHT, FRACSH, FRSHV, KDIRBL, LAISH, LAISL, !Input
& LWIDTH, REFHT, TAIRHR, TCAN, WINDHR, !Input
& RA, RB, USTAR) !Output
IMPLICIT NONE
SAVE
REAL CANHT,D,ETAK,ETAKMX,ETAW,ETAWMX,FRSHV,FRACSH,H,KDIRBL,
& KH,LAISH,LAISL,ZMD,HMD,K1,K2,LWIDTH,PSIM,PSIH,RA,
& RBLF,RBSH,RBSL,RB(3),REFHT,RMAX,TAIRHR,TCAN,TKAIR,
& WINDHR,WINDSP,XHLAI,Z0H,Z0M,ZS0H,ZS0M,LHZ0M,LZZ0H,LZZ0M,DT,
& MO,USTAR,X,A,B,PI,RATIO,WINDH,RBSS,ZERO
PARAMETER (ETAKMX=2.0, ETAWMX=3.0, PI=3.14159, RMAX=1.0E4,
& ZERO=1.0E-6)
C Initialization and calculation of zero plane displacement height and
C canopy surface roughness (Brutsaert, 1982).
XHLAI = LAISH + LAISL
WINDSP = MAX(WINDHR,1.0)
H = CANHT
ZS0M = 0.03 ! m
Z0M = MAX(ZS0M,0.13*H) ! m
IF (FRSHV .GT. 0.0 .AND. FRSHV.LT. 1.0) THEN
Z0M = 4.0 * Z0M
ENDIF
D = 0.77 * H ! m
Z0H = Z0M/5.0 ! m
ZS0H = ZS0M/5.0 ! m
HMD = H - D ! m
ZMD = REFHT - D ! m
LZZ0M = LOG(ZMD/Z0M)
LZZ0H = LOG(ZMD/Z0H)
TKAIR = TAIRHR + 273.
DT = TCAN - TAIRHR
IF (ABS(DT) .LE. ZERO) THEN
USTAR = 0.4 * WINDSP / LZZ0M
RA = LZZ0H / (0.4 * USTAR)
ENDIF
C ETAK = ETAKMX*FRSHV
C ETAW = ETAWMX*FRSHV
ETAK = ETAKMX
ETAW = ETAWMX
C Stability correction for RA (Choudhury et al, 1986).
MO = -0.4*9.81*DT*ZMD / (RA*TKAIR*USTAR**3)
IF (ABS(MO) .LE. ZERO) THEN ! neutral
PSIM = 0.0
PSIH = 0.0
ELSE IF (MO .GT. ZERO) THEN ! stable
PSIM = -5.0 * MO
PSIM = MAX(PSIM,-5.0)
PSIH = PSIM
ELSE IF (MO .LT. -ZERO) THEN ! unstable
X = (1.0-16.0*MO)**0.25
A = ((1.0+X)/2.0)**2
B = (1.0+X**2)/2.0
PSIM = ALOG(A*B) - 2.0*ATAN(X) + PI/2.0
PSIH = 2.0 * ALOG(B)
RATIO = PSIH / PSIM
PSIM = MIN(PSIM,LZZ0M-0.1)
PSIH = RATIO * PSIM
ENDIF
C Aerodynamic resistance.
USTAR = 0.4 * WINDSP / (LZZ0M-PSIM)
RA = (LZZ0H-PSIH) / (0.4 * USTAR) ! s/m
C Canopy calculations.
IF (XHLAI .GT. 0.0) THEN
C Calculate windspeed at the top of the canopy.
LHZ0M = LOG(HMD/Z0M)
WINDH = MAX(WINDSP*LHZ0M/LZZ0M,0.1) ! m/s
C Calculate leaf boundary layer resistances (extension of Choudhury
C and Montieth, 1988).
RBLF = ETAW * SQRT(LWIDTH/WINDH)
& / (2.0*0.01*XHLAI*(1.0-EXP(-ETAW/2.0))) ! s/m
K1 = ETAW/2.0+KDIRBL*XHLAI/FRACSH
RBSL = SQRT(LWIDTH/WINDH) * K1 / (0.01*XHLAI*(1.0-EXP(-K1)))
RBSL = MIN(RBSL,RMAX)
RBSH = 1.0 / (1.0/RBLF-1.0/RBSL)
RBSH = MAX(RBSH,1.0)
RBSL = FRSHV * RBSL
RBSH = FRSHV * RBSL
C Calculate soil aerodynamic resistance (extension of Choudhury
C and Montieth, 1988).
KH = 0.16 * WINDSP * HMD / LZZ0M ! m/s
K2 = EXP(-ETAK*ZS0H/H) - EXP(-ETAK*(D+Z0H)/H)
RBSS = H*EXP(ETAK)*K2 / (ETAK*KH) ! s/m
RBSS = FRSHV * RBSS
C Bare soil.
ELSE
RBSL = RMAX
RBSH = RMAX
RBSS = 0.0
ENDIF
C Set boundary layer array.
RB(1) = RBSL
RB(2) = RBSH
RB(3) = RBSS
RETURN
END SUBROUTINE RESBLR
C=======================================================================
C ETSOLV, Subroutine, S. S. Jagtap, N.B. Pickering
C Computes latent and sensible heat, and temperatures for zonal ET
C model. Solves system of 10 equations using previous time step's
C top-layer soil temperature. Modified from Jagtap (1976), Jagtap
C and Jones (1989). Small errors in papers; equations here correct.
C-----------------------------------------------------------------------
C REVISION HISTORY
C ??/??/89 SSJ Written
C 01/15/91 NBP Modified
C-----------------------------------------------------------------------
C Called from: ETPHR
C Calls: GAUSSJ,MATADD,MATCON,MATPRO,VPSLOP
C=======================================================================
SUBROUTINE ETSOLV(
& DLAYR1, EAIRHR, PSYCON, RL, RNET, RS, STCND1, !Input
& TAIRHR, TSHR1, VHCAIR, VPD, !Input
& ECAN, G, LH, LHEAT, SH, SHEAT, TCAN, TSURF) !Output
IMPLICIT NONE
SAVE
REAL STCND1,DLAYR1,DZ1,EAIRHR,ECAN,G,PSYCON,RBLCN,SH,LH,
& TAIRHR,TCAN,TSHR1,VHCAIR,VP,VPSLOP,VSP,HOLD
REAL CMAT(3,1),IRL(3,3),IRS(3,3),LMAT(1,1),LHEAT(3,1),ONE(1,3),
& RL(3,3),RNET(3,1),RS(3,3),RSL(3,3),RSLRN(3,1),SMAT(1,1),
& SHEAT(3,1),TDIFF(3,1),TSURF(3,1),VHAIRS(3,3),VPD(3,1),
& VPIRLD(3,1),VSPIRL(3,3),XMAT(3,3),YMAT(3,1)
DATA ONE/1.0,1.0,1.0/
C Initialize
VSP = VHCAIR * VPSLOP(TAIRHR) / PSYCON
VP = VHCAIR / PSYCON
DZ1 = DLAYR1 / 2.0
C Invert latent and sensible heat matrices to get [iRL] and [iRS].
C Create temporary [XMAT] = [iRSL] = [C4] then adjust for G estimated from
C last E-zone temperature. [XMAT] inverted to get [RSL].
C VHAIRS and VSPIRL saved for later use.
CALL GAUSSJ(
& RL,3, !Input
& IRL) !Output
CALL GAUSSJ(
& RS,3, !Input
& IRS) !Output
CALL MATCON(
& VSP,IRL,3,3,'*', !Input
& VSPIRL) !Output
CALL MATCON(
& VHCAIR,IRS,3,3,'*', !Input
& VHAIRS) !Output
CALL MATADD(
& VHAIRS,VSPIRL,3,3,'+', !Input
& XMAT) !Output
XMAT(3,3) = XMAT(3,3) + STCND1/DZ1
CALL GAUSSJ(
& XMAT,3, !Input
& RSL) !Output
C Create [CMAT] = VP*[RSL]*[IRL]*[VPD]. Adjust [VPiRLD] = [C3] for G
C estimated from last E-zone temp. Original matrix [VPIRLD] restored
C for later use.
CALL MATPRO(IRL,VPD,3,3,1,YMAT)
CALL MATCON(VP,YMAT,3,1,'*',VPIRLD)
HOLD = VPIRLD(3,1)
VPIRLD(3,1) = HOLD - STCND1/DZ1*(TSHR1-TAIRHR)
CALL MATPRO(RSL,VPIRLD,3,3,1,CMAT)
VPIRLD(3,1) = HOLD
C Solve for TDIFF matrix.
CALL MATPRO(RSL,RNET,3,3,1,RSLRN)
CALL MATADD(VHAIRS,VSPIRL,3,3,'+',XMAT) ! Check this--XMAT used?
CALL MATADD(RSLRN,CMAT,3,1,'-',TDIFF)
CALL MATCON(TAIRHR,TDIFF,3,1,'+',TSURF)
C Solve remaining 6 equations for latent and sensible heat fluxes and
C sum up latent and sensible heats for 3 zones.
CALL MATPRO(VHAIRS,TDIFF,3,3,1,SHEAT)
CALL MATPRO(VSPIRL,TDIFF,3,3,1,YMAT)
CALL MATADD(YMAT,VPIRLD,3,1,'+',LHEAT)
CALL MATPRO(ONE,LHEAT,1,3,1,LMAT)
CALL MATPRO(ONE,SHEAT,1,3,1,SMAT)
LH = LMAT(1,1)
SH = SMAT(1,1)
C Calculate other values. RBLCN is any off-diagonal element of RS.
RBLCN = RS(1,2)
ECAN = (LH*RBLCN/VP) + EAIRHR
TCAN = (SH*RBLCN/VHCAIR) + TAIRHR
G = STCND1 / DZ1 * (TSURF(3,1)-TSHR1)
RETURN
END SUBROUTINE ETSOLV
C=======================================================================
C GAUSSJ, Subroutine, Numerical Recipes, N.B. Pickering
C Inverts N x N matrix using Gauss-Jordan elimination with full
C pivoting (Numerical Recipes, Press et al., 1986, pg. 28).
C The original matrix is AMAT and the inverse AINV.
C-----------------------------------------------------------------------
C REVISION HISTORY
C ??/??/?? Written
C 01/10/91 NBP Modified
C-----------------------------------------------------------------------
C Called from: ETSOLV
C Calls:
C=======================================================================
SUBROUTINE GAUSSJ(
& AMAT,N, !Input
& AINV) !Output