forked from mikedurand/EnKFYampa
-
Notifications
You must be signed in to change notification settings - Fork 0
/
interfacez.f90
337 lines (300 loc) · 13.2 KB
/
interfacez.f90
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
! ----------------------------------------------------------------------
!
SUBROUTINE INTERFACEZ(Y,ZOUT,N_RL,N_R,N_Y,N_yP,N_Yr,N_Z,N_ZP,N_ZR,N_C,N_X,&
CTRL,FREQ,THETA,X,N_A,N_U,FORCING,MONTH,TBRAW,ALBEDO,MEAS,RANK,VEG_MAP,&
IERR,F_VEG,BDRF,N_BDRF,MEAS_SWITCH,VCOVER)
!
! ----------------------------------------------------------------------
!
! THIS SUBROUTINE PERFORMS SEVERAL FUNCTIONS. OVER THE LOOP OF ALL THE
! REPLICATES, 1) THE PREDICTED OBSERVATIONS ARE CALCULATED AT THE STATE
! PIXEL SCALE 2) THE PREDICTED OBSERVATIONS ARE AGGREGATED TO THE
! MEASUREMENT SCALE FOR COMPARISON WITH THE ACTUAL OBSERVATIONS
!
! Y - STATE VECTOR
! Z - OUTPUT PREDICTED OBSERVATIONS
! N_RL - NUMBER OF REPLICATES LOCAL ON THE CALLING PROCESS; THIS ARGUMENT IS
! SET TO 1 WHEN CALLING FROM ENKFGEN, SINCE WE ARE ONLY COMPUTING A SINGLE
! PREDICTED OBSERVATION
! N_R - NUMBER OF REPLICATES GLOBALLY
! N_Y - TOTAL NUMBER OF STATES
! N_rP - TOTAL NUMBER OF running PIXELS AT STATE RESOLUTION
! N_Z - TOTAL NUMBER OF OBSERVATIONS
! N_ZP - ARRAY OF TOTAL NUMBER OF MEASUREMENTS FOR EACH RESOLUTION
! N_ZR - NUMBER OF MEASUREMENTS (AT DIFFERENT RESOLUTIONS) OF EACH PIXEL
! N_C - LENGTH OF CONTROL ARRAY (CTRL)
! N_X - DIMENSION OF X VARIABLE
! N_A - DIMENSION OF ALPHA VARIABLE
! CTRL - THE ARRAY WITH PARAMETERS CONTROLLING SIMULATIONS
! FREQ - FREQUENCY VECTOR
! THETA - MEASUREMENT ANGLE VECTOR
! X - AUXILIARY STATE VARIABLES
! N_U - DIMENSION OF FORCING VARIABLE
! FORCING - FORCING DATA FOR STATE MODEL, NEEDED FOR RTMs
! MONTH AT WHICH RTM IS RUN
! TBRAW - PIXEL BRIGHTNESS TEMPS TO BE AGGREGATED
! ALBEDO - RAW ALBEDO OUTPUT
! MEAS - MEASUREMENT NUMBER
! RANK - MPI PROCESS RANK
! VEG_MAP - THE SSIB VEGETATION COVER CLASSIFICATION FOR EACH PIXEL
! IERR - MPI ERROR VARIABLE
! F_VEG - FACTOR BY WHICH VEGETATION IS SCALED (SPATIALLY UNIFORM)
! BDRF - ARRAY CONTAINING THE ALBEDO LOOKUP TABLE
! N_BDRF - NUMBER OF ROWS (ENTRIES) IN BDRF LUT
! MEAS_SWITCH - SWITCH THAT SPECIFIES WHICH MEASUREMENTS ARE USED. IF EQUAL
! TO 3, THEN NO PM MEASUREMENTS ARE CALCULATED
! VCOVER - VEGETATION COVER PERCENTAGE
!
! CALLS: SSIB_LAYER, RAD_XFER_MOD
! CALLED FROM: ENKF, ENKFGEN
!
! VERSION HISTORY: WRITTEN - MD - MAY 2005
! SET UP FOR VEG UNCERTAINTY - MD - JUN 2006
IMPLICIT NONE
INTEGER,INTENT(IN)::N_Y,N_yP,N_ZR,N_RL,N_R,N_Z,N_C,N_X,N_A,N_U,MONTH,N_YR,RANK,&
IERR,N_BDRF,MEAS_SWITCH
INTEGER,INTENT(IN)::N_ZP(N_ZR),VEG_MAP(N_yP)
INTEGER,INTENT(INOUT) :: CTRL(N_C)
REAL,INTENT(IN)::Y(N_Y,N_RL),FREQ(CTRL(9)),THETA(CTRL(9)),X(N_X*N_yP,N_RL),&
FORCING(N_U,N_RL,N_yP),ALBEDO(N_yP,N_RL), F_VEG(N_R,1),BDRF(N_BDRF,2),&
VCOVER(N_yP)
REAL,INTENT(INOUT)::ZOUT(1:N_Z,N_RL),TBRAW(N_yP*N_ZP(1),N_RL)
INTEGER R,P,YI(3),ZI(4),i,J,K,L,M,N,O,S,T,KGLOBAL
REAL, DIMENSION(:), ALLOCATABLE :: DZ,ATM_IN,CAN_IN,AUX_SNOW_IN
REAL,DIMENSION(:,:),ALLOCATABLE :: SNOW_IN,ZHAT,TB
INTEGER,INTENT(IN) :: MEAS
INTEGER,PARAMETER :: N_VEG=13 !NUMBER OF SSIB VEGETATION TYPES
INTEGER,PARAMETER :: N_MOS=12 !NUMBER OF MONTHS
!SSIB AND SLA VEGETATION DATA
REAL :: POROSITY(N_VEG), LAI(N_VEG,N_MOS,2),SLA(N_VEG),TRAN_VEG_NIR(N_VEG),&
REF_VEG_NIR(N_VEG)
REAL :: RHO_VEG !USED FOR VEGETATION INPUTS
REAL :: RADIUS_IN,R_SNOW !USED FOR INTERPOLATING ALBEDOS AND INTERPOLATED OUTPUT
REAL :: R_VEG,T_VEG,V_C !USED IN NIR RTM
ALLOCATE( ZHAT(N_yP,2),DZ(3),ATM_IN(CTRL(8)),CAN_IN(CTRL(7)),&
AUX_SNOW_IN(CTRL(5)), TB(2,CTRL(9)) )
! DEFINE CTRL PARAMETERS
CTRL(1)=N_YP !TOTAL NUMBER OF CALCULATIONS TO DO IS NUMBER OF STATE PIXELS
!NOTE: ATM RTM SWITCH CTRL(3) NOW SET IN RUN_PARAMS
!NOTE: CANOPY RTM SWITCH CTRL(4) NOW VARIES WITH DATA
! CALL TO GET SSIB VEGETATION DATA
CALL GET_VEG_DATA(POROSITY,LAI,TRAN_VEG_NIR,REF_VEG_NIR)
! CALL TO GET SPECIFIC LEAF AREA DATA
CALL GET_SLA_DATA(SLA)
! DEFINE VEGETATION DENSITY
RHO_VEG=950.
DO R=1,N_RL
CTRL(2)=3!OBTAIN NUMBER OF SNOW LAYERS FROM AUXILIARY SNOW ARRAY
ZOUT(1:12,R)=0. !INITIALIZE ZOUT SUMS FOR COMPUTING AVERAGE TB
IF(N_R.EQ.1)THEN
KGLOBAL=1 !IF THERE IS ONE GLOBAL REPLICATE (I.E., GENERATOR IS RUNNING)
ELSE
KGLOBAL=RANK*N_RL+R !GLOBAL REPLICATE NUMBER, FOR INDEXING F_VEG FOR FILTER
END IF
DO P=1,N_YP
! 0) ISOLATE VEGETATION COVER FRACTION FOR THIS PIXEL
V_C=VCOVER(P)/100.
! 1) SNOW AND SOIL INPUTS
J=N_YR*(P-1)+1 !BEGINNING Y INDEX
K=N_YR*P !ENDING Y INDEX
L=N_X*(P-1)+1 !BEGINNING X INDEX
M=N_X*P !ENDING X INDEX
N=N_ZP(1)*(P-1)+1 !BEGINNING TB RAW INDEX
O=N_ZP(1)*P !ENDING TB RAW INDEX
! DETERMINE NUMBER OF SNOW LAYERS AND ASSIGN SNOW_IN ARRAY
IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
CTRL(2)=1 !NO SNOW; THIS IS SET TO 1 FOR ALLOCATING SNOW_IN
AUX_SNOW_IN(1)=0 ! NO SNOW
ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
SNOW_IN=0. !FOR NO SNOW CASE, SET THE SNOW_IN ARRAY TO A DUMMY VALUE 0.
ELSEIF(Y(J,R).LT.0.05)THEN
CTRL(2)=1 ! ONE LAYER
AUX_SNOW_IN(1)=1 ! ONE LAYER
ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
SNOW_IN(1,1)=Y(J,R) !SNOW DEPTH
SNOW_IN(1,2)=Y(J+1,R) !ONE LAYER SNOW DENSITY
SNOW_IN(1,3)=Y(J+10,R) !ONE LAYER SNOW GRAIN SIZE
SNOW_IN(1,4)=Y(J+8,R) !ONE LAYER SNOW LIQUID WATER CONTENT
SNOW_IN(1,5)=Y(J+5,R) !ONE LAYER SNOW TEMPERATURE
ELSE
CTRL(2)=3 !THREE SNOW LAYERS
AUX_SNOW_IN(1)=3 ! THREE SNOW LAYERS
ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
CALL SSIB_LAYER(Y(J,R),DZ)
SNOW_IN(1:3,1)=DZ(1:3) !THREE LAYER SNOW DEPTH
SNOW_IN(1:3,2)=Y(J+1:J+3,R) !THREE LAYER SNOW DENSITY
SNOW_IN(1:3,3)=Y(J+10:J+12,R) !THREE LAYER SNOW GRAIN SIZE
SNOW_IN(1:3,4)=Y(J+7:J+9,R) !THREE LAYER SNOW LIQUID WATER CONTENT
SNOW_IN(1:3,5)=Y(J+4:J+6,R) !THREE LAYER SNOW TEMPERATURE
! if(meas.eq.85) then
! print *, 'snow profile for m=85'
! print *, snow_in(1,1:5)
! print *, snow_in(2,1:5)
! print *, snow_in(3,1:5)
! end if
END IF
AUX_SNOW_IN(2)=Y(J+13,R) !SOIL TEMPERATURE
AUX_SNOW_IN(3)=X(L+4,R) !SOIL MOISTURE (USE UPPER LAYER)
AUX_SNOW_IN(4)=POROSITY(VEG_MAP(P))!SOIL POROSITY
AUX_SNOW_IN(5)=0.16 !CONSTANT OF PROPORTIONALITY BETWEEN
!GRAIN SIZE AND CORRELATION LENGTH
! 2) DEFINE CANOPY PM RTM PROPERTIES
IF(LAI(VEG_MAP(P),MONTH,1).LT.1.1E-3.OR.SLA(VEG_MAP(P)).EQ.0.0 &
.OR.SLA(VEG_MAP(P)).EQ.-9999.) THEN
! IF LEAF AREA INDEX IS LESS THAN 0.001, OR SLA IS 0.0, THEN DON'T
! APPLY THE VEGETATION MODEL
CTRL(4)=0
ELSE
CTRL(4)=1 !TURN ON VEGETATION RTM SWITCH
CAN_IN(1)=8 !VEGETATION WATER SALINITY: [PPT]
CAN_IN(2)=X(L+1,R) !CANOPY TEMPERATURE [k]
CAN_IN(3)=0.5 !VEGETATION GRAVIMETRIC WATER CONTENT [FRAC]
! NEEDLE DIAMETER (M) IS CALCULATED FROM SPECIFIC LEAF AREA AND VEGETATION
! DENSITY
CAN_IN(4)=1./SLA(VEG_MAP(P))/RHO_VEG
! VEGETATION BIOMASS (KG VEGETATION /M2 GROUND AREA) IS CALCULATED AS THE
! QUOTIENT OF LEAF AREA INDEX ( M2 VEGETATION AREA / M2 GROUND AREA)
! AND SPECIFIC LEAF AREA (M2 VEGETATION AREA / KG VEGETATION MASS)
! LEAF AREA INDEX PERTURBATION (F_VEG) FACTOR IS APPLIED HERE.
CAN_IN(5)=LAI(VEG_MAP(P),MONTH,1)*EXP(F_VEG(KGLOBAL,1))/SLA(VEG_MAP(P))
END IF
! print *, ctrl(4)
! 3) DEFINE ATMOSPHERIC RTM PROPERTIES
ATM_IN(1)=FORCING(5,R,P)/100.!CONVERT SURFACE PRESSURE TO MBAR FROM pA
ATM_IN(2)=FORCING(6,R,P) !SURFACE ATMOSPHERIC TEMPERATURE
ATM_IN(3)=MONTH !MONTH OF YEAR
ATM_IN(4)=1200*FORCING(7,R,P)!CONVERT MOISTURE TO G/M3 FROM G/G : ASSUME
! AN AIR DENSITY OF 1200 G/M3
! 4) COMPUTE PREDICTED BRIGHTNESS TEMPERATURES AT EACH PIXEL ONLY
! IF MEASUREMENT INTERVAL IS ODD, WHICH CORRESPONDS TO THE 1 am
! MEASUREMENTS. DO NOT COMPUTE IF THE SWITCH INDICATES NIR
! MEASUREMENTS ONLY, OR IF IT INDICATES NO MEASUREMENTS
! THIS OPTION OUGHT TO BE AUTOMATED: DURING THE SYNTHETIC TESTS, I ALWAYS
! CALCULATED ONLY THE NIGHTTIME PM MEASUREMENTS, WHICH SAVED ALOT OF TIME.
! THIS OPTION CAN BE DISABLED BY RE-ARRANGING THE COMMENTS IN THE 'IF' BLOCK
! BELOW, WHICH I DID FOR THE CLPX GBMR WORK.
! IF(MOD(MEAS,2).EQ.0)THEN
! ! DON'T COMPUTE BRIGHTNESS TEMPS DURING THE DAYTIME MEASUREMENTS
! TB=-9999.
! ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
! IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
TB=-9999.
ELSE
CALL RAD_XFER_MOD(CTRL,FREQ,THETA,AUX_SNOW_IN,SNOW_IN,CAN_IN,ATM_IN,&
TB,P,R,RANK,MEAS,V_C,IERR)
END IF
! 5) AGGREGATE PREDICTED OBSERVATIONS OVER THE DOMAIN
!AGGREGATE MEASUREMENT 1 TO MEASUREMENT SCALE
!RECORD TB RAW: SEE JOURNAL ON 10/10/2005 ON INDEXING
DO I=1,CTRL(9)
TBRAW(N+I-1,R)=TB(1,I)
TBRAW(CTRL(9)+N+I-1,R)=TB(2,I)
END DO
! I AM COMMENTING OUT THE PART THAT ONLY CALCULATES TB FOR DAYTIME MEAS...
!ADD HORIZONTAL TB TO ZOUT
DO I=1,CTRL(9)
! IF(MOD(MEAS,2).EQ.0)THEN
! !IT'S DAYTIME
! ZOUT(I,R)=-9999.
! ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
ZOUT(I,R)=-9999.
ELSE
ZOUT(I,R)=ZOUT(I,R)+TB(1,I)/N_YP !THIS WILL BE AN AVERAGE
END If
END DO
!ADD VERTICAL TB TO ZOUT
DO I=1,CTRL(9)
! IF(MOD(MEAS,2).EQ.0)THEN
! !IT'S DAYTIME
! ZOUT(I+6,R)=-9999.
! ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
ZOUT(I+6,R)=-9999.
ELSE
ZOUT(I+6,R)=ZOUT(I+6,R)+TB(2,I)/N_YP !THIS WILL BE AN AVERAGE
END IF
END DO
!ALBEDOS...
! THIS WAS THE OLD BROADBAND MEASUREMENT
! ZOUT(N_ZP(1)+P,R)=ALBEDO(P,R) !ALBEDO FROM SSIB...
!6) COMPUTE NARROWBAND ALBEDOS
IF(MOD(MEAS,2).EQ.0)THEN
!IT'S DAYTTIME
!ISOLATE VEGETATION NIR RTM PROPERTIES FOR THIS PIXEL
T_VEG=TRAN_VEG_NIR(VEG_MAP(P))
R_VEG=REF_VEG_NIR(VEG_MAP(P))
IF ( CTRL(2).EQ.3 )THEN
RADIUS_IN=0.5*SNOW_IN(3,3)*1.0E6
CALL INTERPOLATE(BDRF(1:N_BDRF,1),BDRF(1:N_BDRF,2),RADIUS_IN,&
R_SNOW,N_BDRF)
!FRACTIONAL VEGETATION COVERAGE MEASUREMENT RELATIONSHIPS
ZOUT(N_ZP(1)+P,R)=R_SNOW*(1-V_C)+(R_SNOW*T_VEG**2+R_VEG)*V_C
ELSEIF ( CTRL(2).EQ.1 )THEN
IF( AUX_SNOW_IN(1).EQ.1) THEN
RADIUS_IN=0.5*SNOW_IN(1,3)*1.0E6
CALL INTERPOLATE(BDRF(1:N_BDRF,1),BDRF(1:N_BDRF,2),RADIUS_IN,&
R_SNOW,N_BDRF)
!FRACTIONAL VEGETATION COVERAGE MEASUREMENT RELATIONSHIPS
ZOUT(N_ZP(1)+P,R)=R_SNOW*(1-V_C)+(R_SNOW*T_VEG**2+R_VEG)*V_C
ELSEIF (AUX_SNOW_IN(1).EQ.0) THEN
!if there is no snow, set albedo to 0.1?
ZOUT(N_ZP(1)+P,R)=0.1
ELSE
PRINT *, 'INTERFACEZ: ERROR IN ALBEDO ASSIGNMENT CASE STRUCTURE 1.',&
'ABORTING'
CALL MPI_FINALIZE(IERR)
STOP
END IF
ELSE
PRINT *, 'INTERFACEZ: ERROR IN ALBEDO ASSIGNMENT CASE STRUCTURE 2.',&
'ABORTING'
CALL MPI_FINALIZE(IERR)
STOP
END IF
ELSE
!IT'S NIGHTTIME, SET BDRF TO NODATA
ZOUT(N_ZP(1)+P,R)=-9999.
END IF
!7) SURFACE TEMPERATURE MEASUREMENT
IF(CTRL(4).EQ.1) THEN
!OLD WAY: ZOUT(SUM(N_ZP(1:2))+P,R)=X(L+1,R)
!IF THE VEGETATION RTM SWITCH IS ON CTRL(4)=1), THEN WE MUST TAKE THE
! VEGETATION TEMPERATURE INTO ACCOUNT. APPLY SIMPLE FRACTIONAL
! FRACTIONAL VEGETATION COVERAGE MODEL
IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
!USE GROUND TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*AUX_SNOW_IN(2)
ELSEIF(Y(J,R).LT.0.05)THEN
!USE ONE LAYER SNOW TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*SNOW_IN(1,5)
ELSE
!USE SURFACE SNOW TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*SNOW_IN(3,5)
END IF
ELSEIF(CTRL(4).EQ.0) THEN
!THIS CODE IS NOW REFERENCE BASICALLY IF WE HAVE BARE SOIL
IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
!USE GROUND TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)= AUX_SNOW_IN(2)
ELSEIF(Y(J,R).LT.0.05)THEN
!USE ONE LAYER SNOW TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)= SNOW_IN(1,5)
ELSE
!USE SURFACE SNOW TEMPERATURE
ZOUT(SUM(N_ZP(1:2))+P,R)= SNOW_IN(3,5)
END IF
ELSE
PRINT *, 'ERROR IN INTERFACEZ; ILLEGAL VALUE OF CTRL(4). STOPPING.'
CALL MPI_FINALIZE(IERR)
STOP
END IF
DEALLOCATE(SNOW_IN)
END DO
END DO
deallocate( ZHAT,dz,atm_in,can_in,aux_snow_in,tb)
END SUBROUTINE INTERFACEZ