-
Notifications
You must be signed in to change notification settings - Fork 5
/
CALTSXY.for
352 lines (340 loc) · 11.8 KB
/
CALTSXY.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
SUBROUTINE CALTSXY
C
C CHANGE RECORD
C ** SUBROUTINE CALTSXY UPDATES TIME VARIABLE SURFACE WIND STRESS
C
USE GLOBAL
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::CLOUDTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::EVAPTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::PATMTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::RAINTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::RHAT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::SOLSWRTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::SVPAT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::TATMTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::TWETTT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::VPAT
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDE
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDN
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDSXX
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDSXY
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDSYX
REAL,SAVE,ALLOCATABLE,DIMENSION(:)::WINDSYY
C
IF(.NOT.ALLOCATED(CLOUDTT))THEN
ALLOCATE(CLOUDTT(NASERM))
ALLOCATE(EVAPTT(NASERM))
ALLOCATE(PATMTT(NASERM))
ALLOCATE(RAINTT(NASERM))
ALLOCATE(RHAT(NASERM))
ALLOCATE(SOLSWRTT(NASERM))
ALLOCATE(SVPAT(NASERM))
ALLOCATE(TATMTT(NASERM))
ALLOCATE(TWETTT(NASERM))
ALLOCATE(VPAT(NASERM))
ALLOCATE(WINDE(NWSERM))
ALLOCATE(WINDN(NWSERM))
ALLOCATE(WINDSXX(LCM))
ALLOCATE(WINDSXY(LCM))
ALLOCATE(WINDSYX(LCM))
ALLOCATE(WINDSYY(LCM))
CLOUDTT=0.0
EVAPTT=0.0
PATMTT=0.0
RAINTT=0.0
RHAT=0.0
SOLSWRTT=0.0
SVPAT=0.0
TATMTT=0.0
TWETTT=0.0
VPAT=0.0
WINDE=0.0
WINDN=0.0
WINDSXX=0.0
WINDSXY=0.0
WINDSYX=0.0
WINDSYY=0.0
! *** ONE TIME SPATIAL DISTRIBUTION,\
C *** OOPS, REVC & RCHC NOT SAVED FOR EACH SERIES
C DO L=2,LA
C CLEVAP(L)=0.
C CCNHTT(L)=0.
C ENDDO
C DO NA=1,NASER
C CLEVAPT=0.001*ABS(REVC)
C CCNHTTT=0.001*ABS(RCHC)
C DO L=2,LA
C CLEVAP(L)=CLEVAP(L)+ATMWHT(L,NA)*CLEVAPT
C CCNHTT(L)=CCNHTT(L)+ATMWHT(L,NA)*CCNHTTT
C ENDDO
C ENDDO
DO L=2,LA
CLEVAP(L)=0.001*ABS(REVC)
CCNHTT(L)=0.001*ABS(RCHC)
ENDDO
ENDIF
C
C**********************************************************************C
C
C INITIALIZE WIND SHELTERED SURFACE GAS TRANSFER
C
IF(N.EQ.-1.AND.NWSER.GT.0)THEN
IF(DEBUG)THEN
OPEN(1,FILE='WINDSHELT.OUT')
CLOSE(1,STATUS='DELETE')
OPEN(1,FILE='WINDSHELT.OUT')
ENDIF
DO L=2,LA
I=IL(L)
J=JL(L)
IF(WINDSTKA(L).GT.0.0)THEN
! ** IF WINDSTKA > 0 BOTH X AND Y COMPONENTS ARE APPLIED
WINDSXX(L)=CVN(L)
WINDSXY(L)=-CVE(L)
WINDSYX(L)=-CUN(L)
WINDSYY(L)=CUE(L)
ELSE
! ** IF WINDSTKA < 0 SLECTIVELY APPLY X AND Y COMPONENTS
! ** FIRST CASE IS FULLY OPEN WATER
WINDSXX(L)=CVN(L)
WINDSXY(L)=-CVE(L)
WINDSYX(L)=-CUN(L)
WINDSYY(L)=CUE(L)
LS=LSC(L)
LN=LNC(L)
! ** SECOND CASE IS 1D CHANNEL IN COMP X DIRECTION
IF(SVB(L).LT.0.5.AND.IJCT(I,J-1).NE.5)THEN
IF(SVB(LN).LT.0.5.AND.IJCT(I,J+1).NE.5)THEN
WINDSXX(L)=CVN(L)
WINDSXY(L)=-CVE(L)
WINDSYX(L)=-1000.
WINDSYY(L)=0.
ENDIF
ENDIF
! ** THIRD CASE IS 1D CHANNEL IN COMP Y DIRECTION
IF(SUB(L).LT.0.5.AND.IJCT(I-1,J).NE.5)THEN
IF(SUB(L+1).LT.0.5.AND.IJCT(I+1,J).NE.5)THEN
WINDSXX(L)=0.
WINDSXY(L)=-1000.
WINDSYX(L)=-CUN(L)
WINDSYY(L)=CUE(L)
ENDIF
ENDIF
ENDIF
IF(DEBUG)WRITE(1,1111)IL(L),JL(L),WINDSTKA(L),WINDSXX(L),
& WINDSXY(L),WINDSYX(L),WINDSYY(L)
ENDDO
IF(DEBUG)CLOSE(1)
ENDIF
1111 FORMAT(2I5,10F10.6)
C
C**********************************************************************C
C
IF(NWSER.GT.0)THEN
! *** UPDATE THE FORCING WIND DATA TO THE CURRENT TIME
DO NA=1,NWSER
IF(ISDYNSTP.EQ.0)THEN
TIME=DT*FLOAT(N)/TCWSER(NA)+TBEGIN*(TCON/TCWSER(NA))
ELSE
TIME=TIMESEC/TCWSER(NA)
ENDIF
M1=MWTLAST(NA)
MSAVE=M1
200 CONTINUE
M2=M1+1
IF(TIME.GT.TWSER(M2,NA))THEN
M1=M2
GOTO 200
ELSE
MWTLAST(NA)=M1
ENDIF
TDIFF=TWSER(M2,NA)-TWSER(M1,NA)
WTM1=(TWSER(M2,NA)-TIME)/TDIFF
WTM2=(TIME-TWSER(M1,NA))/TDIFF
DEGM1=90.-WINDD(M1,NA)
DEGM2=90.-WINDD(M2,NA)
WINDS1=WTM1*WINDS(M1,NA)+WTM2*WINDS(M2,NA)
WINDS2=WTM1*WINDS(M1,NA)+WTM2*WINDS(M2,NA)
WINDE1=WINDS(M1,NA)*COS(DEGM1/57.29578)
WINDN1=WINDS(M1,NA)*SIN(DEGM1/57.29578)
WINDE2=WINDS(M2,NA)*COS(DEGM2/57.29578)
WINDN2=WINDS(M2,NA)*SIN(DEGM2/57.29578)
WINDE(NA)=WTM1*WINDE1+WTM2*WINDE2
WINDN(NA)=WTM1*WINDN1+WTM2*WINDN2
ENDDO
! *** CALCULATE THE WIND STRESS
IF(NWSER.GT.1)THEN
DO L=2,LA
WNDVELE(L)=0.
WNDVELN(L)=0.
ENDDO
DO NA=1,NWSER
DO L=2,LA
WNDVELE(L)=WNDVELE(L)+WNDWHT(L,NA)*WINDE(NA)
WNDVELN(L)=WNDVELN(L)+WNDWHT(L,NA)*WINDN(NA)
ENDDO
ENDDO
ELSE !IF(NWSER.EQ.1)THEN
DO L=2,LA
WNDVELE(L)=WINDE(1)
WNDVELN(L)=WINDN(1)
ENDDO
ENDIF
C
DO L=2,LA
! ** CASE 0 MAGNITUDE SHELTERING AND NO DIRECTIONAL SHELTERING
IF(WINDSTKA(L).GT.0.0)THEN
WNDFAC=ABS(WINDSTKA(L))
WNDVELE(L)=WNDFAC*WNDVELE(L)
WNDVELN(L)=WNDFAC*WNDVELN(L)
WINDST(L)=SQRT( WNDVELE(L)*WNDVELE(L)
& +WNDVELN(L)*WNDVELN(L) )
C2=1.2E-6*(0.8+0.065*WINDST(L))
TSEAST=C2*WINDST(L)*WNDVELE(L)
TSNORT=C2*WINDST(L)*WNDVELN(L)
TSX(L)=WINDSXX(L)*TSEAST+WINDSXY(L)*TSNORT
TSY(L)=WINDSYX(L)*TSEAST+WINDSYY(L)*TSNORT
ELSEIF(WINDSTKA(L).LT.0.0)THEN
! ** CASE 1 MAGNITUDE SHELTERING AND DIRECTIONAL SHELTERING, OPEN WATER
IF(WINDSYX(L).GT.-99.0.AND.WINDSXY(L).GT.-99.0)THEN
WNDFAC=ABS(WINDSTKA(L))
WNDVELE(L)=WNDFAC*WNDVELE(L)
WNDVELN(L)=WNDFAC*WNDVELN(L)
WINDST(L)=SQRT( WNDVELE(L)*WNDVELE(L)
& +WNDVELN(L)*WNDVELN(L) )
C2=1.2E-6*(0.8+0.065*WINDST(L))
TSEAST=C2*WINDST(L)*WNDVELE(L)
TSNORT=C2*WINDST(L)*WNDVELN(L)
TSX(L)=WINDSXX(L)*TSEAST+WINDSXY(L)*TSNORT
TSY(L)=WINDSYX(L)*TSEAST+WINDSYY(L)*TSNORT
ENDIF
! ** CASE 2 MAGNITUDE SHELTERING AND DIRECTIONAL SHELTERING, X CHANNEL
IF(WINDSYX(L).LT.-99.0)THEN
WINDXX=WINDSXX(L)*WNDVELE(L)+WINDSXY(L)*WNDVELN(L)
WNDFAC=ABS(WINDSTKA(L))
WINDXX=WNDFAC*WNDVELE(L)
WINDST(L)=ABS(WINDXX)
TSX(L)=1.2E-6*(0.8+0.065*WINDST(L))*WINDST(L)*WINDXX
TSY(L)=0.
ENDIF
! ** CASE 3 MAGNITUDE SHELTERING AND DIRECTIONAL SHELTERING, Y CHANNEL
IF(WINDSXY(L).LT.-99.0)THEN
WINDYY=WINDSYX(L)*WNDVELE(L)+WINDSYY(L)*WNDVELN(L)
WNDFAC=ABS(WINDSTKA(L))
WINDYY=WNDFAC*WINDYY
WINDST(L)=ABS(WINDYY)
TSX(L)=0
TSY(L)=1.2E-6*(0.8+0.065*WINDST(L))*WINDST(L)*WINDYY
ENDIF
ENDIF
ENDDO
ENDIF
C
C CFTSX=1.
C CFTSY=1.
C HHUU=2.*HUWET(L)
C HHVV=2.*HVWET(L)
C
IF(NASER.GT.0)THEN
DO NA=1,NASER
IF(ISDYNSTP.EQ.0)THEN
TIME=DT*FLOAT(N)/TCASER(NA)+TBEGIN*(TCON/TCASER(NA))
ELSE
TIME=TIMESEC/TCASER(NA)
ENDIF
M1=MATLAST(NA)
100 CONTINUE
M2=M1+1
IF(TIME.GT.TASER(M2,NA))THEN
M1=M2
GOTO 100
ELSE
MATLAST(NA)=M1
ENDIF
TDIFF=TASER(M2,NA)-TASER(M1,NA)
WTM1=(TASER(M2,NA)-TIME)/TDIFF
WTM2=(TIME-TASER(M1,NA))/TDIFF
PATMTT(NA)=WTM1*PATM(M1,NA)+WTM2*PATM(M2,NA)
TATMTT(NA)=WTM1*TDRY(M1,NA)+WTM2*TDRY(M2,NA)
TWETTT(NA)=WTM1*TWET(M1,NA)+WTM2*TWET(M2,NA)
RAINTT(NA)=WTM1*RAIN(M1,NA)+WTM2*RAIN(M2,NA)
EVAPTT(NA)=WTM1*EVAP(M1,NA)+WTM2*EVAP(M2,NA)
SOLSWRTT(NA)=WTM1*SOLSWR(M1,NA)+WTM2*SOLSWR(M2,NA)
CLOUDTT(NA)=WTM1*CLOUD(M1,NA)+WTM2*CLOUD(M2,NA)
SVPAT(NA)=
& 10.**((0.7859+0.03477*TATMTT(NA))/(1.+0.00412*TATMTT(NA)))
IF(IRELH(NA).EQ.0.AND.ISTOPT(2).NE.4)THEN
C RHAT(NA)=1.
C & -0.00066*(PATMTT(NA)/SVPAT(NA))*(TATMTT(NA)-TWETTT(NA))
! *** DSLLC Begin
! *** (Correct RHA Computation from wet bulb)
TMPVAL=0.00066*(1.0+0.00115*TWETTT(NA))
SVPWET=
& 10.**((0.7859+0.03477*TWETTT(NA))/(1.+0.00412*TWETTT(NA)))
TMPVL1=SVPWET-TMPVAL*PATMTT(NA)*(TATMTT(NA)-TWETTT(NA))
RHAT(NA)=MAX(TMPVL1/ SVPAT(NA),.01)
! *** DSLLC End
ELSE
RHAT(NA)=TWETTT(NA)
ENDIF
VPAT(NA)=RHAT(NA)*SVPAT(NA)
ENDDO
C ENDIF
C IF(NASER.GT.0)THEN
IF(NASER.GT.1)THEN
DO L=2,LA
PATMT(L)=0.
TATMT(L)=0.
RAINT(L)=0.
EVAPT(L)=0.
SOLSWRT(L)=0.
CLOUDT(L)=0.
SVPA(L)=0.
RHA(L)=0.
VPA(L)=0.
ENDDO
DO NA=1,NASER
DO L=2,LA
PATMT(L)=PATMT(L)+ATMWHT(L,NA)*PATMTT(NA)
TATMT(L)=TATMT(L)+ATMWHT(L,NA)*TATMTT(NA)
RAINT(L)=RAINT(L)+ATMWHT(L,NA)*RAINTT(NA)
EVAPT(L)=EVAPT(L)+ATMWHT(L,NA)*EVAPTT(NA)
SOLSWRT(L)=SOLSWRT(L)+ATMWHT(L,NA)*SOLSWRTT(NA)
CLOUDT(L)=CLOUDT(L)+ATMWHT(L,NA)*CLOUDTT(NA)
SVPA(L)=SVPA(L)+ATMWHT(L,NA)*SVPAT(NA)
RHA(L)=RHA(L)+ATMWHT(L,NA)*RHAT(NA)
VPA(L)=VPA(L)+ATMWHT(L,NA)*VPAT(NA)
ENDDO
ENDDO
ELSE
DO L=2,LA
PATMT(L)=PATMTT(1)
TATMT(L)=TATMTT(1)
RAINT(L)=RAINTT(1)
EVAPT(L)=EVAPTT(1)
SOLSWRT(L)=SOLSWRTT(1)
CLOUDT(L)=CLOUDTT(1)
SVPA(L)=SVPAT(1)
RHA(L)=RHAT(1)
VPA(L)=VPAT(1)
ENDDO
ENDIF
! *** PMC - MOVED ALL TIME INVARIANT PARAMETERS TO KEEP FROM COMPUTING EVERY TIME
IF(REVC.LT.0.)THEN
CLEVAPTMP=0.001*ABS(REVC)
DO L=2,LA
CLEVAP(L)=1.E-3*(0.8+0.065*WINDST(L))
CLEVAP(L)=MAX(CLEVAP(L),CLEVAPTMP)
ENDDO
ENDIF
IF(RCHC.LT.0.)THEN
CCNHTTTMP=0.001*ABS(RCHC)
DO L=2,LA
CCNHTT(L)=1.E-3*(0.8+0.065*WINDST(L))
CCNHTT(L)=MAX(CCNHTT(L),CCNHTTTMP)
ENDDO
ENDIF
ENDIF
C
RETURN
END