| blob | 3132b4678faf0b5d2b9a56ebb4883e6b18029048 |
1 MODULE module_soil_pre
3 USE module_date_time
4 USE module_state_description
6 CHARACTER (LEN=3) :: num_cat_count
7 INTEGER , PARAMETER , DIMENSION(0:300) :: ints = &
8 (/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, &
9 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, &
10 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, &
11 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, &
12 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, &
13 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, &
14 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, &
15 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, &
16 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, &
17 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, &
18 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, &
19 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, &
20 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, &
21 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, &
22 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, &
23 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, &
24 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, &
25 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, &
26 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, &
27 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, &
28 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, &
29 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, &
30 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, &
31 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, &
32 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, &
33 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, &
34 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, &
35 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, &
36 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, &
37 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300 /)
39 ! Excluded middle processing
41 LOGICAL , SAVE :: hold_ups
42 INTEGER , SAVE :: em_width
44 LOGICAL , EXTERNAL :: skip_middle_points_t
46 CONTAINS
48 SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
49 xland , landusef , isltyp , soilcat , soilctop , &
50 soilcbot , tmn , &
51 seaice_threshold , &
52 fractional_seaice, &
53 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
54 iswater , isice , &
55 scheme , &
56 ids , ide , jds , jde , kds , kde , &
57 ims , ime , jms , jme , kms , kme , &
58 its , ite , jts , jte , kts , kte )
60 IMPLICIT NONE
62 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
63 ims , ime , jms , jme , kms , kme , &
64 its , ite , jts , jte , kts , kte , &
65 iswater , isice
66 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
68 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
69 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
70 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
71 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
72 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
73 vegcat, xland , soilcat , tmn
74 REAL , INTENT(IN) :: seaice_threshold
76 INTEGER :: i , j , num_seaice_changes , loop
77 CHARACTER (LEN=132) :: message
79 INTEGER, INTENT(IN) :: fractional_seaice
80 REAL :: XICE_THRESHOLD
82 IF ( FRACTIONAL_SEAICE == 0 ) THEN
83 xice_threshold = 0.5
84 ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
85 xice_threshold = 0.02
86 ENDIF
88 num_seaice_changes = 0
89 fix_seaice : SELECT CASE ( scheme )
91 CASE ( SLABSCHEME )
92 DO j = jts , MIN(jde-1,jte)
93 DO i = its , MIN(ide-1,ite)
94 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
95 IF ( xice(i,j) .GT. 200.0 ) THEN
96 xice(i,j) = 0.
97 num_seaice_changes = num_seaice_changes + 1
98 END IF
99 END DO
100 END DO
101 IF ( num_seaice_changes .GT. 0 ) THEN
102 WRITE ( message , FMT='(A,I6)' ) &
103 'Total pre number of sea ice locations removed (due to FLAG values) = ', &
104 num_seaice_changes
105 CALL wrf_debug ( 0 , message )
106 END IF
107 num_seaice_changes = 0
108 DO j = jts , MIN(jde-1,jte)
109 DO i = its , MIN(ide-1,ite)
110 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
111 IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
112 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
113 IF ( FRACTIONAL_SEAICE == 0 ) THEN
114 xice(i,j) = 1.0
115 ENDIF
116 num_seaice_changes = num_seaice_changes + 1
117 if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
118 vegcat(i,j)=isice
119 ivgtyp(i,j)=isice
120 lu_index(i,j)=isice
121 landmask(i,j)=1.
122 xland(i,j)=1.
123 DO loop=1,num_veg_cat
124 landusef(i,loop,j)=0.
125 END DO
126 landusef(i,ivgtyp(i,j),j)=1.
128 isltyp(i,j) = 16
129 soilcat(i,j)=isltyp(i,j)
130 DO loop=1,num_soil_top_cat
131 soilctop(i,loop,j)=0
132 END DO
133 DO loop=1,num_soil_bot_cat
134 soilcbot(i,loop,j)=0
135 END DO
136 soilctop(i,isltyp(i,j),j)=1.
137 soilcbot(i,isltyp(i,j),j)=1.
138 ELSE
139 xice(i,j) = 0.0
140 END IF
141 END DO
142 END DO
143 IF ( num_seaice_changes .GT. 0 ) THEN
144 WRITE ( message , FMT='(A,I6)' ) &
145 'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
146 CALL wrf_debug ( 0 , message )
147 END IF
149 CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME ,CLMSCHEME,SSIBSCHEME) !mchen add for ssib
150 num_seaice_changes = 0
151 DO j = jts , MIN(jde-1,jte)
152 DO i = its , MIN(ide-1,ite)
153 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
154 IF ( landmask(i,j) .GT. 0.5 ) THEN
155 if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
156 xice(i,j) = 0.
157 END IF
158 END DO
159 END DO
160 IF ( num_seaice_changes .GT. 0 ) THEN
161 WRITE ( message , FMT='(A,I6)' ) &
162 'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
163 CALL wrf_debug ( 0 , message )
164 END IF
166 END SELECT fix_seaice
168 END SUBROUTINE adjust_for_seaice_pre
170 SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
171 xland , landusef , isltyp , soilcat , soilctop , &
172 soilcbot , tmn , vegfra , &
173 tslb , smois , sh2o , &
174 seaice_threshold , &
175 sst , flag_sst , &
176 fractional_seaice, &
177 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
178 num_soil_layers , &
179 iswater , isice , &
180 scheme , &
181 ids , ide , jds , jde , kds , kde , &
182 ims , ime , jms , jme , kms , kme , &
183 its , ite , jts , jte , kts , kte )
185 IMPLICIT NONE
187 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
188 ims , ime , jms , jme , kms , kme , &
189 its , ite , jts , jte , kts , kte , &
190 iswater , isice
191 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
192 INTEGER , INTENT(IN) :: num_soil_layers
194 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
195 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
196 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
197 REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
198 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN):: sst
199 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
200 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
201 vegcat, xland , soilcat , tmn , &
202 tsk_old , vegfra
203 INTEGER , INTENT(IN) :: flag_sst
204 REAL , INTENT(IN) :: seaice_threshold
205 REAL :: total_depth , mid_point_depth
207 INTEGER :: i , j , num_seaice_changes , loop
208 CHARACTER (LEN=132) :: message
211 INTEGER, INTENT(IN) :: fractional_seaice
212 real :: xice_threshold
214 IF ( FRACTIONAL_SEAICE == 0 ) THEN
215 xice_threshold = 0.5
216 ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
217 xice_threshold = 0.02
218 ENDIF
220 num_seaice_changes = 0
221 fix_seaice : SELECT CASE ( scheme )
223 CASE ( SLABSCHEME )
225 #ifdef WRF_USE_CTSM
226 CASE ( LSMSCHEME , NOAHMPSCHEME , CLMSCHEME, CTSMSCHEME, SSIBSCHEME ) !mchen add for ssib
227 #else
228 CASE ( LSMSCHEME , NOAHMPSCHEME , CLMSCHEME, SSIBSCHEME ) !mchen add for ssib
229 #endif
231 DO j = jts , MIN(jde-1,jte)
232 DO i = its , MIN(ide-1,ite)
233 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
234 IF ( xice(i,j) .GT. 200.0 ) THEN
235 xice(i,j) = 0.
236 num_seaice_changes = num_seaice_changes + 1
237 END IF
238 END DO
239 END DO
240 IF ( num_seaice_changes .GT. 0 ) THEN
241 WRITE ( message , FMT='(A,I6)' ) &
242 'Total post number of sea ice locations removed (due to FLAG values) = ', &
243 num_seaice_changes
244 CALL wrf_debug ( 0 , message )
245 END IF
246 num_seaice_changes = 0
247 DO j = jts , MIN(jde-1,jte)
248 DO i = its , MIN(ide-1,ite)
249 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
250 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
251 ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
252 tsk(i,j) = tsk_old(i,j)
253 END IF
254 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
255 ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
256 print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
257 CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
258 ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
259 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
260 IF ( FRACTIONAL_SEAICE == 0 ) THEN
261 xice(i,j) = 1.0
262 ENDIF
263 num_seaice_changes = num_seaice_changes + 1
264 if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
265 vegcat(i,j)=isice
266 ivgtyp(i,j)=isice
267 lu_index(i,j)=isice
268 landmask(i,j)=1.
269 xland(i,j)=1.
270 vegfra(i,j)=0.
271 DO loop=1,num_veg_cat
272 landusef(i,loop,j)=0.
273 END DO
274 landusef(i,ivgtyp(i,j),j)=1.
276 tsk_old(i,j) = tsk(i,j)
278 isltyp(i,j) = 16
279 soilcat(i,j)=isltyp(i,j)
280 DO loop=1,num_soil_top_cat
281 soilctop(i,loop,j)=0
282 END DO
283 DO loop=1,num_soil_bot_cat
284 soilcbot(i,loop,j)=0
285 END DO
286 soilctop(i,isltyp(i,j),j)=1.
287 soilcbot(i,isltyp(i,j),j)=1.
289 total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
290 DO loop = 1,num_soil_layers
291 mid_point_depth=(total_depth/num_soil_layers)/2. + &
292 (loop-1)*(total_depth/num_soil_layers)
293 tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
294 mid_point_depth*tmn(i,j) ) / total_depth
295 END DO
297 DO loop=1,num_soil_layers
298 smois(i,loop,j) = 1.0
299 sh2o(i,loop,j) = 0.0
300 END DO
301 ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
302 xice(i,j) = 0.
303 END IF
304 END DO
305 END DO
306 IF ( num_seaice_changes .GT. 0 ) THEN
307 WRITE ( message , FMT='(A,I6)' ) &
308 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
309 CALL wrf_debug ( 0 , message )
310 END IF
312 CASE ( RUCLSMSCHEME )
313 DO j = jts , MIN(jde-1,jte)
314 DO i = its , MIN(ide-1,ite)
315 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
316 IF ( xice(i,j) .GT. 200.0 ) THEN
317 xice(i,j) = 0.
318 num_seaice_changes = num_seaice_changes + 1
319 END IF
320 END DO
321 END DO
322 IF ( num_seaice_changes .GT. 0 ) THEN
323 WRITE ( message , FMT='(A,I6)' ) &
324 'Total post number of sea ice locations removed (due to FLAG values) = ', &
325 num_seaice_changes
326 CALL wrf_debug ( 0 , message )
327 END IF
328 num_seaice_changes = 0
329 DO j = jts , MIN(jde-1,jte)
330 DO i = its , MIN(ide-1,ite)
331 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
332 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
333 ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
334 tsk(i,j) = tsk_old(i,j)
335 END IF
336 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
337 ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
338 print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
339 CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
340 ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
341 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
342 IF ( FRACTIONAL_SEAICE == 0 ) THEN
343 xice(i,j) = 1.0
344 ELSE
345 xice(i,j)=max(0.25,xice(i,j))
346 ENDIF
347 num_seaice_changes = num_seaice_changes + 1
348 if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
349 vegcat(i,j)=isice
350 ivgtyp(i,j)=isice
351 lu_index(i,j)=isice
352 landmask(i,j)=1.
353 xland(i,j)=1.
354 vegfra(i,j)=0.
355 DO loop=1,num_veg_cat
356 landusef(i,loop,j)=0.
357 END DO
358 landusef(i,ivgtyp(i,j),j)=1.
360 !tgs - compute blended sea ice/water skin temperature
361 if(flag_sst.eq.1) then
362 tsk(i,j) = xice(i,j)*(min(271.4,tsk(i,j))) &
363 +(1-xice(i,j))*sst(i,j)
364 else
365 tsk(i,j) = xice(i,j)*(min(271.4,tsk(i,j))) &
366 +(1-xice(i,j))*tsk(i,j)
367 endif
368 tsk_old(i,j) = tsk(i,j)
370 isltyp(i,j) = 16
371 soilcat(i,j)=isltyp(i,j)
372 DO loop=1,num_soil_top_cat
373 soilctop(i,loop,j)=0
374 END DO
375 DO loop=1,num_soil_bot_cat
376 soilcbot(i,loop,j)=0
377 END DO
378 soilctop(i,isltyp(i,j),j)=1.
379 soilcbot(i,isltyp(i,j),j)=1.
381 total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
382 tslb(i,1,j) = tsk(i,j)
383 tslb(i,num_soil_layers,j) = tmn(i,j)
384 DO loop = 2,num_soil_layers-1
385 mid_point_depth=(total_depth/num_soil_layers)/4. + &
386 (loop-2)*(total_depth/num_soil_layers)
387 tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
388 mid_point_depth*tmn(i,j) ) / total_depth
389 END DO
391 DO loop=1,num_soil_layers
392 smois(i,loop,j) = 1.0
393 sh2o(i,loop,j) = 0.0
394 END DO
395 ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
396 xice(i,j) = 0.
397 END IF
398 END DO
399 END DO
400 IF ( num_seaice_changes .GT. 0 ) THEN
401 WRITE ( message , FMT='(A,I6)' ) &
402 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
403 CALL wrf_debug ( 0 , message )
404 END IF
407 END SELECT fix_seaice
409 END SUBROUTINE adjust_for_seaice_post
411 SUBROUTINE process_percent_cat_new ( landmask , &
412 landuse_frac , soil_top_cat , soil_bot_cat , &
413 isltyp , ivgtyp , &
414 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
415 ids , ide , jds , jde , kds , kde , &
416 ims , ime , jms , jme , kms , kme , &
417 its , ite , jts , jte , kts , kte , &
418 iswater )
420 IMPLICIT NONE
422 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
423 ims , ime , jms , jme , kms , kme , &
424 its , ite , jts , jte , kts , kte , &
425 iswater
426 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
427 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
428 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
429 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
430 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
431 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
433 INTEGER :: i , j , l , ll, dominant_index
434 REAL :: dominant_value
436 #if ( WRF_CHEM == 1 )
437 ! REAL :: lwthresh = .99
438 REAL :: lwthresh = .50
439 #else
440 REAL :: lwthresh = .50
441 #endif
443 INTEGER , PARAMETER :: iswater_soil = 14
444 INTEGER :: iforce
445 CHARACTER (LEN=132) :: message
446 CHARACTER(LEN=256) :: mminlu
447 LOGICAL :: aggregate_lu
448 integer :: change_water , change_land
449 change_water = 0
450 change_land = 0
452 ! Sanity check on the 50/50 points
454 DO j = jts , MIN(jde-1,jte)
455 DO i = its , MIN(ide-1,ite)
456 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
457 dominant_value = landuse_frac(i,iswater,j)
458 IF ( dominant_value .EQ. lwthresh ) THEN
459 DO l = 1 , num_veg_cat
460 IF ( l .EQ. iswater ) CYCLE
461 IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
462 PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
463 landuse_frac(i,l,j) = lwthresh - .01
464 landuse_frac(i,iswater,j) = lwthresh + 0.01
465 ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
466 PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
467 landuse_frac(i,l,j) = lwthresh + .01
468 landuse_frac(i,iswater,j) = lwthresh - 0.01
469 END IF
470 END DO
471 END IF
472 END DO
473 END DO
475 ! Compute the aggregate of the vegetation/land use categories. Lump all of the grasses together,
476 ! all of the shrubs, all of the trees, etc. Choose the correct set of available land use
477 ! categories. Also, make sure that the user wants to actually avail themselves of aforementioned
478 ! opportunity, as mayhaps they don't.
480 CALL nl_get_mminlu ( 1 , mminlu )
481 CALL nl_get_aggregate_lu ( 1 , aggregate_lu )
482 IF ( aggregate_lu ) THEN
483 DO j = jts , MIN(jde-1,jte)
484 DO i = its , MIN(ide-1,ite)
485 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
486 CALL aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu(1:4) )
487 END DO
488 END DO
489 END IF
491 ! Compute the dominant VEGETATION INDEX.
493 DO j = jts , MIN(jde-1,jte)
494 DO i = its , MIN(ide-1,ite)
495 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
496 dominant_value = landuse_frac(i,1,j)
497 dominant_index = 1
498 DO l = 2 , num_veg_cat
499 IF ( l .EQ. iswater ) THEN
500 ! wait a bit
501 ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
502 dominant_value = landuse_frac(i,l,j)
503 dominant_index = l
504 END IF
505 END DO
506 IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
507 dominant_value = landuse_frac(i,iswater,j)
508 dominant_index = iswater
509 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
510 ( landmask(i,j) .LT. 0.5) .AND. &
511 ( dominant_value .EQ. lwthresh) ) THEN
512 dominant_value = landuse_frac(i,iswater,j)
513 dominant_index = iswater
514 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
515 ( landmask(i,j) .GT. 0.5) .AND. &
516 ( dominant_value .EQ. lwthresh) ) THEN
517 !no op
518 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
519 ( dominant_value .LT. lwthresh ) ) THEN
520 dominant_value = landuse_frac(i,iswater,j)
521 dominant_index = iswater
522 END IF
523 IF ( dominant_index .EQ. iswater ) THEN
524 if(landmask(i,j).gt.lwthresh) then
525 !print *,'changing to water at point ',i,j
526 !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
527 !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
528 !CALL wrf_debug(1,message)
529 !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
530 !CALL wrf_debug(1,message)
531 change_water=change_water+1
532 endif
533 landmask(i,j) = 0
534 ELSE IF ( dominant_index .NE. iswater ) THEN
535 if(landmask(i,j).lt.lwthresh) then
536 !print *,'changing to land at point ',i,j
537 !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
538 !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
539 !CALL wrf_debug(1,message)
540 !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
541 !CALL wrf_debug(1,message)
542 change_land=change_land+1
543 endif
544 landmask(i,j) = 1
545 END IF
546 ivgtyp(i,j) = dominant_index
547 END DO
548 END DO
550 ! Compute the dominant SOIL TEXTURE INDEX, TOP.
552 iforce = 0
553 DO i = its , MIN(ide-1,ite)
554 DO j = jts , MIN(jde-1,jte)
555 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
556 dominant_value = soil_top_cat(i,1,j)
557 dominant_index = 1
558 IF ( landmask(i,j) .GT. lwthresh ) THEN
559 DO l = 2 , num_soil_top_cat
560 IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
561 dominant_value = soil_top_cat(i,l,j)
562 dominant_index = l
563 END IF
564 END DO
565 IF ( dominant_value .LT. 0.01 ) THEN
566 iforce = iforce + 1
567 WRITE ( message , FMT = '(A,I4,I4)' ) &
568 'based on landuse, changing soil to land at point ',i,j
569 CALL wrf_debug(1,message)
570 WRITE ( num_cat_count , FMT = '(I3)' ) num_soil_top_cat
571 WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) (ints(l),l=1,num_soil_top_cat)
572 CALL wrf_debug(1,message)
573 WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) &
574 ((nint(soil_top_cat(i,ints(l),j)*100)), l=1,num_soil_top_cat)
575 CALL wrf_debug(1,message)
576 dominant_index = 8
577 END IF
578 ELSE
579 dominant_index = iswater_soil
580 END IF
581 isltyp(i,j) = dominant_index
582 END DO
583 END DO
585 if(iforce.ne.0)then
586 WRITE(message,FMT='(A,I4,A,I6)' ) &
587 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
588 (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
589 CALL wrf_debug(0,message)
590 endif
591 print *,'LAND CHANGE = ',change_land
592 print *,'WATER CHANGE = ',change_water
594 END SUBROUTINE process_percent_cat_new
596 SUBROUTINE process_soil_real ( tsk , tmn , tavgsfc, &
597 landmask , sst , ht, toposoil, &
598 st_input , sm_input , sw_input , &
599 st_levels_input , sm_levels_input , sw_levels_input , &
600 zs , dzs , flag_sm_adj, tslb , smois , sh2o , &
601 flag_sst , flag_tavgsfc, flag_soilhgt, &
602 flag_soil_layers, flag_soil_levels, &
603 ids , ide , jds , jde , kds , kde , &
604 ims , ime , jms , jme , kms , kme , &
605 its , ite , jts , jte , kts , kte , &
606 sf_surface_physics , num_soil_layers , real_data_init_type , &
607 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
608 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
610 IMPLICIT NONE
612 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
613 ims , ime , jms , jme , kms , kme , &
614 its , ite , jts , jte , kts , kte , &
615 sf_surface_physics , num_soil_layers , real_data_init_type , &
616 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
617 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
619 INTEGER , INTENT(IN) :: flag_sst, flag_tavgsfc, flag_sm_adj
620 INTEGER , INTENT(IN) :: flag_soil_layers, flag_soil_levels, flag_soilhgt
622 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
624 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
625 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
626 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
627 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
628 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
629 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
631 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
632 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
633 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tavgsfc, ht, toposoil
634 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk, tmn
636 INTEGER :: i , j , k, l , dominant_index , num_soil_cat , num_veg_cat, closest_layer
637 REAL :: dominant_value, closest_depth, diff_cm
638 REAL , ALLOCATABLE , DIMENSION(:) :: depth ! Soil layer thicknesses (cm)
639 REAL, PARAMETER :: get_temp_closest_to = 30. ! use soil temperature closest to this depth (cm)
640 REAL, PARAMETER :: something_big = 1.e6 ! Initialize closest depth as something big (cm)
641 INTEGER :: something_far = 1000 ! Soil array index far away
642 CHARACTER (LEN=132) :: message
644 ! Case statement for tmn initialization
645 ! Need to have a reasonable default value for annual mean deeeeep soil temperature
646 ! For sf_surface_physics = 1, we want to use close to a 30 cm value
647 ! for the bottom level of the soil temps.
648 ! NOTE: We are assuming that soil_layers are the same for each grid point
650 fix_bottom_level_for_temp : SELECT CASE ( sf_surface_physics )
651 CASE (SLABSCHEME)
652 IF ( flag_tavgsfc .EQ. 1 ) THEN
653 CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
654 DO j = jts , MIN(jde-1,jte)
655 DO i = its , MIN(ide-1,ite)
656 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
657 tmn(i,j) = tavgsfc(i,j)
658 END DO
659 END DO
660 ELSE
661 ! Look for soil temp close to 30 cm
662 closest_layer = something_far
663 closest_depth = something_big
664 DO k = 1, num_st_levels_input
665 diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
666 IF ( diff_cm < closest_depth ) THEN
667 closest_depth = diff_cm
668 closest_layer = k
669 END IF
670 END DO
671 IF ( closest_layer == something_far ) THEN
672 CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
673 CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
674 ELSE
675 write(message, FMT='(A,F7.2,A,I3)')&
676 'Soil temperature closest to ',get_temp_closest_to, &
677 ' at level ',st_levels_input(closest_layer)
678 CALL wrf_debug ( 0 , message )
679 DO j = jts , MIN(jde-1,jte)
680 DO i = its , MIN(ide-1,ite)
681 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
682 tmn(i,j) = st_input(i,closest_layer+1,j)
683 END DO
684 END DO
685 END IF
686 END IF
688 CASE (LSMSCHEME)
690 CASE (NOAHMPSCHEME)
692 CASE (CLMSCHEME)
694 CASE (RUCLSMSCHEME)
696 CASE (PXLSMSCHEME)
697 ! When the input data from met_em is in layers, there is an additional level added to the beginning
698 ! of the array to define the surface, which is why we add the extra value (closest_layer+1)
699 IF ( flag_tavgsfc .EQ. 1 ) THEN
700 CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
701 DO j = jts , MIN(jde-1,jte)
702 DO i = its , MIN(ide-1,ite)
703 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
704 tmn(i,j) = tavgsfc(i,j)
705 END DO
706 END DO
707 ELSE
708 ! Look for soil temp close to 30 cm
709 closest_layer = num_st_levels_input+1
710 closest_depth = something_big
711 DO k = 1, num_st_levels_input
712 diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
713 IF ( diff_cm < closest_depth ) THEN
714 closest_depth = diff_cm
715 closest_layer = k
716 END IF
717 END DO
718 IF ( closest_layer == num_st_levels_input + 1 ) THEN
719 CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
720 CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
721 ELSE
722 write(message, FMT='(A,F7.2,A,I3)')&
723 'Soil temperature closest to ',get_temp_closest_to, &
724 ' at level ',st_levels_input(closest_layer)
725 CALL wrf_debug ( 0 , message )
726 DO j = jts , MIN(jde-1,jte)
727 DO i = its , MIN(ide-1,ite)
728 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
729 tmn(i,j) = st_input(i,closest_layer+1,j)
730 END DO
731 END DO
732 END IF
733 END IF
735 #if 0
736 ! Loop over layers and do a weighted mean
737 IF ( ALLOCATED ( depth ) ) DEALLOCATE ( depth )
738 ALLOCATE ( depth(num_st_levels_input) )
739 IF ( flag_soil_layers == 1 ) THEN
740 DO k = num_st_levels_input, 2, -1
741 depth(k) = st_levels_input(k) - st_levels_input(k-1)
742 END DO
743 depth(1) = st_levels_input(1)
744 ELSE IF ( flag_soil_levels == 1 ) THEN
745 DO k = 2, num_st_levels_input
746 depth(k) = st_levels_input(k) - st_levels_input(k-1)
747 END DO
748 depth(1) = 0.
749 END IF
751 DO j = jts , MIN(jde-1,jte)
752 DO i = its , MIN(ide-1,ite)
753 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
754 tmn(i,j) = 0.
755 DO k = 1, num_st_levels_input
756 tmn(i,j) = tmn(i,j) + depth(k) * st_input(i,k,j)
757 END DO
758 END DO
759 END DO
760 DEALLOCATE ( depth )
762 #endif
764 END SELECT fix_bottom_level_for_temp
766 ! Adjust the various soil temperature values depending on the difference in
767 ! elevation between the current model's elevation and the incoming data's
768 ! orography.
770 adjust_soil : SELECT CASE ( sf_surface_physics )
772 CASE ( SLABSCHEME,LSMSCHEME,NOAHMPSCHEME,RUCLSMSCHEME,PXLSMSCHEME,CLMSCHEME,SSIBSCHEME )
773 CALL adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ht , &
774 toposoil , landmask , st_input, st_levels_input, &
775 flag_soilhgt , flag_tavgsfc , &
776 flag_soil_layers , flag_soil_levels, &
777 num_st_levels_input, num_st_levels_alloc, &
778 ids , ide , jds , jde , kds , kde , &
779 ims , ime , jms , jme , kms , kme , &
780 its , ite , jts , jte , kts , kte )
782 END SELECT adjust_soil
784 ! Initialize the soil depth, and the soil temperature and moisture.
786 IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
787 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
788 CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
789 landmask , sst , flag_sst , &
790 ids , ide , jds , jde , kds , kde , &
791 ims , ime , jms , jme , kms , kme , &
792 its , ite , jts , jte , kts , kte )
794 ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
795 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
796 CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
797 st_input , sm_input , sw_input , landmask , sst , &
798 zs , dzs , &
799 st_levels_input , sm_levels_input , sw_levels_input , &
800 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
801 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
802 flag_sst , flag_soil_layers , flag_soil_levels , &
803 ids , ide , jds , jde , kds , kde , &
804 ims , ime , jms , jme , kms , kme , &
805 its , ite , jts , jte , kts , kte )
806 ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
807 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
808 CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
809 st_input , sm_input , sw_input , landmask , sst , &
810 zs , dzs , &
811 st_levels_input , sm_levels_input , sw_levels_input , &
812 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
813 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
814 flag_sst , flag_soil_layers , flag_soil_levels , &
815 ids , ide , jds , jde , kds , kde , &
816 ims , ime , jms , jme , kms , kme , &
817 its , ite , jts , jte , kts , kte )
818 ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
819 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
820 CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
821 st_input , sm_input , landmask , sst , &
822 zs , dzs , flag_sm_adj, &
823 st_levels_input , sm_levels_input , &
824 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
825 num_st_levels_alloc , num_sm_levels_alloc , &
826 flag_sst , flag_soil_layers , flag_soil_levels , &
827 ids , ide , jds , jde , kds , kde , &
828 ims , ime , jms , jme , kms , kme , &
829 its , ite , jts , jte , kts , kte )
830 !CLM -- Jiming Jin 10/17/2012
831 ELSE IF ( ( sf_surface_physics .EQ. CLMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
832 CALL init_soil_depth_4 ( zs , dzs , num_soil_layers )
833 CALL init_soil_4_real ( tsk , tmn , smois , sh2o , tslb , &
834 st_input , sm_input , sw_input , landmask , sst , &
835 zs , dzs , &
836 st_levels_input , sm_levels_input , sw_levels_input , &
837 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
838 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
839 flag_sst , flag_soil_layers , flag_soil_levels , &
840 ids , ide , jds , jde , kds , kde , &
841 ims , ime , jms , jme , kms , kme , &
842 its , ite , jts , jte , kts , kte )
843 ELSE IF ( ( sf_surface_physics .EQ. PXLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
844 CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
845 CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
846 st_input , sm_input , sw_input, landmask , sst , &
847 zs , dzs , &
848 st_levels_input , sm_levels_input , sw_levels_input, &
849 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
850 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
851 flag_sst , flag_soil_layers , flag_soil_levels , &
852 ids , ide , jds , jde , kds , kde , &
853 ims , ime , jms , jme , kms , kme , &
854 its , ite , jts , jte , kts , kte )
855 ELSE IF ( ( sf_surface_physics .EQ. SSIBSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
856 CALL init_soil_depth_8 ( zs , dzs , num_soil_layers )
857 CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
858 st_input , sm_input , sw_input , landmask , sst , &
859 zs , dzs , &
860 st_levels_input , sm_levels_input , sw_levels_input , &
861 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
862 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
863 flag_sst , flag_soil_layers , flag_soil_levels , &
864 ids , ide , jds , jde , kds , kde , &
865 ims , ime , jms , jme , kms , kme , &
866 its , ite , jts , jte , kts , kte )
867 END IF
869 END SUBROUTINE process_soil_real
871 SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
872 ivgtyp,isltyp,tslb,smois, &
873 tsk,tmn,zs,dzs, &
874 num_soil_layers, &
875 sf_surface_physics , &
876 ids,ide, jds,jde, kds,kde,&
877 ims,ime, jms,jme, kms,kme,&
878 its,ite, jts,jte, kts,kte )
880 IMPLICIT NONE
882 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
883 ims,ime, jms,jme, kms,kme, &
884 its,ite, jts,jte, kts,kte
886 INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
888 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
890 REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
892 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
893 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
894 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
896 ! Local variables.
898 INTEGER :: itf,jtf
900 itf=MIN(ite,ide-1)
901 jtf=MIN(jte,jde-1)
903 IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
904 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
905 CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
906 ivgtyp,zs,dzs,num_soil_layers, &
907 ids,ide, jds,jde, kds,kde, &
908 ims,ime, jms,jme, kms,kme, &
909 its,ite, jts,jte, kts,kte )
910 ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
911 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
912 CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
913 ivgtyp,isltyp,tslb,smois,tmn, &
914 num_soil_layers, &
915 ids,ide, jds,jde, kds,kde, &
916 ims,ime, jms,jme, kms,kme, &
917 its,ite, jts,jte, kts,kte )
918 ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
919 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
920 CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
921 ivgtyp,isltyp,tslb,smois,tmn, &
922 num_soil_layers, &
923 ids,ide, jds,jde, kds,kde, &
924 ims,ime, jms,jme, kms,kme, &
925 its,ite, jts,jte, kts,kte )
926 ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
927 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
929 END IF
931 END SUBROUTINE process_soil_ideal
933 SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ter , &
934 toposoil , landmask , st_input , st_levels_input, &
935 flag_toposoil , flag_tavgsfc , &
936 flag_soil_layers , flag_soil_levels, &
937 num_st_levels_input, num_st_levels_alloc, &
938 ids , ide , jds , jde , kds , kde , &
939 ims , ime , jms , jme , kms , kme , &
940 its , ite , jts , jte , kts , kte )
942 IMPLICIT NONE
944 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
945 ims , ime , jms , jme , kms , kme , &
946 its , ite , jts , jte , kts , kte
947 INTEGER , INTENT(IN) :: num_st_levels_input, num_st_levels_alloc
949 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
950 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
951 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
952 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(IN) :: st_levels_input
954 INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
955 INTEGER , INTENT(IN) :: flag_soil_layers , flag_soil_levels
957 INTEGER :: i , j, k , st_near_sfc
959 REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
961 ! Adjust the annual mean temperature as if it is based on from a sea-level elevation
962 ! if the value used is from the actual annula mean data set. If the input field to
963 ! be used for tmn is one of the first-guess input temp fields, need to do an adjustment
964 ! only on the diff in topo from the model terrain and the first-guess terrain.
966 SELECT CASE ( sf_surface_physics )
968 CASE ( LSMSCHEME , NOAHMPSCHEME,CLMSCHEME )
969 DO j = jts , MIN(jde-1,jte)
970 DO i = its , MIN(ide-1,ite)
971 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
972 IF (landmask(i,j) .GT. 0.5 ) THEN
973 tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
974 END IF
975 END DO
976 END DO
978 CASE (RUCLSMSCHEME)
979 DO j = jts , MIN(jde-1,jte)
980 DO i = its , MIN(ide-1,ite)
981 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
982 IF (landmask(i,j) .GT. 0.5 ) THEN
983 tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
984 END IF
985 END DO
986 END DO
988 END SELECT
991 ! Do we have a soil field with which to modify soil temperatures?
993 IF ( flag_toposoil .EQ. 1 ) THEN
995 DO j = jts , MIN(jde-1,jte)
996 DO i = its , MIN(ide-1,ite)
998 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
999 ! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
1000 ! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
1001 ! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
1002 ! the difference between the soil elevation and the terrain is greater than 3 km means
1003 ! that the soil data is either all zeros or that the data are inconsistent. Any of these
1004 ! three conditions is grievous enough to induce a WRF fatality. However, if they are at
1005 ! a water point, then we can safely ignore them.
1007 soil_elev_min_val = toposoil(i,j)
1008 soil_elev_max_val = toposoil(i,j)
1009 soil_elev_min_dif = ter(i,j) - toposoil(i,j)
1010 soil_elev_max_dif = ter(i,j) - toposoil(i,j)
1012 IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
1013 CYCLE
1014 ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
1015 !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
1016 cycle
1017 ! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
1018 ENDIF
1020 IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
1021 CYCLE
1022 ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
1023 print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
1024 cycle
1025 CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
1026 ENDIF
1028 IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
1029 ( landmask(i,j) .LT. 0.5 ) ) THEN
1030 CYCLE
1031 ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
1032 ( landmask(i,j) .GT. 0.5 ) ) THEN
1033 print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
1034 cycle
1035 CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
1036 ENDIF
1038 ! For each of the fields that we would like to modify, check to see if it came in from the SI.
1039 ! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
1040 ! are not at a water point.
1042 IF (landmask(i,j) .GT. 0.5 ) THEN
1043 IF ( sf_surface_physics .EQ. SLABSCHEME ) THEN
1044 st_near_sfc = 0 ! Check if there is a soil layer above 40 cm
1045 DO k = 1, num_st_levels_input
1046 IF ( st_levels_input(k) .LE. 40 ) THEN
1047 st_near_sfc = 1
1048 END IF
1049 END DO
1050 IF ( ( flag_tavgsfc == 1 ) .OR. ( st_near_sfc == 1 ) ) THEN
1051 tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1052 ELSE
1053 tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
1054 END IF
1055 END IF
1057 tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1059 IF ( flag_soil_layers == 1 ) THEN
1060 DO k = 2, num_st_levels_input+1
1061 st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1062 END DO
1063 ELSE
1064 DO k = 1, num_st_levels_input
1065 st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1066 END DO
1067 END IF
1069 END IF
1070 END DO
1071 END DO
1073 END IF
1075 END SUBROUTINE adjust_soil_temp_new
1078 SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
1080 IMPLICIT NONE
1082 INTEGER, INTENT(IN) :: num_soil_layers
1084 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1086 INTEGER :: l
1088 ! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
1089 ! The distance from the ground level to the midpoint of the layer is given by zs.
1091 ! ------- Ground Level ---------- || || || ||
1092 ! || || || || zs(1) = 0.005 m
1093 ! -- -- -- -- -- -- -- -- -- || || || \/
1094 ! || || ||
1095 ! ----------------------------------- || || || \/ dzs(1) = 0.01 m
1096 ! || || ||
1097 ! || || || zs(2) = 0.02
1098 ! -- -- -- -- -- -- -- -- -- || || \/
1099 ! || ||
1100 ! || ||
1101 ! ----------------------------------- || || \/ dzs(2) = 0.02 m
1102 ! || ||
1103 ! || ||
1104 ! || ||
1105 ! || || zs(3) = 0.05
1106 ! -- -- -- -- -- -- -- -- -- || \/
1107 ! ||
1108 ! ||
1109 ! ||
1110 ! ||
1111 ! ----------------------------------- \/ dzs(3) = 0.04 m
1113 IF ( num_soil_layers .NE. 5 ) THEN
1114 PRINT '(A)','The SLAB thermal diffusion LSM uses 5 layers. Set num_soil_layers=5 in the namelist (&physics).'
1115 CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
1116 END IF
1118 dzs(1)=.01
1119 zs(1)=.5*dzs(1)
1121 DO l=2,num_soil_layers
1122 dzs(l)=2*dzs(l-1)
1123 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
1124 ENDDO
1126 END SUBROUTINE init_soil_depth_1
1128 SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
1130 IMPLICIT NONE
1132 INTEGER, INTENT(IN) :: num_soil_layers
1134 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1136 INTEGER :: l
1138 dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
1140 IF ( num_soil_layers .NE. 4 ) THEN
1141 PRINT '(A)','The Noah and NoahMP LSMs use 4 layers. Set num_soil_layers=4 in the namelist (&physics).'
1142 CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
1143 END IF
1145 zs(1)=.5*dzs(1)
1147 DO l=2,num_soil_layers
1148 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
1149 ENDDO
1151 END SUBROUTINE init_soil_depth_2
1153 SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
1155 IMPLICIT NONE
1157 INTEGER, INTENT(IN) :: num_soil_layers
1159 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1160 REAL, DIMENSION(1:num_soil_layers) :: zs2
1162 INTEGER :: l
1164 CHARACTER (LEN=132) :: message
1166 ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
1167 ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
1168 ! Other options with number of levels are possible, but
1169 ! WRF users should change consistently the namelist entry with the
1170 ! ZS array in this subroutine.
1172 IF ( num_soil_layers .EQ. 6) THEN
1173 zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
1174 ELSEIF ( num_soil_layers .EQ. 9) THEN
1175 zs = (/ 0.00 , 0.01 , 0.04 , 0.10 , 0.30, 0.60, 1.00 , 1.60, 3.00 /)
1176 ! zs = (/ 0.00 , 0.005 , 0.02 , 0.10 , 0.30, 0.60, 1.00 , 1.60, 3.00 /)
1177 ENDIF
1179 zs2(1) = 0.
1180 zs2(2) = (zs(2) + zs(1))*0.5
1181 dzs(1) = zs2(2) - zs2(1)
1182 do l = 2, num_soil_layers - 1
1183 zs2(l) = (zs(l+1) + zs(l)) * 0.5
1184 dzs(l) = zs2(l) - zs2(l-1)
1185 enddo
1186 zs2(num_soil_layers) = zs(num_soil_layers)
1187 dzs(num_soil_layers) = zs2(num_soil_layers) - zs2(num_soil_layers-1)
1189 IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
1190 write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels. Set num_soil_layers in the namelist (&physics).'
1191 CALL wrf_error_fatal ( message )
1192 END IF
1194 END SUBROUTINE init_soil_depth_3
1196 SUBROUTINE init_soil_depth_4 ( zs , dzs , num_soil_layers )
1198 IMPLICIT NONE
1200 INTEGER, INTENT(IN) :: num_soil_layers
1202 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1203 REAL, PARAMETER :: scalez = 0.025
1205 INTEGER :: l
1207 ! Define layers (top layer = 0.01 m). Double the thicknesses at each
1208 ! step (dzs values).
1209 ! The distance from the ground level to the midpoint of the layer is
1210 ! given by zs.
1212 ! ------- Ground Level ---------- || || || ||
1213 ! || || || || zs(1) =
1214 ! 0.005 m
1215 ! -- -- -- -- -- -- -- -- -- || || || \/
1216 ! || || ||
1217 ! ----------------------------------- || || || \/ dzs(1) =
1218 ! 0.01 m
1219 ! || || ||
1220 ! || || || zs(2) = 0.02
1221 ! -- -- -- -- -- -- -- -- -- || || \/
1222 ! || ||
1223 ! || ||
1224 ! ----------------------------------- || || \/ dzs(2) = 0.02 m
1225 ! || ||
1226 ! || ||
1227 ! || ||
1228 ! || || zs(3) = 0.05
1229 ! -- -- -- -- -- -- -- -- -- || \/
1230 ! ||
1231 ! ||
1232 ! ||
1233 ! ||
1234 ! ----------------------------------- \/ dzs(3) = 0.04 m
1236 do l = 1, num_soil_layers
1237 zs(l) = scalez*(exp(0.5*(l-0.5))-1.) !node depths
1238 enddo
1239 dzs(1) = 0.5*(zs(1)+zs(2)) !thickness b/n two interfaces
1240 do l = 2,num_soil_layers-1
1241 dzs(l)= 0.5*(zs(l+1)-zs(l-1))
1242 enddo
1243 dzs(num_soil_layers) = zs(num_soil_layers)-zs(num_soil_layers-1)
1245 END SUBROUTINE init_soil_depth_4
1247 SUBROUTINE init_soil_depth_7 ( zs , dzs , num_soil_layers )
1249 IMPLICIT NONE
1251 INTEGER, INTENT(IN) :: num_soil_layers
1253 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1255 INTEGER :: l
1257 dzs = (/ 0.01 , 0.99 /)
1259 IF ( num_soil_layers .NE. 2 ) THEN
1260 PRINT '(A)','The PX LSM uses 2 layers. Set num_soil_layers=2 in the namelist (&physics).'
1261 CALL wrf_error_fatal ( 'PXLSM_uses_2_layers' )
1262 END IF
1264 zs(1) = 0.5 * dzs(1)
1265 zs(2) = dzs(1) + 0.5 * dzs(2)
1267 END SUBROUTINE init_soil_depth_7
1269 SUBROUTINE init_soil_depth_8 ( zs , dzs , num_soil_layers )
1271 IMPLICIT NONE
1273 INTEGER, INTENT(IN) :: num_soil_layers
1275 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1277 zs(1) = 0.00
1278 zs(2) = 0.05
1279 zs(3) = 1.05
1281 dzs(1) = 0.05
1282 dzs(2) = 1.00
1283 dzs(3) = 0.00
1285 END SUBROUTINE init_soil_depth_8
1287 SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
1288 num_soil_layers , real_data_init_type , &
1289 landmask , sst , flag_sst , &
1290 ids , ide , jds , jde , kds , kde , &
1291 ims , ime , jms , jme , kms , kme , &
1292 its , ite , jts , jte , kts , kte )
1294 IMPLICIT NONE
1296 INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
1297 ids , ide , jds , jde , kds , kde , &
1298 ims , ime , jms , jme , kms , kme , &
1299 its , ite , jts , jte , kts , kte
1301 INTEGER , INTENT(IN) :: flag_sst
1303 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1304 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
1306 REAL , DIMENSION(num_soil_layers) :: zs , dzs
1308 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
1310 INTEGER :: i , j , l
1312 ! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
1313 ! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
1314 ! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
1315 ! annual mean temperature.
1317 DO j = jts , MIN(jde-1,jte)
1318 DO i = its , MIN(ide-1,ite)
1319 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1320 IF ( landmask(i,j) .GT. 0.5 ) THEN
1321 DO l = 1 , num_soil_layers
1322 tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
1323 tmn(i,j) * ( zs( l) - zs(1) ) ) / &
1324 ( zs(num_soil_layers) - zs(1) )
1325 END DO
1326 ELSE
1327 IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
1328 DO l = 1 , num_soil_layers
1329 tslb(i,l,j)= sst(i,j)
1330 END DO
1331 ELSE
1332 DO l = 1 , num_soil_layers
1333 tslb(i,l,j)= tsk(i,j)
1334 END DO
1335 END IF
1336 END IF
1337 END DO
1338 END DO
1340 END SUBROUTINE init_soil_1_real
1342 SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
1343 ivgtyp,ZS,DZS,num_soil_layers, &
1344 ids,ide, jds,jde, kds,kde, &
1345 ims,ime, jms,jme, kms,kme, &
1346 its,ite, jts,jte, kts,kte )
1348 IMPLICIT NONE
1350 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1351 ims,ime, jms,jme, kms,kme, &
1352 its,ite, jts,jte, kts,kte
1354 INTEGER, INTENT(IN ) :: num_soil_layers
1356 REAL, DIMENSION( ims:ime , num_soil_layers , jms:jme ), INTENT(OUT) :: tslb
1357 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
1358 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
1360 REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
1362 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
1364 ! Lcal variables.
1366 INTEGER :: l,j,i,itf,jtf
1368 itf=MIN(ite,ide-1)
1369 jtf=MIN(jte,jde-1)
1371 IF (num_soil_layers.NE.1)THEN
1372 DO j=jts,jtf
1373 DO l=1,num_soil_layers
1374 DO i=its,itf
1375 tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
1376 ( zs(num_soil_layers)-zs(1) )
1377 ENDDO
1378 ENDDO
1379 ENDDO
1380 ENDIF
1382 ! DO j=jts,jtf
1383 ! DO i=its,itf
1384 ! xland(i,j) = 2
1385 ! ivgtyp(i,j) = 7
1386 ! ENDDO
1387 ! ENDDO
1389 END SUBROUTINE init_soil_1_ideal
1391 SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
1392 st_input , sm_input , sw_input , landmask , sst , &
1393 zs , dzs , &
1394 st_levels_input , sm_levels_input , sw_levels_input , &
1395 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1396 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1397 flag_sst , flag_soil_layers , flag_soil_levels , &
1398 ids , ide , jds , jde , kds , kde , &
1399 ims , ime , jms , jme , kms , kme , &
1400 its , ite , jts , jte , kts , kte )
1402 IMPLICIT NONE
1404 INTEGER , INTENT(IN) :: num_soil_layers , &
1405 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1406 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1407 ids , ide , jds , jde , kds , kde , &
1408 ims , ime , jms , jme , kms , kme , &
1409 its , ite , jts , jte , kts , kte
1411 INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
1413 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1414 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1415 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1417 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1418 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1419 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1420 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1422 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1423 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1424 REAL , DIMENSION(num_soil_layers) :: zs , dzs
1426 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1428 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1430 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1431 REAL :: temp
1432 LOGICAL :: found_levels
1434 CHARACTER (LEN=132) :: message
1436 ! Are there any soil temp and moisture levels - ya know, they are mandatory.
1438 num = num_st_levels_input * num_sm_levels_input
1440 IF ( num .GE. 1 ) THEN
1442 ! Ordered levels that we have data for.
1444 !tgs add option to initialize from RUCLSM
1445 IF ( flag_soil_levels == 1 ) THEN
1446 write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1447 CALL wrf_message ( message )
1448 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
1449 ELSE
1450 write(message, FMT='(A)') ' Assume Noah LSM input'
1451 CALL wrf_message ( message )
1452 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1453 END IF
1456 ! Sort the levels for temperature.
1458 outert : DO lout = 1 , num_st_levels_input-1
1459 innert : DO lin = lout+1 , num_st_levels_input
1460 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1461 temp = st_levels_input(lout)
1462 st_levels_input(lout) = st_levels_input(lin)
1463 st_levels_input(lin) = NINT(temp)
1464 DO j = jts , MIN(jde-1,jte)
1465 DO i = its , MIN(ide-1,ite)
1466 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1467 temp = st_input(i,lout+1,j)
1468 st_input(i,lout+1,j) = st_input(i,lin+1,j)
1469 st_input(i,lin+1,j) = temp
1470 END DO
1471 END DO
1472 END IF
1473 END DO innert
1474 END DO outert
1475 !tgs add IF
1476 IF ( flag_soil_layers == 1 ) THEN
1477 DO j = jts , MIN(jde-1,jte)
1478 DO i = its , MIN(ide-1,ite)
1479 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1480 st_input(i,1,j) = tsk(i,j)
1481 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1482 END DO
1483 END DO
1484 ENDIF
1486 ! Sort the levels for moisture.
1488 outerm: DO lout = 1 , num_sm_levels_input-1
1489 innerm : DO lin = lout+1 , num_sm_levels_input
1490 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1491 temp = sm_levels_input(lout)
1492 sm_levels_input(lout) = sm_levels_input(lin)
1493 sm_levels_input(lin) = NINT(temp)
1494 DO j = jts , MIN(jde-1,jte)
1495 DO i = its , MIN(ide-1,ite)
1496 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1497 temp = sm_input(i,lout+1,j)
1498 sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1499 sm_input(i,lin+1,j) = temp
1500 END DO
1501 END DO
1502 END IF
1503 END DO innerm
1504 END DO outerm
1505 !tgs add IF
1506 IF ( flag_soil_layers == 1 ) THEN
1507 DO j = jts , MIN(jde-1,jte)
1508 DO i = its , MIN(ide-1,ite)
1509 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1510 sm_input(i,1,j) = sm_input(i,2,j)
1511 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1512 END DO
1513 END DO
1514 ENDIF
1516 ! Sort the levels for liquid moisture.
1518 outerw: DO lout = 1 , num_sw_levels_input-1
1519 innerw : DO lin = lout+1 , num_sw_levels_input
1520 IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1521 temp = sw_levels_input(lout)
1522 sw_levels_input(lout) = sw_levels_input(lin)
1523 sw_levels_input(lin) = NINT(temp)
1524 DO j = jts , MIN(jde-1,jte)
1525 DO i = its , MIN(ide-1,ite)
1526 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1527 temp = sw_input(i,lout+1,j)
1528 sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1529 sw_input(i,lin+1,j) = temp
1530 END DO
1531 END DO
1532 END IF
1533 END DO innerw
1534 END DO outerw
1535 IF ( num_sw_levels_input .GT. 1 ) THEN
1536 DO j = jts , MIN(jde-1,jte)
1537 DO i = its , MIN(ide-1,ite)
1538 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1539 sw_input(i,1,j) = sw_input(i,2,j)
1540 sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1541 END DO
1542 END DO
1543 END IF
1545 found_levels = .TRUE.
1547 ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
1549 found_levels = .FALSE.
1551 ELSE
1552 CALL wrf_error_fatal ( &
1553 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
1554 END IF
1556 ! Is it OK to continue?
1558 IF ( found_levels ) THEN
1560 !tgs: Here are the levels that we have from the input for temperature.
1562 IF ( flag_soil_levels == 1 ) THEN
1563 DO l = 1 , num_st_levels_input
1564 zhave(l) = st_levels_input(l) / 100.
1565 END DO
1567 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1569 z_wantt : DO lwant = 1 , num_soil_layers
1570 z_havet : DO lhave = 1 , num_st_levels_input -1
1571 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1572 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1573 DO j = jts , MIN(jde-1,jte)
1574 DO i = its , MIN(ide-1,ite)
1575 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1576 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1577 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1578 ( zhave(lhave+1) - zhave(lhave) )
1579 END DO
1580 END DO
1581 EXIT z_havet
1582 END IF
1583 END DO z_havet
1584 END DO z_wantt
1586 ELSE
1588 ! Here are the levels that we have from the input for temperature. The input levels plus
1589 ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1591 zhave(1) = 0.
1592 DO l = 1 , num_st_levels_input
1593 zhave(l+1) = st_levels_input(l) / 100.
1594 END DO
1595 zhave(num_st_levels_input+2) = 300. / 100.
1597 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1599 z_wantt_2: DO lwant = 1 , num_soil_layers
1600 z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
1601 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1602 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1603 DO j = jts , MIN(jde-1,jte)
1604 DO i = its , MIN(ide-1,ite)
1605 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1606 tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1607 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1608 ( zhave(lhave+1) - zhave(lhave) )
1609 END DO
1610 END DO
1611 EXIT z_havet_2
1612 END IF
1613 END DO z_havet_2
1614 END DO z_wantt_2
1616 END IF
1618 !tgs:
1619 IF ( flag_soil_levels == 1 ) THEN
1620 DO l = 1 , num_sm_levels_input
1621 zhave(l) = sm_levels_input(l) / 100.
1622 END DO
1624 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1626 z_wantm : DO lwant = 1 , num_soil_layers
1627 z_havem : DO lhave = 1 , num_sm_levels_input -1
1628 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1629 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1630 DO j = jts , MIN(jde-1,jte)
1631 DO i = its , MIN(ide-1,ite)
1632 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1633 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1634 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1635 ( zhave(lhave+1) - zhave(lhave) )
1636 END DO
1637 END DO
1638 EXIT z_havem
1639 END IF
1640 END DO z_havem
1641 END DO z_wantm
1643 ELSE
1644 ! Here are the levels that we have from the input for moisture. The input levels plus
1645 ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
1646 ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
1647 ! as the most deep layer's value.
1649 zhave(1) = 0.
1650 DO l = 1 , num_sm_levels_input
1651 zhave(l+1) = sm_levels_input(l) / 100.
1652 END DO
1653 zhave(num_sm_levels_input+2) = 300. / 100.
1655 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1657 z_wantm_2 : DO lwant = 1 , num_soil_layers
1658 z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
1659 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1660 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1661 DO j = jts , MIN(jde-1,jte)
1662 DO i = its , MIN(ide-1,ite)
1663 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1664 smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1665 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1666 ( zhave(lhave+1) - zhave(lhave) )
1667 END DO
1668 END DO
1669 EXIT z_havem_2
1670 END IF
1671 END DO z_havem_2
1672 END DO z_wantm_2
1673 ENDIF
1675 ! Any liquid soil moisture to worry about?
1677 IF ( num_sw_levels_input .GT. 1 ) THEN
1678 ! Interpolate between the layers we have (zhave) and those that we want
1679 ! (zs).
1680 IF ( flag_soil_levels == 1 ) THEN
1681 !tgs: for RUC LSM
1682 z_wantw : DO lwant = 1 , num_soil_layers
1683 z_havew : DO lhave = 1 , num_sw_levels_input -1
1684 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1685 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1686 DO j = jts , MIN(jde-1,jte)
1687 DO i = its , MIN(ide-1,ite)
1688 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1689 sh2o(i,lwant,j)= ( sw_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1690 sw_input(i,lhave+1,j) * ( zs (lwant) - zhave(lhave) ) ) / &
1691 ( zhave(lhave+1) - zhave(lhave) )
1692 END DO
1693 END DO
1694 EXIT z_havew
1695 END IF
1696 END DO z_havew
1697 END DO z_wantw
1699 ELSE
1700 zhave(1) = 0.
1701 DO l = 1 , num_sw_levels_input
1702 zhave(l+1) = sw_levels_input(l) / 100.
1703 END DO
1704 zhave(num_sw_levels_input+2) = 300. / 100.
1706 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1708 z_wantw_2 : DO lwant = 1 , num_soil_layers
1709 z_havew_2 : DO lhave = 1 , num_sw_levels_input +2 -1
1710 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1711 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1712 DO j = jts , MIN(jde-1,jte)
1713 DO i = its , MIN(ide-1,ite)
1714 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1715 sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1716 sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1717 ( zhave(lhave+1) - zhave(lhave) )
1718 END DO
1719 END DO
1720 EXIT z_havew_2
1721 END IF
1722 END DO z_havew_2
1723 END DO z_wantw_2
1724 ENDIF
1726 END IF
1729 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
1730 ! used, but they will make a more continuous plot.
1732 IF ( flag_sst .EQ. 1 ) THEN
1733 DO j = jts , MIN(jde-1,jte)
1734 DO i = its , MIN(ide-1,ite)
1735 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1736 IF ( landmask(i,j) .LT. 0.5 ) THEN
1737 DO l = 1 , num_soil_layers
1738 tslb(i,l,j)= sst(i,j)
1739 smois(i,l,j)= 1.0
1740 sh2o (i,l,j)= 1.0
1741 END DO
1742 END IF
1743 END DO
1744 END DO
1745 ELSE
1746 DO j = jts , MIN(jde-1,jte)
1747 DO i = its , MIN(ide-1,ite)
1748 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1749 IF ( landmask(i,j) .LT. 0.5 ) THEN
1750 DO l = 1 , num_soil_layers
1751 tslb(i,l,j)= tsk(i,j)
1752 smois(i,l,j)= 1.0
1753 sh2o (i,l,j)= 1.0
1754 END DO
1755 END IF
1756 END DO
1757 END DO
1758 END IF
1760 DEALLOCATE (zhave)
1762 END IF
1764 END SUBROUTINE init_soil_2_real
1766 SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
1767 ivgtyp,isltyp,tslb,smois,tmn, &
1768 num_soil_layers, &
1769 ids,ide, jds,jde, kds,kde, &
1770 ims,ime, jms,jme, kms,kme, &
1771 its,ite, jts,jte, kts,kte )
1773 IMPLICIT NONE
1775 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
1776 ims,ime, jms,jme, kms,kme, &
1777 its,ite, jts,jte, kts,kte
1779 INTEGER, INTENT(IN) ::num_soil_layers
1781 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
1783 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: xland
1784 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: snow, canwat, xice, vegfra, tmn
1786 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1788 INTEGER :: icm,jcm,itf,jtf
1789 INTEGER :: i,j,l
1791 itf=min0(ite,ide-1)
1792 jtf=min0(jte,jde-1)
1794 icm = ide/2
1795 jcm = jde/2
1797 DO j=jts,jtf
1798 DO l=1,num_soil_layers
1799 DO i=its,itf
1800 if (xland(i,j) .lt. 1.5) then
1801 smois(i,1,j)=0.30
1802 smois(i,2,j)=0.30
1803 smois(i,3,j)=0.30
1804 smois(i,4,j)=0.30
1806 tslb(i,1,j)=290.
1807 tslb(i,2,j)=290.
1808 tslb(i,3,j)=290.
1809 tslb(i,4,j)=290.
1810 endif
1811 ENDDO
1812 ENDDO
1813 ENDDO
1815 END SUBROUTINE init_soil_2_ideal
1817 SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
1818 st_input , sm_input , landmask, sst, &
1819 zs , dzs , flag_sm_adj, &
1820 st_levels_input , sm_levels_input , &
1821 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
1822 num_st_levels_alloc , num_sm_levels_alloc , &
1823 flag_sst , flag_soil_layers , flag_soil_levels , &
1824 ids , ide , jds , jde , kds , kde , &
1825 ims , ime , jms , jme , kms , kme , &
1826 its , ite , jts , jte , kts , kte )
1828 IMPLICIT NONE
1830 INTEGER , INTENT(IN) :: num_soil_layers , flag_sm_adj, &
1831 num_st_levels_input , num_sm_levels_input , &
1832 num_st_levels_alloc , num_sm_levels_alloc , &
1833 ids , ide , jds , jde , kds , kde , &
1834 ims , ime , jms , jme , kms , kme , &
1835 its , ite , jts , jte , kts , kte
1837 INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
1839 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1840 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1842 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1843 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1844 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1846 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1847 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1848 REAL , DIMENSION(ims:ime,jms:jme) :: smtotn, smtotr, smtotn_1m, smtotr_1m
1849 REAL , DIMENSION(num_soil_layers) :: zs , dzs, factorsm
1851 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
1853 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1855 INTEGER :: i , j , l , lout , lin , lwant , lhave, k
1856 REAL :: temp
1858 CHARACTER (LEN=132) :: message
1860 ! Allocate the soil layer array used for interpolating.
1862 IF ( ( num_st_levels_input .LE. 0 ) .OR. &
1863 ( num_sm_levels_input .LE. 0 ) ) THEN
1864 write (message, FMT='(A)')&
1865 'No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
1866 CALL wrf_error_fatal ( message )
1867 ELSE
1868 IF ( flag_soil_levels == 1 ) THEN
1869 write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1870 CALL wrf_message ( message )
1871 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
1872 ELSE
1873 write(message, FMT='(A)') ' Assume non-RUC LSM input'
1874 CALL wrf_message ( message )
1875 ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
1876 END IF
1877 END IF
1879 ! Sort the levels for temperature.
1881 outert : DO lout = 1 , num_st_levels_input-1
1882 innert : DO lin = lout+1 , num_st_levels_input
1883 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1884 temp = st_levels_input(lout)
1885 st_levels_input(lout) = st_levels_input(lin)
1886 st_levels_input(lin) = NINT(temp)
1887 DO j = jts , MIN(jde-1,jte)
1888 DO i = its , MIN(ide-1,ite)
1889 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1890 temp = st_input(i,lout,j)
1891 st_input(i,lout,j) = st_input(i,lin,j)
1892 st_input(i,lin,j) = temp
1893 END DO
1894 END DO
1895 END IF
1896 END DO innert
1897 END DO outert
1899 IF ( flag_soil_layers == 1 ) THEN
1900 DO j = jts , MIN(jde-1,jte)
1901 DO i = its , MIN(ide-1,ite)
1902 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1903 st_input(i,1,j) = tsk(i,j)
1904 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1905 END DO
1906 END DO
1907 END IF
1909 ! Sort the levels for moisture.
1911 outerm: DO lout = 1 , num_sm_levels_input-1
1912 innerm : DO lin = lout+1 , num_sm_levels_input
1913 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1914 temp = sm_levels_input(lout)
1915 sm_levels_input(lout) = sm_levels_input(lin)
1916 sm_levels_input(lin) = NINT(temp)
1917 DO j = jts , MIN(jde-1,jte)
1918 DO i = its , MIN(ide-1,ite)
1919 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1920 temp = sm_input(i,lout,j)
1921 sm_input(i,lout,j) = sm_input(i,lin,j)
1922 sm_input(i,lin,j) = temp
1923 END DO
1924 END DO
1925 END IF
1926 END DO innerm
1927 END DO outerm
1929 if( flag_soil_levels .ne. 1) then
1930 write(message, FMT='(A)') ' from Noah to RUC - compute Noah bucket'
1931 CALL wrf_message ( message )
1932 ! Compute Noah sil moisture bucket
1933 DO j = jts , MIN(jde-1,jte)
1934 DO i = its , MIN(ide-1,ite)
1935 smtotn(i,j)=sm_input(i,2,j)*0.1 + sm_input(i,3,j)*0.2 + sm_input(i,4,j)*0.7 + sm_input(i,5,j)*1.
1936 smtotn_1m(i,j)=sm_input(i,2,j)*0.1 + sm_input(i,3,j)*0.2 + sm_input(i,4,j)*0.7
1937 END DO
1938 END DO
1939 endif
1941 IF ( flag_soil_layers == 1 ) THEN
1942 DO j = jts , MIN(jde-1,jte)
1943 DO i = its , MIN(ide-1,ite)
1944 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1945 sm_input(i,1,j) = (sm_input(i,2,j)-sm_input(i,3,j))/ &
1946 (st_levels_input(2)-st_levels_input(1))*st_levels_input(1)+ &
1947 sm_input(i,2,j)
1948 ! sm_input(i,1,j) = sm_input(i,2,j)
1949 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1950 END DO
1951 END DO
1952 END IF
1954 ! Here are the levels that we have from the input for temperature.
1956 IF ( flag_soil_levels == 1 ) THEN
1957 DO l = 1 , num_st_levels_input
1958 zhave(l) = st_levels_input(l) / 100.
1959 END DO
1961 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1964 z_wantt : DO lwant = 1 , num_soil_layers
1965 z_havet : DO lhave = 1 , num_st_levels_input -1
1966 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1967 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1968 DO j = jts , MIN(jde-1,jte)
1969 DO i = its , MIN(ide-1,ite)
1970 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1971 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1972 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1973 ( zhave(lhave+1) - zhave(lhave) )
1974 END DO
1975 END DO
1976 EXIT z_havet
1977 END IF
1978 END DO z_havet
1979 END DO z_wantt
1981 ELSE
1983 zhave(1) = 0.
1984 DO l = 1 , num_st_levels_input
1985 zhave(l+1) = st_levels_input(l) / 100.
1986 END DO
1987 zhave(num_st_levels_input+2) = 300. / 100.
1989 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1991 z_wantt_2 : DO lwant = 1 , num_soil_layers
1992 z_havet_2 : DO lhave = 1 , num_st_levels_input +2
1993 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1994 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1995 DO j = jts , MIN(jde-1,jte)
1996 DO i = its , MIN(ide-1,ite)
1997 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1998 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1999 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2000 ( zhave(lhave+1) - zhave(lhave) )
2001 END DO
2002 END DO
2003 EXIT z_havet_2
2004 END IF
2005 END DO z_havet_2
2006 END DO z_wantt_2
2008 END IF
2010 ! Here are the levels that we have from the input for moisture.
2012 IF ( flag_soil_levels .EQ. 1 ) THEN
2013 DO l = 1 , num_sm_levels_input
2014 zhave(l) = sm_levels_input(l) / 100.
2015 END DO
2017 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2019 z_wantm : DO lwant = 1 , num_soil_layers
2020 z_havem : DO lhave = 1 , num_sm_levels_input -1
2021 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2022 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2023 DO j = jts , MIN(jde-1,jte)
2024 DO i = its , MIN(ide-1,ite)
2025 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2026 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2027 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2028 ( zhave(lhave+1) - zhave(lhave) )
2029 END DO
2030 END DO
2031 EXIT z_havem
2032 END IF
2033 END DO z_havem
2034 END DO z_wantm
2036 ELSE
2038 zhave(1) = 0.
2039 DO l = 1 , num_sm_levels_input
2040 zhave(l+1) = sm_levels_input(l) / 100.
2041 END DO
2042 zhave(num_sm_levels_input+2) = 300. / 100.
2044 z_wantm_2 : DO lwant = 1 , num_soil_layers
2045 z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
2046 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2047 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2048 DO j = jts , MIN(jde-1,jte)
2049 DO i = its , MIN(ide-1,ite)
2050 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2051 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2052 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2053 ( zhave(lhave+1) - zhave(lhave) )
2054 END DO
2055 END DO
2056 EXIT z_havem_2
2057 END IF
2058 END DO z_havem_2
2059 END DO z_wantm_2
2061 END IF
2063 IF(flag_sm_adj == 1) then
2064 write(message, FMT='(A)') ' Soil moisture adjustment from Noah to RUC is ON'
2065 CALL wrf_message ( message )
2066 !-- perform soil moisture adjustmen when flag_sm_adj == 1
2067 if( flag_soil_levels .ne. 1) then
2068 write(message, FMT='(A)') ' from Noah to RUC - compute RUC bucket'
2069 CALL wrf_message ( message )
2070 DO j = jts , MIN(jde-1,jte)
2071 DO i = its , MIN(ide-1,ite)
2072 smtotr(i,j)=0.
2073 smtotr_1m(i,j)=0.
2074 do k=1,num_soil_layers-1
2075 smtotr(i,j)=smtotr(i,j) + smois(i,k,j) *dzs(k)
2076 enddo
2077 do k=1,num_soil_layers-2
2078 smtotr_1m(i,j)=smtotr_1m(i,j) + smois(i,k,j) *dzs(k)
2079 enddo
2081 IF ( landmask(i,j) > 0.5) then
2082 ! land
2083 ! initialize factor
2084 do k=1,num_soil_layers
2085 factorsm(k)=1.
2086 enddo
2087 ! Correct RUC soil moisture to match Noah bucket
2088 ! print *,'1-m Buckets: RUC and Noah',i,j,smtotr_1m(i,j),0.77*smtotr_1m(i,j),smtotn_1m(i,j)
2089 ! print *,'Buckets: RUC and Noah',i,j,smtotr(i,j),smtotn(i,j)
2090 ! print *,'before smois=',i,j,smois(i,:,j)
2091 do k=1,num_soil_layers-1
2092 smois(i,k,j) = max(0.02,smois(i,k,j)*smtotn(i,j)/(0.9*smtotr(i,j)))
2093 enddo
2094 ! print *,'after smois=',i,j,smois(i,:,j)
2095 ! Compute RUC bucket after correction
2096 smtotr(i,j) = 0.
2097 do k=1,num_soil_layers-1
2098 smtotr(i,j)=smtotr(i,j) + smois(i,k,j) *dzs(k)
2099 enddo
2100 ! print *,'Buckets after correction: RUC and Noah at',i,j,smtotr(i,j),smtotn(i,j)
2101 if( smois(i,2,j) > smois(i,1,j) .and. smois(i,3,j) > smois(i,2,j)) then
2102 ! typical for daytime, no recent precip
2103 factorsm(1) = 0.75
2104 factorsm(2) = 0.8
2105 factorsm(3) = 0.85
2106 factorsm(4) = 0.9
2107 factorsm(5) = 0.95
2108 endif
2109 do k=1,num_soil_layers-1
2110 smois(i,k,j) = factorsm(k) * smois(i,k,j)
2111 enddo
2113 ! print *,'1 - after smois=',i,j,smois(i,:,j)
2114 smtotr(i,j) = 0.
2115 do k=1,num_soil_layers-1
2116 smtotr(i,j)=smtotr(i,j) + smois(i,k,j) *dzs(k)
2117 enddo
2118 ! print *,'1 - Buckets after correction: RUC and Noah at',i,j,smtotr(i,j),smtotn(i,j)
2119 ENDIF ! land
2120 END DO
2121 END DO
2123 endif ! flag_soil_levels
2124 ENDIF ! flag_sm_adj
2126 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2127 ! used, but they will make a more continuous plot.
2129 IF ( flag_sst .EQ. 1 ) THEN
2130 DO j = jts , MIN(jde-1,jte)
2131 DO i = its , MIN(ide-1,ite)
2132 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2133 IF ( landmask(i,j) .LT. 0.5 ) THEN
2134 DO l = 1 , num_soil_layers
2135 tslb(i,l,j) = sst(i,j)
2136 tsk(i,j) = sst(i,j)
2137 smois(i,l,j)= 1.0
2138 END DO
2139 END IF
2140 END DO
2141 END DO
2142 ELSE
2143 DO j = jts , MIN(jde-1,jte)
2144 DO i = its , MIN(ide-1,ite)
2145 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2146 IF ( landmask(i,j) .LT. 0.5 ) THEN
2147 DO l = 1 , num_soil_layers
2148 tslb(i,l,j)= tsk(i,j)
2149 smois(i,l,j)= 1.0
2150 END DO
2151 END IF
2152 END DO
2153 END DO
2154 END IF
2156 DEALLOCATE (zhave)
2158 END SUBROUTINE init_soil_3_real
2160 SUBROUTINE init_soil_4_real ( tsk , tmn , smois , sh2o , tslb , &
2161 st_input , sm_input , sw_input , landmask , sst , &
2162 zs , dzs , &
2163 st_levels_input , sm_levels_input , sw_levels_input , &
2164 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2165 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2166 flag_sst , flag_soil_layers , flag_soil_levels , &
2167 ids , ide , jds , jde , kds , kde , &
2168 ims , ime , jms , jme , kms , kme , &
2169 its , ite , jts , jte , kts , kte )
2171 IMPLICIT NONE
2173 INTEGER , INTENT(IN) :: num_soil_layers , &
2174 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2175 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2176 ids , ide , jds , jde , kds , kde , &
2177 ims , ime , jms , jme , kms , kme , &
2178 its , ite , jts , jte , kts , kte
2180 INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
2182 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2183 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2184 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2186 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2187 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2188 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2189 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2191 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2192 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2193 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2195 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2197 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2199 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2200 REAL :: temp
2201 LOGICAL :: found_levels
2203 CHARACTER (LEN=132) :: message
2205 ! Are there any soil temp and moisture levels - ya know, they are
2206 ! mandatory.
2208 num = num_st_levels_input * num_sm_levels_input
2210 IF ( num .GE. 1 ) THEN
2212 ! Ordered levels that we have data for.
2214 !tgs add option to initialize from RUCLSM
2215 IF ( flag_soil_levels == 1 ) THEN
2216 write(message, FMT='(A)') ' Assume CLM input'
2217 CALL wrf_message ( message )
2218 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
2219 ELSE
2220 write(message, FMT='(A)') ' Assume non-CLM input'
2221 CALL wrf_message ( message )
2222 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2223 ! ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
2224 END IF
2227 ! Sort the levels for temperature.
2229 outert : DO lout = 1 , num_st_levels_input-1
2230 innert : DO lin = lout+1 , num_st_levels_input
2231 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2232 temp = st_levels_input(lout)
2233 st_levels_input(lout) = st_levels_input(lin)
2234 st_levels_input(lin) = NINT(temp)
2235 DO j = jts , MIN(jde-1,jte)
2236 DO i = its , MIN(ide-1,ite)
2237 temp = st_input(i,lout+1,j)
2238 st_input(i,lout+1,j) = st_input(i,lin+1,j)
2239 st_input(i,lin+1,j) = temp
2240 END DO
2241 END DO
2242 END IF
2243 END DO innert
2244 END DO outert
2245 !tgs add IF
2246 IF ( flag_soil_layers == 1 ) THEN
2247 DO j = jts , MIN(jde-1,jte)
2248 DO i = its , MIN(ide-1,ite)
2249 st_input(i,1,j) = tsk(i,j)
2250 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2251 END DO
2252 END DO
2253 ENDIF
2254 ! Sort the levels for moisture.
2256 outerm: DO lout = 1 , num_sm_levels_input-1
2257 innerm : DO lin = lout+1 , num_sm_levels_input
2258 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2259 temp = sm_levels_input(lout)
2260 sm_levels_input(lout) = sm_levels_input(lin)
2261 sm_levels_input(lin) = NINT(temp)
2262 DO j = jts , MIN(jde-1,jte)
2263 DO i = its , MIN(ide-1,ite)
2264 temp = sm_input(i,lout+1,j)
2265 sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2266 sm_input(i,lin+1,j) = temp
2267 END DO
2268 END DO
2269 END IF
2270 END DO innerm
2271 END DO outerm
2272 !tgs add IF
2273 IF ( flag_soil_layers == 1 ) THEN
2274 DO j = jts , MIN(jde-1,jte)
2275 DO i = its , MIN(ide-1,ite)
2276 sm_input(i,1,j) = sm_input(i,2,j)
2277 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2278 END DO
2279 END DO
2280 ENDIF
2282 ! Sort the levels for liquid moisture.
2284 outerw: DO lout = 1 , num_sw_levels_input-1
2285 innerw : DO lin = lout+1 , num_sw_levels_input
2286 IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2287 temp = sw_levels_input(lout)
2288 sw_levels_input(lout) = sw_levels_input(lin)
2289 sw_levels_input(lin) = NINT(temp)
2290 DO j = jts , MIN(jde-1,jte)
2291 DO i = its , MIN(ide-1,ite)
2292 temp = sw_input(i,lout+1,j)
2293 sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2294 sw_input(i,lin+1,j) = temp
2295 END DO
2296 END DO
2297 END IF
2298 END DO innerw
2299 END DO outerw
2300 IF ( num_sw_levels_input .GT. 1 ) THEN
2301 DO j = jts , MIN(jde-1,jte)
2302 DO i = its , MIN(ide-1,ite)
2303 sw_input(i,1,j) = sw_input(i,2,j)
2304 sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2305 END DO
2306 END DO
2307 END IF
2309 found_levels = .TRUE.
2311 ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2313 found_levels = .FALSE.
2315 ELSE
2316 CALL wrf_error_fatal ( &
2317 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
2318 END IF
2320 ! Is it OK to continue?
2322 IF ( found_levels ) THEN
2324 !tgs: Here are the levels that we have from the input for temperature.
2326 IF ( flag_soil_levels == 1 ) THEN
2327 DO l = 1 , num_st_levels_input
2328 zhave(l) = st_levels_input(l) / 100.
2329 END DO
2331 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2333 z_wantt : DO lwant = 1 , num_soil_layers
2334 z_havet : DO lhave = 1 , num_st_levels_input -1
2335 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2336 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2337 DO j = jts , MIN(jde-1,jte)
2338 DO i = its , MIN(ide-1,ite)
2339 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2340 st_input(i,lhave+1,j) * ( zs (lwant) - zhave(lhave) ) ) / &
2341 ( zhave(lhave+1) - zhave(lhave) )
2342 END DO
2343 END DO
2344 EXIT z_havet
2345 END IF
2346 END DO z_havet
2347 END DO z_wantt
2349 ELSE
2351 ! Here are the levels that we have from the input for temperature. The input levels plus
2352 ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2354 zhave(1) = 0.
2355 DO l = 1 , num_st_levels_input
2356 zhave(l+1) = st_levels_input(l) / 100.
2357 END DO
2358 zhave(num_st_levels_input+2) = 300. / 100.
2360 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2362 z_wantt_2: DO lwant = 1 , num_soil_layers
2363 z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
2364 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2365 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2366 DO j = jts , MIN(jde-1,jte)
2367 DO i = its , MIN(ide-1,ite)
2368 tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2369 st_input(i,lhave+1,j) * ( zs (lwant) - zhave(lhave) ) ) / &
2370 ( zhave(lhave+1) - zhave(lhave) )
2371 END DO
2372 END DO
2373 EXIT z_havet_2
2374 END IF
2375 END DO z_havet_2
2376 END DO z_wantt_2
2378 END IF
2380 !tgs:
2381 IF ( flag_soil_levels == 1 ) THEN
2382 DO l = 1 , num_sm_levels_input
2383 zhave(l) = sm_levels_input(l) / 100.
2384 END DO
2386 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2388 z_wantm : DO lwant = 1 , num_soil_layers
2389 z_havem : DO lhave = 1 , num_sm_levels_input -1
2390 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2391 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2392 DO j = jts , MIN(jde-1,jte)
2393 DO i = its , MIN(ide-1,ite)
2394 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2395 sm_input(i,lhave+1,j) * ( zs (lwant) - zhave(lhave) ) ) / &
2396 ( zhave(lhave+1) - zhave(lhave) )
2397 if(smois(i,lwant,j)<=0.0) smois(i,lwant,j) = 0.005
2398 END DO
2399 END DO
2400 EXIT z_havem
2401 END IF
2402 END DO z_havem
2403 END DO z_wantm
2405 ELSE
2406 ! Here are the levels that we have from the input for moisture. The input levels plus
2407 ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2408 ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2409 ! as the most deep layer's value.
2411 zhave(1) = 0.
2412 DO l = 1 , num_sm_levels_input
2413 zhave(l+1) = sm_levels_input(l) / 100.
2414 END DO
2415 zhave(num_sm_levels_input+2) = 300. / 100.
2417 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2419 z_wantm_2 : DO lwant = 1 , num_soil_layers
2420 z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
2421 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2422 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2423 DO j = jts , MIN(jde-1,jte)
2424 DO i = its , MIN(ide-1,ite)
2425 smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2426 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2427 ( zhave(lhave+1) - zhave(lhave) )
2428 if(smois(i,lwant,j)<=0.0) smois(i,lwant,j) = 0.005
2429 END DO
2430 END DO
2431 EXIT z_havem_2
2432 END IF
2433 END DO z_havem_2
2434 END DO z_wantm_2
2435 ENDIF
2437 ! Any liquid soil moisture to worry about?
2439 IF ( num_sw_levels_input .GT. 1 ) THEN
2441 zhave(1) = 0.
2442 DO l = 1 , num_sw_levels_input
2443 zhave(l+1) = sw_levels_input(l) / 100.
2444 END DO
2445 zhave(num_sw_levels_input+2) = 300. / 100.
2447 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2449 z_wantw : DO lwant = 1 , num_soil_layers
2450 z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
2451 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2452 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2453 DO j = jts , MIN(jde-1,jte)
2454 DO i = its , MIN(ide-1,ite)
2455 sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2456 sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2457 ( zhave(lhave+1) - zhave(lhave) )
2458 END DO
2459 END DO
2460 EXIT z_havew
2461 END IF
2462 END DO z_havew
2463 END DO z_wantw
2465 END IF
2468 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2469 ! used, but they will make a more continuous plot.
2471 IF ( flag_sst .EQ. 1 ) THEN
2472 DO j = jts , MIN(jde-1,jte)
2473 DO i = its , MIN(ide-1,ite)
2474 IF ( landmask(i,j) .LT. 0.5 ) THEN
2475 DO l = 1 , num_soil_layers
2476 tslb(i,l,j)= sst(i,j)
2477 smois(i,l,j)= 1.0
2478 sh2o (i,l,j)= 1.0
2479 END DO
2480 END IF
2481 END DO
2482 END DO
2483 ELSE
2484 DO j = jts , MIN(jde-1,jte)
2485 DO i = its , MIN(ide-1,ite)
2486 IF ( landmask(i,j) .LT. 0.5 ) THEN
2487 DO l = 1 , num_soil_layers
2488 tslb(i,l,j)= tsk(i,j)
2489 smois(i,l,j)= 1.0
2490 sh2o (i,l,j)= 1.0
2491 END DO
2492 END IF
2493 END DO
2494 END DO
2495 END IF
2497 DEALLOCATE (zhave)
2499 END IF
2501 END SUBROUTINE init_soil_4_real
2503 SUBROUTINE init_soil_7_real ( tsk , tmn , smois , sh2o , tslb , &
2504 st_input , sm_input , sw_input , landmask , sst , &
2505 zs , dzs , &
2506 st_levels_input , sm_levels_input , sw_levels_input , &
2507 num_soil_layers , num_st_levels_input , &
2508 num_sm_levels_input , num_sw_levels_input , &
2509 num_st_levels_alloc , num_sm_levels_alloc , &
2510 num_sw_levels_alloc , &
2511 flag_sst , flag_soil_layers , flag_soil_levels , &
2512 ids , ide , jds , jde , kds , kde , &
2513 ims , ime , jms , jme , kms , kme , &
2514 its , ite , jts , jte , kts , kte )
2516 ! for soil temperature and moisture initialization for the PX LSM
2518 IMPLICIT NONE
2520 INTEGER , INTENT(IN) :: num_soil_layers , &
2521 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2522 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2523 ids , ide , jds , jde , kds , kde , &
2524 ims , ime , jms , jme , kms , kme , &
2525 its , ite , jts , jte , kts , kte
2527 INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
2529 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2530 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2531 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2533 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2534 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2535 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2536 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2538 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2539 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2540 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2542 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2544 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2546 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2547 REAL :: temp
2548 LOGICAL :: found_levels
2550 ! Are there any soil temp and moisture levels - ya know, they are mandatory.
2552 num = num_st_levels_input * num_sm_levels_input
2554 IF ( num .GE. 1 ) THEN
2556 ! Ordered levels that we have data for.
2558 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2560 ! Sort the levels for temperature.
2561 outert : DO lout = 1 , num_st_levels_input-1
2562 innert : DO lin = lout+1 , num_st_levels_input
2563 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2564 temp = st_levels_input(lout)
2565 st_levels_input(lout) = st_levels_input(lin)
2566 st_levels_input(lin) = NINT(temp)
2567 DO j = jts , MIN(jde-1,jte)
2568 DO i = its , MIN(ide-1,ite)
2569 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2570 temp = st_input(i,lout+1,j)
2571 st_input(i,lout+1,j) = st_input(i,lin+1,j)
2572 st_input(i,lin+1,j) = temp
2573 END DO
2574 END DO
2575 END IF
2576 END DO innert
2577 END DO outert
2578 DO j = jts , MIN(jde-1,jte)
2579 DO i = its , MIN(ide-1,ite)
2580 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2581 st_input(i,1,j) = tsk(i,j)
2582 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2583 END DO
2584 END DO
2586 ! Sort the levels for moisture.
2588 outerm: DO lout = 1 , num_sm_levels_input-1
2589 innerm : DO lin = lout+1 , num_sm_levels_input
2590 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2591 temp = sm_levels_input(lout)
2592 sm_levels_input(lout) = sm_levels_input(lin)
2593 sm_levels_input(lin) = NINT(temp)
2594 DO j = jts , MIN(jde-1,jte)
2595 DO i = its , MIN(ide-1,ite)
2596 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2597 temp = sm_input(i,lout+1,j)
2598 sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2599 sm_input(i,lin+1,j) = temp
2600 END DO
2601 END DO
2602 END IF
2603 END DO innerm
2604 END DO outerm
2605 DO j = jts , MIN(jde-1,jte)
2606 DO i = its , MIN(ide-1,ite)
2607 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2608 sm_input(i,1,j) = sm_input(i,2,j)
2609 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2610 END DO
2611 END DO
2613 ! Sort the levels for liquid moisture.
2615 outerw: DO lout = 1 , num_sw_levels_input-1
2616 innerw : DO lin = lout+1 , num_sw_levels_input
2617 IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2618 temp = sw_levels_input(lout)
2619 sw_levels_input(lout) = sw_levels_input(lin)
2620 sw_levels_input(lin) = NINT(temp)
2621 DO j = jts , MIN(jde-1,jte)
2622 DO i = its , MIN(ide-1,ite)
2623 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2624 temp = sw_input(i,lout+1,j)
2625 sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2626 sw_input(i,lin+1,j) = temp
2627 END DO
2628 END DO
2629 END IF
2630 END DO innerw
2631 END DO outerw
2632 IF ( num_sw_levels_input .GT. 1 ) THEN
2633 DO j = jts , MIN(jde-1,jte)
2634 DO i = its , MIN(ide-1,ite)
2635 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2636 sw_input(i,1,j) = sw_input(i,2,j)
2637 sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2638 END DO
2639 END DO
2640 END IF
2642 found_levels = .TRUE.
2644 ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2646 found_levels = .FALSE.
2648 ELSE
2649 CALL wrf_error_fatal ( &
2650 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for PX LSM.' )
2651 END IF
2653 ! Is it OK to continue?
2655 IF ( found_levels ) THEN
2657 ! Here are the levels that we have from the input for temperature. The input levels plus
2658 ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2660 zhave(1) = 0.
2661 DO l = 1 , num_st_levels_input
2662 zhave(l+1) = st_levels_input(l) / 100.
2663 END DO
2664 zhave(num_st_levels_input+2) = 300. / 100.
2666 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2668 z_wantt : DO lwant = 1 , num_soil_layers
2669 z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2670 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2671 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2672 DO j = jts , MIN(jde-1,jte)
2673 DO i = its , MIN(ide-1,ite)
2674 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2675 tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2676 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2677 ( zhave(lhave+1) - zhave(lhave) )
2678 END DO
2679 END DO
2680 EXIT z_havet
2681 END IF
2682 END DO z_havet
2683 END DO z_wantt
2685 ! Here are the levels that we have from the input for moisture. The input levels plus
2686 ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2687 ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2688 ! as the most deep layer's value.
2690 zhave(1) = 0.
2691 DO l = 1 , num_sm_levels_input
2692 zhave(l+1) = sm_levels_input(l) / 100.
2693 END DO
2694 zhave(num_sm_levels_input+2) = 300. / 100.
2696 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2698 z_wantm : DO lwant = 1 , num_soil_layers
2699 z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2700 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2701 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2702 DO j = jts , MIN(jde-1,jte)
2703 DO i = its , MIN(ide-1,ite)
2704 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2705 smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2706 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2707 ( zhave(lhave+1) - zhave(lhave) )
2708 END DO
2709 END DO
2710 EXIT z_havem
2711 END IF
2712 END DO z_havem
2713 END DO z_wantm
2715 ! Any liquid soil moisture to worry about?
2717 IF ( num_sw_levels_input .GT. 1 ) THEN
2719 zhave(1) = 0.
2720 DO l = 1 , num_sw_levels_input
2721 zhave(l+1) = sw_levels_input(l) / 100.
2722 END DO
2723 zhave(num_sw_levels_input+2) = 300. / 100.
2725 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2727 z_wantw : DO lwant = 1 , num_soil_layers
2728 z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
2729 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2730 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2731 DO j = jts , MIN(jde-1,jte)
2732 DO i = its , MIN(ide-1,ite)
2733 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2734 sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2735 sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2736 ( zhave(lhave+1) - zhave(lhave) )
2737 END DO
2738 END DO
2739 EXIT z_havew
2740 END IF
2741 END DO z_havew
2742 END DO z_wantw
2744 END IF
2747 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2748 ! used, but they will make a more continuous plot.
2750 DO j = jts , MIN(jde-1,jte)
2751 DO i = its , MIN(ide-1,ite)
2752 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2753 tslb(i,1,j)= tsk(i,j)
2754 tslb(i,2,j)= tmn(i,j)
2755 END DO
2756 END DO
2758 IF ( flag_sst .EQ. 1 ) THEN
2759 DO j = jts , MIN(jde-1,jte)
2760 DO i = its , MIN(ide-1,ite)
2761 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2762 IF ( landmask(i,j) .LT. 0.5 ) THEN
2763 DO l = 1 , num_soil_layers
2764 tslb(i,l,j)= sst(i,j)
2765 smois(i,l,j)= 1.0
2766 sh2o (i,l,j)= 1.0
2767 END DO
2768 END IF
2769 END DO
2770 END DO
2771 ELSE
2772 DO j = jts , MIN(jde-1,jte)
2773 DO i = its , MIN(ide-1,ite)
2774 IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2775 IF ( landmask(i,j) .LT. 0.5 ) THEN
2776 DO l = 1 , num_soil_layers
2777 tslb(i,l,j)= tsk(i,j)
2778 smois(i,l,j)= 1.0
2779 sh2o (i,l,j)= 1.0
2780 END DO
2781 END IF
2782 END DO
2783 END DO
2784 END IF
2786 DEALLOCATE (zhave)
2788 END IF
2790 END SUBROUTINE init_soil_7_real
2791 !------------SSIB (fds 06/2010)-----------------------------
2792 SUBROUTINE init_soil_8_real ( tsk , tmn , smois , sh2o , tslb , &
2793 st_input , sm_input , sw_input , landmask , sst , &
2794 zs , dzs , &
2795 st_levels_input , sm_levels_input , sw_levels_input , &
2796 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2797 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2798 flag_sst , flag_soil_layers , flag_soil_levels , &
2799 st000010, st010200, sm000010, sm010200, &
2800 st010040 , st040100 , st100200, &
2801 sm010040 , sm040100 , sm100200, &
2802 ids , ide , jds , jde , kds , kde , &
2803 ims , ime , jms , jme , kms , kme , &
2804 its , ite , jts , jte , kts , kte )
2806 IMPLICIT NONE
2808 INTEGER , INTENT(IN) :: num_soil_layers , &
2809 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2810 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2811 ids , ide , jds , jde , kds , kde , &
2812 ims , ime , jms , jme , kms , kme , &
2813 its , ite , jts , jte , kts , kte
2815 INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers , flag_soil_levels
2817 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2818 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2819 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2821 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2822 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2823 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2824 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2825 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: st000010, st010200, st010040 , st040100 , st100200
2826 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: sm000010, sm010200, sm010040 , sm040100 , sm100200
2828 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2829 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2830 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2832 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2834 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2836 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2837 REAL :: temp
2838 LOGICAL :: found_levels
2840 found_levels = .false.
2842 ! write(message, FMT='(A)') ' Prepare SSiB LSM soil fields'
2843 !fds (06/2010)
2844 !-- SSiB needs NCEP/NCAR reanalysis soil temperature and moisture at 0-10cm and 10-200cm --
2845 IF ( flag_soil_layers == 1 ) THEN
2846 DO j = jts , MIN(jde-1,jte)
2847 DO i = its , MIN(ide-1,ite)
2848 tslb(i,1,j) = st000010(i,j)
2849 tslb(i,2,j) = st010040(i,j)
2850 tslb(i,3,j) = st040100(i,j)
2851 smois(i,1,j) = sm000010(i,j)
2852 smois(i,2,j) = sm010040(i,j)
2853 smois(i,3,j) = sm040100(i,j)
2854 ! temporary fix (JD) for missing 4-layer input - use 2 instead
2855 if(tslb(i,2,j) .lt. 10.)tslb(i,2,j) = st000010(i,j)
2856 if(smois(i,2,j) .lt. 0.01)smois(i,2,j) = sm000010(i,j)
2857 if(tslb(i,3,j) .lt. 10.)tslb(i,3,j) = st010200(i,j)
2858 if(smois(i,3,j) .lt. 0.01)smois(i,3,j) = sm010200(i,j)
2859 END DO
2860 END DO
2861 ELSE
2862 CALL wrf_debug ( 0 , 'SSiB LSM needs 0-10 cm soil t')
2863 CALL wrf_debug ( 0 , 'SSiB LSM needs 0-10 cm soil m')
2864 CALL wrf_debug ( 0 , 'SSiB LSM needs 10-200 cm soil t')
2865 CALL wrf_debug ( 0 , 'SSiB LSM needs 10-200 cm soil m')
2866 ENDIF
2868 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2869 ! used, but they will make a more continuous plot.
2871 IF ( flag_sst .EQ. 1 ) THEN
2872 DO j = jts , MIN(jde-1,jte)
2873 DO i = its , MIN(ide-1,ite)
2874 IF ( landmask(i,j) .LT. 0.5 ) THEN
2875 DO l = 1 , num_soil_layers
2876 tslb(i,l,j)= sst(i,j)
2877 smois(i,l,j)= 1.0
2878 sh2o (i,l,j)= 1.0
2879 END DO
2880 END IF
2881 END DO
2882 END DO
2883 ELSE
2884 DO j = jts , MIN(jde-1,jte)
2885 DO i = its , MIN(ide-1,ite)
2886 IF ( landmask(i,j) .LT. 0.5 ) THEN
2887 DO l = 1 , num_soil_layers
2888 tslb(i,l,j)= tsk(i,j)
2889 smois(i,l,j)= 1.0
2890 sh2o (i,l,j)= 1.0
2891 END DO
2892 END IF
2893 END DO
2894 END DO
2895 END IF
2897 END SUBROUTINE init_soil_8_real
2898 !--------------------------------------------------------
2900 SUBROUTINE aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu )
2901 IMPLICIT NONE
2902 INTEGER , INTENT(IN) :: iswater
2903 INTEGER , INTENT(IN) :: num_veg_cat
2904 CHARACTER (LEN=4) , INTENT(IN) :: mminlu
2905 REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
2907 ! Local variables
2909 INTEGER , PARAMETER :: num_special_bins = 3 ! grass, shrubs, trees; add a new line in the "cib" array if updating this
2910 INTEGER , PARAMETER :: max_cats_per_bin = 8 ! larger than max num of cats inside each of the special bins
2911 ! So, no more than 8 total tree categories that we need to consider.
2912 INTEGER , PARAMETER :: fl = -1 ! Flag value so that we know this is not a valid category to consider.
2913 INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_usgs = &
2914 (/ 2 , 3 , 4 , 5 , 6 , 7 , 10 , fl , & ! grass
2915 8 , 9 , fl , fl , fl , fl , fl , fl , & ! shrubs
2916 11 , 12 , 13 , 14 , 15 , fl , fl , fl /) ! trees
2917 INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_modis = &
2918 (/ 9 , 10 , 12 , 14 , fl , fl , fl , fl , & ! grass
2919 6 , 7 , 8 , fl , fl , fl , fl , fl , & ! shrubs
2920 1 , 2 , 3 , 4 , 5 , fl , fl , fl /) ! trees
2922 IF ( mminlu(1:4) .EQ. 'USGS' ) THEN
2923 CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2924 num_special_bins , max_cats_per_bin , fl , cib_usgs )
2925 ELSE IF ( mminlu(1:4) .EQ. 'MODI' ) THEN
2926 CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2927 num_special_bins , max_cats_per_bin , fl , cib_modis )
2928 END IF
2930 END SUBROUTINE aggregate_categories_part1
2932 SUBROUTINE aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2933 num_special_bins , max_cats_per_bin , fl , cib )
2934 IMPLICIT NONE
2936 INTEGER , INTENT(IN) :: iswater
2937 INTEGER , INTENT(IN) :: num_veg_cat
2938 REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
2939 INTEGER , INTENT(IN) :: num_special_bins , max_cats_per_bin , fl
2940 INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) , INTENT(IN) :: cib
2942 ! Local variables
2944 REAL , DIMENSION(1:num_veg_cat) :: landuse_frac_work ! copy of the input array, allows us to be wreckless
2945 INTEGER , DIMENSION ( max_cats_per_bin , num_special_bins ) :: cats_in_bin
2946 INTEGER :: ib , ic ! indexes for the bin and the category loops
2947 REAL , DIMENSION ( num_special_bins ) :: bin_max_val , bin_sum ! max category value in this bin, sum of all cats in this bin
2948 INTEGER , DIMENSION ( num_special_bins ) :: bin_max_idx ! index of this maximum category in this bin
2949 INTEGER :: bin_work , bin_orig ! the bin from whence the maximum hails, respectively for the aggregated and the original data
2950 INTEGER :: max_cat_orig , max_cat_work
2951 REAL :: max_val_orig , max_val_work
2953 ! Find the max in the original. If it is >= 50%, no need to even be in here.
2955 DO ic = 1 , num_veg_cat
2956 IF ( landuse_frac(ic) .GE. 0.5 ) THEN
2957 RETURN
2958 END IF
2959 END DO
2961 ! Put the categories in the bin into a 2d array.
2963 cats_in_bin = RESHAPE ( cib , (/ max_cats_per_bin , num_special_bins /) )
2965 ! Make a copy of the incoming array so that we can eventually diff our working copy with
2966 ! the original.
2968 landuse_frac_work = landuse_frac
2970 ! Loop over each of the special bins that we know about. Find the max values and the locations of such.
2972 DO ib = 1 , num_special_bins
2974 ! For this bin, we know about specific categories. We get the sum of all of those
2975 ! categories that we know about for this bin, and we keep track of which category
2976 ! is dominant.
2978 bin_sum (ib) = 0
2979 bin_max_val(ib) = fl
2981 cat_loop_accumulate : DO ic = 1 , max_cats_per_bin
2983 ! Have we run out of valid categories in this bin? For example, we would not necessarily have
2984 ! the same number of tree categories as there are grass categories.
2986 IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
2987 EXIT cat_loop_accumulate
2988 END IF
2990 bin_sum(ib) = bin_sum(ib) + landuse_frac(cats_in_bin(ic,ib))
2992 IF ( landuse_frac(cats_in_bin(ic,ib)) .GT. bin_max_val(ib) ) THEN
2993 bin_max_val(ib) = landuse_frac(cats_in_bin(ic,ib))
2994 bin_max_idx(ib) = cats_in_bin(ic,ib)
2995 END IF
2997 END DO cat_loop_accumulate
2999 cat_loop_assign : DO ic = 1 , max_cats_per_bin
3001 ! Plow through each cat in the bin. If we find the dominant one, he gets the total sum. If we land on
3002 ! the other guys in the bin, they get set back to zero to maintain the original aggregate influence.
3004 IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
3005 EXIT cat_loop_assign
3006 ELSE IF ( cats_in_bin(ic,ib) .EQ. bin_max_idx(ib) ) THEN
3007 landuse_frac_work(cats_in_bin(ic,ib)) = bin_sum(ib)
3008 ELSE
3009 landuse_frac_work(cats_in_bin(ic,ib)) = 0
3010 END IF
3012 END DO cat_loop_assign
3014 END DO
3016 ! Now we loop through the categorical data, and get the max+location for the original input data and the
3017 ! modified work data. Water is not allowed to be a "max" category unless it is greater than 50%. We hit
3018 ! that test up at the top already, so we can toss out water willy nilly here.
3020 max_cat_orig = fl
3021 max_val_orig = 0
3022 max_cat_work = fl
3023 max_val_work = 0
3025 DO ic = 1 , num_veg_cat
3026 IF ( ic .EQ. iswater ) THEN
3027 CYCLE
3028 END IF
3029 IF ( landuse_frac(ic) .GT. max_val_orig ) THEN
3030 max_val_orig = landuse_frac(ic)
3031 max_cat_orig = ic
3032 END IF
3033 IF ( landuse_frac_work(ic) .GT. max_val_work ) THEN
3034 max_val_work = landuse_frac_work(ic)
3035 max_cat_work = ic
3036 END IF
3037 END DO
3039 ! Find the bin for the maximimum value of the original data.
3041 bin_orig = -1
3042 bin_loop_orig : DO ib = 1 , num_special_bins
3043 cat_loop_orig : DO ic = 1 , max_cats_per_bin
3044 IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
3045 EXIT cat_loop_orig
3046 ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_orig ) THEN
3047 bin_orig = ib
3048 EXIT bin_loop_orig
3049 END IF
3050 END DO cat_loop_orig
3051 END DO bin_loop_orig
3053 ! Find the bin for the maximimum value of the aggregated data.
3055 bin_work = -1
3056 bin_loop_work : DO ib = 1 , num_special_bins
3057 cat_loop_work : DO ic = 1 , max_cats_per_bin
3058 IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
3059 EXIT cat_loop_work
3060 ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_work ) THEN
3061 bin_work = ib
3062 EXIT bin_loop_work
3063 END IF
3064 END DO cat_loop_work
3065 END DO bin_loop_work
3067 ! So the big question is, did the aggregation change the bin? If the aggregation does not change the resulting
3068 ! bin, then we leave everything alone. However, if there would be a change in the eventual bin chosen by
3069 ! the simple dominant-category method, then we become pro-active interventionists and rectify this heinous
3070 ! injustice foisted upon the lowly flora.
3072 IF ( bin_work .EQ. bin_orig ) THEN
3073 ! No op, we do nothing.
3074 ELSE
3075 DO ic = 1 , max_cats_per_bin
3076 landuse_frac(cats_in_bin(ic,bin_work)) = landuse_frac_work(cats_in_bin(ic,bin_work))
3077 END DO
3078 END IF
3080 END SUBROUTINE aggregate_categories_part2
3082 END MODULE module_soil_pre
3084 FUNCTION skip_middle_points_t ( ids , ide , jds , jde , &
3085 i_in , j_in , width , &
3086 subtleties_exist ) &
3087 RESULT ( skip_it )
3089 IMPLICIT NONE
3091 INTEGER , INTENT(IN) :: ids , ide , jds , jde
3092 INTEGER , INTENT(IN) :: i_in , j_in , width
3093 LOGICAL , INTENT(IN) :: subtleties_exist
3095 LOGICAL :: skip_it
3097 INTEGER , PARAMETER :: slop = 0
3099 IF ( ( subtleties_exist ) .OR. ( ide .EQ. 3 ) .OR. ( jde .EQ. 3 ) ) THEN
3100 skip_it = .FALSE.
3101 ELSE
3102 IF ( ( i_in .GE. ids+width+slop ) .AND. ( i_in .LE. ide-1-width-slop ) .AND. &
3103 ( j_in .GE. jds+width+slop ) .AND. ( j_in .LE. jde-1-width-slop ) ) THEN
3104 skip_it = .TRUE.
3105 ELSE
3106 skip_it = .FALSE.
3107 END IF
3108 END IF
3110 END FUNCTION skip_middle_points_t