module_soil_pre.F
References to this file elsewhere.
1 #if ( ! NMM_CORE == 1 )
2
3 MODULE module_soil_pre
4
5 USE module_date_time
6 USE module_state_description
7
8 CONTAINS
9
10 SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
11 xland , landusef , isltyp , soilcat , soilctop , &
12 soilcbot , tmn , &
13 seaice_threshold , &
14 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
15 iswater , isice , &
16 scheme , &
17 ids , ide , jds , jde , kds , kde , &
18 ims , ime , jms , jme , kms , kme , &
19 its , ite , jts , jte , kts , kte )
20
21 IMPLICIT NONE
22
23 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
24 ims , ime , jms , jme , kms , kme , &
25 its , ite , jts , jte , kts , kte , &
26 iswater , isice
27 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
28
29 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
30 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
31 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
32 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
33 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
34 vegcat, xland , soilcat , tmn
35 REAL , INTENT(IN) :: seaice_threshold
36
37 INTEGER :: i , j , num_seaice_changes , loop
38 CHARACTER (LEN=132) :: message
39
40 num_seaice_changes = 0
41 fix_seaice : SELECT CASE ( scheme )
42
43 CASE ( SLABSCHEME )
44 DO j = jts , MIN(jde-1,jte)
45 DO i = its , MIN(ide-1,ite)
46 IF ( xice(i,j) .GT. 200.0 ) THEN
47 xice(i,j) = 0.
48 num_seaice_changes = num_seaice_changes + 1
49 END IF
50 END DO
51 END DO
52 IF ( num_seaice_changes .GT. 0 ) THEN
53 WRITE ( message , FMT='(A,I6)' ) &
54 'Total pre number of sea ice locations removed (due to FLAG values) = ', &
55 num_seaice_changes
56 CALL wrf_debug ( 0 , message )
57 END IF
58 num_seaice_changes = 0
59 DO j = jts , MIN(jde-1,jte)
60 DO i = its , MIN(ide-1,ite)
61 IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
62 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
63 xice(i,j) = 1.
64 num_seaice_changes = num_seaice_changes + 1
65 if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
66 vegcat(i,j)=isice
67 ivgtyp(i,j)=isice
68 lu_index(i,j)=isice
69 landmask(i,j)=1.
70 xland(i,j)=1.
71 DO loop=1,num_veg_cat
72 landusef(i,loop,j)=0.
73 END DO
74 landusef(i,ivgtyp(i,j),j)=1.
75
76 isltyp(i,j) = 16
77 soilcat(i,j)=isltyp(i,j)
78 DO loop=1,num_soil_top_cat
79 soilctop(i,loop,j)=0
80 END DO
81 DO loop=1,num_soil_bot_cat
82 soilcbot(i,loop,j)=0
83 END DO
84 soilctop(i,isltyp(i,j),j)=1.
85 soilcbot(i,isltyp(i,j),j)=1.
86 END IF
87 END DO
88 END DO
89 IF ( num_seaice_changes .GT. 0 ) THEN
90 WRITE ( message , FMT='(A,I6)' ) &
91 'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
92 CALL wrf_debug ( 0 , message )
93 END IF
94
95 CASE ( LSMSCHEME , RUCLSMSCHEME )
96 num_seaice_changes = 0
97 DO j = jts , MIN(jde-1,jte)
98 DO i = its , MIN(ide-1,ite)
99 IF ( landmask(i,j) .GT. 0.5 ) THEN
100 if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
101 xice(i,j) = 0.
102 END IF
103 END DO
104 END DO
105 IF ( num_seaice_changes .GT. 0 ) THEN
106 WRITE ( message , FMT='(A,I6)' ) &
107 'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
108 CALL wrf_debug ( 0 , message )
109 END IF
110
111 END SELECT fix_seaice
112
113 END SUBROUTINE adjust_for_seaice_pre
114
115 SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
116 xland , landusef , isltyp , soilcat , soilctop , &
117 soilcbot , tmn , vegfra , &
118 tslb , smois , sh2o , &
119 seaice_threshold , &
120 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
121 num_soil_layers , &
122 iswater , isice , &
123 scheme , &
124 ids , ide , jds , jde , kds , kde , &
125 ims , ime , jms , jme , kms , kme , &
126 its , ite , jts , jte , kts , kte )
127
128 IMPLICIT NONE
129
130 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
131 ims , ime , jms , jme , kms , kme , &
132 its , ite , jts , jte , kts , kte , &
133 iswater , isice
134 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
135 INTEGER , INTENT(IN) :: num_soil_layers
136
137 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
138 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
139 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
140 REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
141 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
142 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
143 vegcat, xland , soilcat , tmn , &
144 tsk_old , vegfra
145 REAL , INTENT(IN) :: seaice_threshold
146 REAL :: total_depth , mid_point_depth
147
148 INTEGER :: i , j , num_seaice_changes , loop
149 CHARACTER (LEN=132) :: message
150
151 num_seaice_changes = 0
152 fix_seaice : SELECT CASE ( scheme )
153
154 CASE ( SLABSCHEME )
155
156 CASE ( LSMSCHEME , RUCLSMSCHEME )
157 DO j = jts , MIN(jde-1,jte)
158 DO i = its , MIN(ide-1,ite)
159 IF ( xice(i,j) .GT. 200.0 ) THEN
160 xice(i,j) = 0.
161 num_seaice_changes = num_seaice_changes + 1
162 END IF
163 END DO
164 END DO
165 IF ( num_seaice_changes .GT. 0 ) THEN
166 WRITE ( message , FMT='(A,I6)' ) &
167 'Total post number of sea ice locations removed (due to FLAG values) = ', &
168 num_seaice_changes
169 CALL wrf_debug ( 0 , message )
170 END IF
171 num_seaice_changes = 0
172 DO j = jts , MIN(jde-1,jte)
173 DO i = its , MIN(ide-1,ite)
174 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
175 ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
176 tsk(i,j) = tsk_old(i,j)
177 END IF
178 IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
179 ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
180 print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
181 CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
182 ELSE IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
183 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
184 xice(i,j) = 1.
185 num_seaice_changes = num_seaice_changes + 1
186 if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
187 vegcat(i,j)=isice
188 ivgtyp(i,j)=isice
189 lu_index(i,j)=isice
190 landmask(i,j)=1.
191 xland(i,j)=1.
192 vegfra(i,j)=0.
193 DO loop=1,num_veg_cat
194 landusef(i,loop,j)=0.
195 END DO
196 landusef(i,ivgtyp(i,j),j)=1.
197
198 tsk_old(i,j) = tsk(i,j)
199
200 isltyp(i,j) = 16
201 soilcat(i,j)=isltyp(i,j)
202 DO loop=1,num_soil_top_cat
203 soilctop(i,loop,j)=0
204 END DO
205 DO loop=1,num_soil_bot_cat
206 soilcbot(i,loop,j)=0
207 END DO
208 soilctop(i,isltyp(i,j),j)=1.
209 soilcbot(i,isltyp(i,j),j)=1.
210
211 total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
212 DO loop = 1,num_soil_layers
213 mid_point_depth=(total_depth/num_soil_layers)/2. + &
214 (loop-1)*(total_depth/num_soil_layers)
215 tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
216 mid_point_depth*tmn(i,j) ) / total_depth
217 END DO
218
219 DO loop=1,num_soil_layers
220 smois(i,loop,j) = 1.0
221 sh2o(i,loop,j) = 0.0
222 END DO
223 ELSE IF ( xice(i,j) .LT. 0.5 ) THEN
224 xice(i,j) = 0.
225 END IF
226 END DO
227 END DO
228 IF ( num_seaice_changes .GT. 0 ) THEN
229 WRITE ( message , FMT='(A,I6)' ) &
230 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
231 CALL wrf_debug ( 0 , message )
232 END IF
233
234 END SELECT fix_seaice
235
236 END SUBROUTINE adjust_for_seaice_post
237
238 SUBROUTINE process_percent_cat_new ( landmask , &
239 landuse_frac , soil_top_cat , soil_bot_cat , &
240 isltyp , ivgtyp , &
241 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
242 ids , ide , jds , jde , kds , kde , &
243 ims , ime , jms , jme , kms , kme , &
244 its , ite , jts , jte , kts , kte , &
245 iswater )
246
247 IMPLICIT NONE
248
249 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
250 ims , ime , jms , jme , kms , kme , &
251 its , ite , jts , jte , kts , kte , &
252 iswater
253 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
254 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
255 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
256 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
257 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
258 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
259
260 INTEGER :: i , j , l , ll, dominant_index
261 REAL :: dominant_value
262
263 #ifdef WRF_CHEM
264 ! REAL :: lwthresh = .99
265 REAL :: lwthresh = .50
266 #else
267 REAL :: lwthresh = .50
268 #endif
269
270 INTEGER , PARAMETER :: iswater_soil = 14
271 INTEGER :: iforce
272 CHARACTER (LEN=132) :: message
273 integer :: change_water , change_land
274 change_water = 0
275 change_land = 0
276
277 ! Sanity check on the 50/50 points
278
279 DO j = jts , MIN(jde-1,jte)
280 DO i = its , MIN(ide-1,ite)
281 dominant_value = landuse_frac(i,iswater,j)
282 IF ( dominant_value .EQ. lwthresh ) THEN
283 DO l = 1 , num_veg_cat
284 IF ( l .EQ. iswater ) CYCLE
285 IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
286 PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
287 landuse_frac(i,l,j) = lwthresh - .01
288 landuse_frac(i,iswater,j) = lwthresh + 0.01
289 ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
290 PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
291 landuse_frac(i,l,j) = lwthresh + .01
292 landuse_frac(i,iswater,j) = lwthresh - 0.01
293 END IF
294 END DO
295 END IF
296 END DO
297 END DO
298
299 ! Compute the dominant VEGETATION INDEX.
300
301 DO j = jts , MIN(jde-1,jte)
302 DO i = its , MIN(ide-1,ite)
303 dominant_value = landuse_frac(i,1,j)
304 dominant_index = 1
305 DO l = 2 , num_veg_cat
306 IF ( l .EQ. iswater ) THEN
307 ! wait a bit
308 ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
309 dominant_value = landuse_frac(i,l,j)
310 dominant_index = l
311 END IF
312 END DO
313 IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
314 dominant_value = landuse_frac(i,iswater,j)
315 dominant_index = iswater
316 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
317 ( landmask(i,j) .LT. 0.5) .AND. &
318 ( dominant_value .EQ. lwthresh) ) THEN
319 dominant_value = landuse_frac(i,iswater,j)
320 dominant_index = iswater
321 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
322 ( landmask(i,j) .GT. 0.5) .AND. &
323 ( dominant_value .EQ. lwthresh) ) THEN
324 !no op
325 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
326 ( dominant_value .LT. lwthresh ) ) THEN
327 dominant_value = landuse_frac(i,iswater,j)
328 dominant_index = iswater
329 END IF
330 IF ( dominant_index .EQ. iswater ) THEN
331 if(landmask(i,j).gt.lwthresh) then
332 !print *,'changing to water at point ',i,j
333 !print '(24(i3,1x))',1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16, 17,18,19,20,21, 22, 23,24
334 !print '(24(i3,1x))',nint(landuse_frac(i,:,j)*100)
335 change_water=change_water+1
336 endif
337 landmask(i,j) = 0
338 ELSE IF ( dominant_index .NE. iswater ) THEN
339 if(landmask(i,j).lt.lwthresh) then
340 !print *,'changing to land at point ',i,j
341 !print '(24(i3,1x))',1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16, 17,18,19,20,21, 22, 23,24
342 !print '(24(i3,1x))',nint(landuse_frac(i,:,j)*100)
343 change_land=change_land+1
344 endif
345 landmask(i,j) = 1
346 END IF
347 ivgtyp(i,j) = dominant_index
348 END DO
349 END DO
350
351 ! Compute the dominant SOIL TEXTURE INDEX, TOP.
352
353 iforce = 0
354 DO i = its , MIN(ide-1,ite)
355 DO j = jts , MIN(jde-1,jte)
356 dominant_value = soil_top_cat(i,1,j)
357 dominant_index = 1
358 IF ( landmask(i,j) .GT. lwthresh ) THEN
359 DO l = 2 , num_soil_top_cat
360 IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
361 dominant_value = soil_top_cat(i,l,j)
362 dominant_index = l
363 END IF
364 END DO
365 IF ( dominant_value .LT. 0.01 ) THEN
366 iforce = iforce + 1
367 WRITE ( message , FMT = '(A,I4,I4)' ) &
368 'based on landuse, changing soil to land at point ',i,j
369 CALL wrf_debug(1,message)
370 WRITE ( message , FMT = '(16(i3,1x))' ) &
371 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
372 CALL wrf_debug(1,message)
373 WRITE ( message , FMT = '(16(i3,1x))' ) &
374 nint(soil_top_cat(i,:,j)*100)
375 CALL wrf_debug(1,message)
376 dominant_index = 8
377 END IF
378 ELSE
379 dominant_index = iswater_soil
380 END IF
381 isltyp(i,j) = dominant_index
382 END DO
383 END DO
384
385 if(iforce.ne.0)then
386 WRITE(message,FMT='(A,I4,A,I6)' ) &
387 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
388 (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
389 CALL wrf_debug(0,message)
390 endif
391 print *,'LAND CHANGE = ',change_land
392 print *,'WATER CHANGE = ',change_water
393
394 END SUBROUTINE process_percent_cat_new
395
396 SUBROUTINE process_soil_real ( tsk , tmn , &
397 landmask , sst , &
398 st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
399 zs , dzs , tslb , smois , sh2o , &
400 flag_sst , flag_soilt000, flag_soilm000, &
401 ids , ide , jds , jde , kds , kde , &
402 ims , ime , jms , jme , kms , kme , &
403 its , ite , jts , jte , kts , kte , &
404 sf_surface_physics , num_soil_layers , real_data_init_type , &
405 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
406 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
407
408 IMPLICIT NONE
409
410 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
411 ims , ime , jms , jme , kms , kme , &
412 its , ite , jts , jte , kts , kte , &
413 sf_surface_physics , num_soil_layers , real_data_init_type , &
414 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
415 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
416
417 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
418
419 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
420
421 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
422 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
423 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
424 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
425 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
426 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
427
428 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
429 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
430 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
431 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
432
433 INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
434 REAL :: dominant_value
435
436 ! Initialize the soil depth, and the soil temperature and moisture.
437
438 IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
439 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
440 CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
441 landmask , sst , flag_sst , &
442 ids , ide , jds , jde , kds , kde , &
443 ims , ime , jms , jme , kms , kme , &
444 its , ite , jts , jte , kts , kte )
445
446 ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
447 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
448 CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
449 st_input , sm_input , sw_input , landmask , sst , &
450 zs , dzs , &
451 st_levels_input , sm_levels_input , sw_levels_input , &
452 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
453 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
454 flag_sst , flag_soilt000 , flag_soilm000 , &
455 ids , ide , jds , jde , kds , kde , &
456 ims , ime , jms , jme , kms , kme , &
457 its , ite , jts , jte , kts , kte )
458 ! CALL init_soil_old ( tsk , tmn , &
459 ! smois , tslb , zs , dzs , num_soil_layers , &
460 ! st000010_input , st010040_input , st040100_input , st100200_input , &
461 ! st010200_input , &
462 ! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
463 ! sm010200_input , &
464 ! landmask_input , sst_input , &
465 ! ids , ide , jds , jde , kds , kde , &
466 ! ims , ime , jms , jme , kms , kme , &
467 ! its , ite , jts , jte , kts , kte )
468 ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
469 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
470 CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
471 st_input , sm_input , landmask , sst , &
472 zs , dzs , &
473 st_levels_input , sm_levels_input , &
474 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
475 num_st_levels_alloc , num_sm_levels_alloc , &
476 flag_sst , flag_soilt000 , flag_soilm000 , &
477 ids , ide , jds , jde , kds , kde , &
478 ims , ime , jms , jme , kms , kme , &
479 its , ite , jts , jte , kts , kte )
480 END IF
481
482 END SUBROUTINE process_soil_real
483
484 SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
485 ivgtyp,isltyp,tslb,smois, &
486 tsk,tmn,zs,dzs, &
487 num_soil_layers, &
488 sf_surface_physics , &
489 ids,ide, jds,jde, kds,kde,&
490 ims,ime, jms,jme, kms,kme,&
491 its,ite, jts,jte, kts,kte )
492
493 IMPLICIT NONE
494
495 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
496 ims,ime, jms,jme, kms,kme, &
497 its,ite, jts,jte, kts,kte
498
499 INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
500
501 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
502
503 REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
504
505 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
506 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
507 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
508
509 ! Local variables.
510
511 INTEGER :: itf,jtf
512
513 itf=MIN(ite,ide-1)
514 jtf=MIN(jte,jde-1)
515
516 IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
517 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
518 CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
519 ivgtyp,zs,dzs,num_soil_layers, &
520 ids,ide, jds,jde, kds,kde, &
521 ims,ime, jms,jme, kms,kme, &
522 its,ite, jts,jte, kts,kte )
523 ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
524 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
525 CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
526 ivgtyp,isltyp,tslb,smois,tmn, &
527 num_soil_layers, &
528 ids,ide, jds,jde, kds,kde, &
529 ims,ime, jms,jme, kms,kme, &
530 its,ite, jts,jte, kts,kte )
531 ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
532 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
533
534 END IF
535
536 END SUBROUTINE process_soil_ideal
537
538 SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
539 tsk , ter , toposoil , landmask , flag_toposoil , &
540 st000010 , st010040 , st040100 , st100200 , st010200 , &
541 flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
542 st000007 , st007028 , st028100 , st100255 , &
543 flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
544 soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
545 flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
546 ids , ide , jds , jde , kds , kde , &
547 ims , ime , jms , jme , kms , kme , &
548 its , ite , jts , jte , kts , kte )
549
550 IMPLICIT NONE
551
552 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
553 ims , ime , jms , jme , kms , kme , &
554 its , ite , jts , jte , kts , kte
555
556 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
557 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
558 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
559 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000007 , st007028 , st028100 , st100255
560 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
561
562 INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
563 INTEGER , INTENT(IN) :: flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255
564 INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
565 INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
566
567 INTEGER :: i , j
568
569 REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
570
571 ! Do we have a soil field with which to modify soil temperatures?
572
573 IF ( flag_toposoil .EQ. 1 ) THEN
574
575 DO j = jts , MIN(jde-1,jte)
576 DO i = its , MIN(ide-1,ite)
577
578 ! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
579 ! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
580 ! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
581 ! the difference between the soil elevation and the terrain is greater than 3 km means
582 ! that the soil data is either all zeros or that the data are inconsistent. Any of these
583 ! three conditions is grievous enough to induce a WRF fatality. However, if they are at
584 ! a water point, then we can safely ignore them.
585
586 soil_elev_min_val = toposoil(i,j)
587 soil_elev_max_val = toposoil(i,j)
588 soil_elev_min_dif = ter(i,j) - toposoil(i,j)
589 soil_elev_max_dif = ter(i,j) - toposoil(i,j)
590
591 IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
592 CYCLE
593 ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
594 !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
595 cycle
596 ! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
597 ENDIF
598
599 IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
600 CYCLE
601 ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
602 print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
603 cycle
604 CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
605 ENDIF
606
607 IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
608 ( landmask(i,j) .LT. 0.5 ) ) THEN
609 CYCLE
610 ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
611 ( landmask(i,j) .GT. 0.5 ) ) THEN
612 print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
613 cycle
614 CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
615 ENDIF
616
617 ! For each of the fields that we would like to modify, check to see if it came in from the SI.
618 ! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
619 ! are not at a water point.
620
621 IF (landmask(i,j) .GT. 0.5 ) THEN
622
623 IF ( sf_surface_physics .EQ. 1 ) THEN
624 tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
625 END IF
626
627 tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
628
629 IF ( flag_st000010 .EQ. 1 ) THEN
630 st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
631 END IF
632 IF ( flag_st010040 .EQ. 1 ) THEN
633 st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
634 END IF
635 IF ( flag_st040100 .EQ. 1 ) THEN
636 st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
637 END IF
638 IF ( flag_st100200 .EQ. 1 ) THEN
639 st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
640 END IF
641 IF ( flag_st010200 .EQ. 1 ) THEN
642 st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
643 END IF
644
645 IF ( flag_st000007 .EQ. 1 ) THEN
646 st000007(i,j) = st000007(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
647 END IF
648 IF ( flag_st007028 .EQ. 1 ) THEN
649 st007028(i,j) = st007028(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
650 END IF
651 IF ( flag_st028100 .EQ. 1 ) THEN
652 st028100(i,j) = st028100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
653 END IF
654 IF ( flag_st100255 .EQ. 1 ) THEN
655 st100255(i,j) = st100255(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
656 END IF
657
658 IF ( flag_soilt000 .EQ. 1 ) THEN
659 soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
660 END IF
661 IF ( flag_soilt005 .EQ. 1 ) THEN
662 soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
663 END IF
664 IF ( flag_soilt020 .EQ. 1 ) THEN
665 soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
666 END IF
667 IF ( flag_soilt040 .EQ. 1 ) THEN
668 soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
669 END IF
670 IF ( flag_soilt160 .EQ. 1 ) THEN
671 soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
672 END IF
673 IF ( flag_soilt300 .EQ. 1 ) THEN
674 soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
675 END IF
676
677 END IF
678 END DO
679 END DO
680
681 END IF
682
683 END SUBROUTINE adjust_soil_temp_new
684
685
686 SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
687
688 IMPLICIT NONE
689
690 INTEGER, INTENT(IN) :: num_soil_layers
691
692 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
693
694 INTEGER :: l
695
696 ! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
697 ! The distance from the ground level to the midpoint of the layer is given by zs.
698
699 ! ------- Ground Level ---------- || || || ||
700 ! || || || || zs(1) = 0.005 m
701 ! -- -- -- -- -- -- -- -- -- || || || \/
702 ! || || ||
703 ! ----------------------------------- || || || \/ dzs(1) = 0.01 m
704 ! || || ||
705 ! || || || zs(2) = 0.02
706 ! -- -- -- -- -- -- -- -- -- || || \/
707 ! || ||
708 ! || ||
709 ! ----------------------------------- || || \/ dzs(2) = 0.02 m
710 ! || ||
711 ! || ||
712 ! || ||
713 ! || || zs(3) = 0.05
714 ! -- -- -- -- -- -- -- -- -- || \/
715 ! ||
716 ! ||
717 ! ||
718 ! ||
719 ! ----------------------------------- \/ dzs(3) = 0.04 m
720
721 IF ( num_soil_layers .NE. 5 ) THEN
722 PRINT '(A)','Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
723 CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
724 END IF
725
726 dzs(1)=.01
727 zs(1)=.5*dzs(1)
728
729 DO l=2,num_soil_layers
730 dzs(l)=2*dzs(l-1)
731 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
732 ENDDO
733
734 END SUBROUTINE init_soil_depth_1
735
736 SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
737
738 IMPLICIT NONE
739
740 INTEGER, INTENT(IN) :: num_soil_layers
741
742 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
743
744 INTEGER :: l
745
746 dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
747
748 IF ( num_soil_layers .NE. 4 ) THEN
749 PRINT '(A)','Usually, the LSM uses 4 layers. Change this in the namelist.'
750 CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
751 END IF
752
753 zs(1)=.5*dzs(1)
754
755 DO l=2,num_soil_layers
756 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
757 ENDDO
758
759 END SUBROUTINE init_soil_depth_2
760
761 SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
762
763 IMPLICIT NONE
764
765 INTEGER, INTENT(IN) :: num_soil_layers
766
767 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
768
769 INTEGER :: l
770
771 CHARACTER (LEN=132) :: message
772
773 ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
774 ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
775 ! Other options with number of levels are possible, but
776 ! WRF users should change consistently the namelist entry with the
777 ! ZS array in this subroutine.
778
779 IF ( num_soil_layers .EQ. 6) THEN
780 zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
781 ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
782 ELSEIF ( num_soil_layers .EQ. 9) THEN
783 zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
784 ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
785 ENDIF
786
787 IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
788 write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
789 CALL wrf_error_fatal ( message )
790 END IF
791
792 END SUBROUTINE init_soil_depth_3
793
794 SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
795 num_soil_layers , real_data_init_type , &
796 landmask , sst , flag_sst , &
797 ids , ide , jds , jde , kds , kde , &
798 ims , ime , jms , jme , kms , kme , &
799 its , ite , jts , jte , kts , kte )
800
801 IMPLICIT NONE
802
803 INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
804 ids , ide , jds , jde , kds , kde , &
805 ims , ime , jms , jme , kms , kme , &
806 its , ite , jts , jte , kts , kte
807
808 INTEGER , INTENT(IN) :: flag_sst
809
810 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
811 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
812
813 REAL , DIMENSION(num_soil_layers) :: zs , dzs
814
815 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
816
817 INTEGER :: i , j , l
818
819 ! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
820 ! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
821 ! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
822 ! annual mean temperature.
823
824 DO j = jts , MIN(jde-1,jte)
825 DO i = its , MIN(ide-1,ite)
826 IF ( landmask(i,j) .GT. 0.5 ) THEN
827 DO l = 1 , num_soil_layers
828 tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
829 tmn(i,j) * ( zs( l) - zs(1) ) ) / &
830 ( zs(num_soil_layers) - zs(1) )
831 END DO
832 ELSE
833 IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
834 DO l = 1 , num_soil_layers
835 tslb(i,l,j)= sst(i,j)
836 END DO
837 ELSE
838 DO l = 1 , num_soil_layers
839 tslb(i,l,j)= tsk(i,j)
840 END DO
841 END IF
842 END IF
843 END DO
844 END DO
845
846 END SUBROUTINE init_soil_1_real
847
848 SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
849 ivgtyp,ZS,DZS,num_soil_layers, &
850 ids,ide, jds,jde, kds,kde, &
851 ims,ime, jms,jme, kms,kme, &
852 its,ite, jts,jte, kts,kte )
853
854 IMPLICIT NONE
855
856 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
857 ims,ime, jms,jme, kms,kme, &
858 its,ite, jts,jte, kts,kte
859
860 INTEGER, INTENT(IN ) :: num_soil_layers
861
862 REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
863 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
864 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
865
866 REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
867
868 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
869
870 ! Lcal variables.
871
872 INTEGER :: l,j,i,itf,jtf
873
874 itf=MIN(ite,ide-1)
875 jtf=MIN(jte,jde-1)
876
877 IF (num_soil_layers.NE.1)THEN
878 DO j=jts,jtf
879 DO l=1,num_soil_layers
880 DO i=its,itf
881 tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
882 ( zs(num_soil_layers)-zs(1) )
883 ENDDO
884 ENDDO
885 ENDDO
886 ENDIF
887 DO j=jts,jtf
888 DO i=its,itf
889 xland(i,j) = 2
890 ivgtyp(i,j) = 7
891 ENDDO
892 ENDDO
893
894 END SUBROUTINE init_soil_1_ideal
895
896 SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
897 st_input , sm_input , sw_input , landmask , sst , &
898 zs , dzs , &
899 st_levels_input , sm_levels_input , sw_levels_input , &
900 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
901 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
902 flag_sst , flag_soilt000 , flag_soilmt000 , &
903 ids , ide , jds , jde , kds , kde , &
904 ims , ime , jms , jme , kms , kme , &
905 its , ite , jts , jte , kts , kte )
906
907 IMPLICIT NONE
908
909 INTEGER , INTENT(IN) :: num_soil_layers , &
910 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
911 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
912 ids , ide , jds , jde , kds , kde , &
913 ims , ime , jms , jme , kms , kme , &
914 its , ite , jts , jte , kts , kte
915
916 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilmt000
917
918 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
919 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
920 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
921
922 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
923 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
924 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
925 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
926
927 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
928 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
929 REAL , DIMENSION(num_soil_layers) :: zs , dzs
930
931 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
932
933 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
934
935 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
936 REAL :: temp
937 LOGICAL :: found_levels
938
939 ! Are there any soil temp and moisture levels - ya know, they are mandatory.
940
941 num = num_st_levels_input * num_sm_levels_input
942
943 IF ( num .GE. 1 ) THEN
944
945 ! Ordered levels that we have data for.
946
947 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
948
949 ! Sort the levels for temperature.
950
951 outert : DO lout = 1 , num_st_levels_input-1
952 innert : DO lin = lout+1 , num_st_levels_input
953 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
954 temp = st_levels_input(lout)
955 st_levels_input(lout) = st_levels_input(lin)
956 st_levels_input(lin) = NINT(temp)
957 DO j = jts , MIN(jde-1,jte)
958 DO i = its , MIN(ide-1,ite)
959 temp = st_input(i,lout+1,j)
960 st_input(i,lout+1,j) = st_input(i,lin+1,j)
961 st_input(i,lin+1,j) = temp
962 END DO
963 END DO
964 END IF
965 END DO innert
966 END DO outert
967 DO j = jts , MIN(jde-1,jte)
968 DO i = its , MIN(ide-1,ite)
969 st_input(i,1,j) = tsk(i,j)
970 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
971 END DO
972 END DO
973
974 ! Sort the levels for moisture.
975
976 outerm: DO lout = 1 , num_sm_levels_input-1
977 innerm : DO lin = lout+1 , num_sm_levels_input
978 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
979 temp = sm_levels_input(lout)
980 sm_levels_input(lout) = sm_levels_input(lin)
981 sm_levels_input(lin) = NINT(temp)
982 DO j = jts , MIN(jde-1,jte)
983 DO i = its , MIN(ide-1,ite)
984 temp = sm_input(i,lout+1,j)
985 sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
986 sm_input(i,lin+1,j) = temp
987 END DO
988 END DO
989 END IF
990 END DO innerm
991 END DO outerm
992 DO j = jts , MIN(jde-1,jte)
993 DO i = its , MIN(ide-1,ite)
994 sm_input(i,1,j) = sm_input(i,2,j)
995 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
996 END DO
997 END DO
998
999 ! Sort the levels for liquid moisture.
1000
1001 outerw: DO lout = 1 , num_sw_levels_input-1
1002 innerw : DO lin = lout+1 , num_sw_levels_input
1003 IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1004 temp = sw_levels_input(lout)
1005 sw_levels_input(lout) = sw_levels_input(lin)
1006 sw_levels_input(lin) = NINT(temp)
1007 DO j = jts , MIN(jde-1,jte)
1008 DO i = its , MIN(ide-1,ite)
1009 temp = sw_input(i,lout+1,j)
1010 sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1011 sw_input(i,lin+1,j) = temp
1012 END DO
1013 END DO
1014 END IF
1015 END DO innerw
1016 END DO outerw
1017 IF ( num_sw_levels_input .GT. 1 ) THEN
1018 DO j = jts , MIN(jde-1,jte)
1019 DO i = its , MIN(ide-1,ite)
1020 sw_input(i,1,j) = sw_input(i,2,j)
1021 sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1022 END DO
1023 END DO
1024 END IF
1025
1026 found_levels = .TRUE.
1027
1028 ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
1029
1030 found_levels = .FALSE.
1031
1032 ELSE
1033 CALL wrf_error_fatal ( &
1034 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
1035 END IF
1036
1037 ! Is it OK to continue?
1038
1039 IF ( found_levels ) THEN
1040
1041 ! Here are the levels that we have from the input for temperature. The input levels plus
1042 ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1043
1044 zhave(1) = 0.
1045 DO l = 1 , num_st_levels_input
1046 zhave(l+1) = st_levels_input(l) / 100.
1047 END DO
1048 zhave(num_st_levels_input+2) = 300. / 100.
1049
1050 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1051
1052 z_wantt : DO lwant = 1 , num_soil_layers
1053 z_havet : DO lhave = 1 , num_st_levels_input +2 -1
1054 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1055 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1056 DO j = jts , MIN(jde-1,jte)
1057 DO i = its , MIN(ide-1,ite)
1058 tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1059 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1060 ( zhave(lhave+1) - zhave(lhave) )
1061 END DO
1062 END DO
1063 EXIT z_havet
1064 END IF
1065 END DO z_havet
1066 END DO z_wantt
1067
1068 ! Here are the levels that we have from the input for moisture. The input levels plus
1069 ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
1070 ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
1071 ! as the most deep layer's value.
1072
1073 zhave(1) = 0.
1074 DO l = 1 , num_sm_levels_input
1075 zhave(l+1) = sm_levels_input(l) / 100.
1076 END DO
1077 zhave(num_sm_levels_input+2) = 300. / 100.
1078
1079 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1080
1081 z_wantm : DO lwant = 1 , num_soil_layers
1082 z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
1083 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1084 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1085 DO j = jts , MIN(jde-1,jte)
1086 DO i = its , MIN(ide-1,ite)
1087 smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1088 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1089 ( zhave(lhave+1) - zhave(lhave) )
1090 END DO
1091 END DO
1092 EXIT z_havem
1093 END IF
1094 END DO z_havem
1095 END DO z_wantm
1096
1097 ! Any liquid soil moisture to worry about?
1098
1099 IF ( num_sw_levels_input .GT. 1 ) THEN
1100
1101 zhave(1) = 0.
1102 DO l = 1 , num_sw_levels_input
1103 zhave(l+1) = sw_levels_input(l) / 100.
1104 END DO
1105 zhave(num_sw_levels_input+2) = 300. / 100.
1106
1107 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1108
1109 z_wantw : DO lwant = 1 , num_soil_layers
1110 z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1111 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1112 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1113 DO j = jts , MIN(jde-1,jte)
1114 DO i = its , MIN(ide-1,ite)
1115 sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1116 sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1117 ( zhave(lhave+1) - zhave(lhave) )
1118 END DO
1119 END DO
1120 EXIT z_havew
1121 END IF
1122 END DO z_havew
1123 END DO z_wantw
1124
1125 END IF
1126
1127
1128 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
1129 ! used, but they will make a more continuous plot.
1130
1131 IF ( flag_sst .EQ. 1 ) THEN
1132 DO j = jts , MIN(jde-1,jte)
1133 DO i = its , MIN(ide-1,ite)
1134 IF ( landmask(i,j) .LT. 0.5 ) THEN
1135 DO l = 1 , num_soil_layers
1136 tslb(i,l,j)= sst(i,j)
1137 smois(i,l,j)= 1.0
1138 sh2o (i,l,j)= 1.0
1139 END DO
1140 END IF
1141 END DO
1142 END DO
1143 ELSE
1144 DO j = jts , MIN(jde-1,jte)
1145 DO i = its , MIN(ide-1,ite)
1146 IF ( landmask(i,j) .LT. 0.5 ) THEN
1147 DO l = 1 , num_soil_layers
1148 tslb(i,l,j)= tsk(i,j)
1149 smois(i,l,j)= 1.0
1150 sh2o (i,l,j)= 1.0
1151 END DO
1152 END IF
1153 END DO
1154 END DO
1155 END IF
1156
1157 DEALLOCATE (zhave)
1158
1159 END IF
1160
1161 END SUBROUTINE init_soil_2_real
1162
1163 SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
1164 ivgtyp,isltyp,tslb,smois,tmn, &
1165 num_soil_layers, &
1166 ids,ide, jds,jde, kds,kde, &
1167 ims,ime, jms,jme, kms,kme, &
1168 its,ite, jts,jte, kts,kte )
1169
1170 IMPLICIT NONE
1171
1172 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
1173 ims,ime, jms,jme, kms,kme, &
1174 its,ite, jts,jte, kts,kte
1175
1176 INTEGER, INTENT(IN) ::num_soil_layers
1177
1178 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
1179
1180 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
1181
1182 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1183
1184 INTEGER :: icm,jcm,itf,jtf
1185 INTEGER :: i,j,l
1186
1187 itf=min0(ite,ide-1)
1188 jtf=min0(jte,jde-1)
1189
1190 icm = ide/2
1191 jcm = jde/2
1192
1193 DO j=jts,jtf
1194 DO l=1,num_soil_layers
1195 DO i=its,itf
1196
1197 smois(i,1,j)=0.10
1198 smois(i,2,j)=0.10
1199 smois(i,3,j)=0.10
1200 smois(i,4,j)=0.10
1201
1202 tslb(i,1,j)=295.
1203 tslb(i,2,j)=297.
1204 tslb(i,3,j)=293.
1205 tslb(i,4,j)=293.
1206
1207 ENDDO
1208 ENDDO
1209 ENDDO
1210
1211 DO j=jts,jtf
1212 DO i=its,itf
1213 xland(i,j) = 2
1214 tmn(i,j) = 294.
1215 xice(i,j) = 0.
1216 vegfra(i,j) = 0.
1217 snow(i,j) = 0.
1218 canwat(i,j) = 0.
1219 ivgtyp(i,j) = 7
1220 isltyp(i,j) = 8
1221 ENDDO
1222 ENDDO
1223
1224 END SUBROUTINE init_soil_2_ideal
1225
1226 SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
1227 st_input , sm_input , landmask, sst, &
1228 zs , dzs , &
1229 st_levels_input , sm_levels_input , &
1230 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
1231 num_st_levels_alloc , num_sm_levels_alloc , &
1232 flag_sst , flag_soilt000 , flag_soilm000 , &
1233 ids , ide , jds , jde , kds , kde , &
1234 ims , ime , jms , jme , kms , kme , &
1235 its , ite , jts , jte , kts , kte )
1236
1237 IMPLICIT NONE
1238
1239 INTEGER , INTENT(IN) :: num_soil_layers , &
1240 num_st_levels_input , num_sm_levels_input , &
1241 num_st_levels_alloc , num_sm_levels_alloc , &
1242 ids , ide , jds , jde , kds , kde , &
1243 ims , ime , jms , jme , kms , kme , &
1244 its , ite , jts , jte , kts , kte
1245
1246 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
1247
1248 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1249 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1250
1251 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1252 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1253 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1254
1255 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1256 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1257 REAL , DIMENSION(num_soil_layers) :: zs , dzs
1258
1259 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
1260
1261 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1262
1263 INTEGER :: i , j , l , lout , lin , lwant , lhave
1264 REAL :: temp
1265
1266 CHARACTER (LEN=132) :: message
1267
1268 ! Allocate the soil layer array used for interpolating.
1269
1270 IF ( ( num_st_levels_input .LE. 0 ) .OR. &
1271 ( num_sm_levels_input .LE. 0 ) ) THEN
1272 write (message, FMT='(A)')&
1273 'No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
1274 CALL wrf_error_fatal ( message )
1275 ELSE
1276 IF ( flag_soilt000 .eq. 1 ) THEN
1277 write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1278 CALL wrf_message ( message )
1279 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
1280 ELSE
1281 write(message, FMT='(A)') ' Assume non-RUC LSM input'
1282 CALL wrf_message ( message )
1283 ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
1284 END IF
1285 END IF
1286
1287 ! Sort the levels for temperature.
1288
1289 outert : DO lout = 1 , num_st_levels_input-1
1290 innert : DO lin = lout+1 , num_st_levels_input
1291 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1292 temp = st_levels_input(lout)
1293 st_levels_input(lout) = st_levels_input(lin)
1294 st_levels_input(lin) = NINT(temp)
1295 DO j = jts , MIN(jde-1,jte)
1296 DO i = its , MIN(ide-1,ite)
1297 temp = st_input(i,lout,j)
1298 st_input(i,lout,j) = st_input(i,lin,j)
1299 st_input(i,lin,j) = temp
1300 END DO
1301 END DO
1302 END IF
1303 END DO innert
1304 END DO outert
1305
1306 IF ( flag_soilt000 .NE. 1 ) THEN
1307 DO j = jts , MIN(jde-1,jte)
1308 DO i = its , MIN(ide-1,ite)
1309 st_input(i,1,j) = tsk(i,j)
1310 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1311 END DO
1312 END DO
1313 END IF
1314
1315 ! Sort the levels for moisture.
1316
1317 outerm: DO lout = 1 , num_sm_levels_input-1
1318 innerm : DO lin = lout+1 , num_sm_levels_input
1319 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1320 temp = sm_levels_input(lout)
1321 sm_levels_input(lout) = sm_levels_input(lin)
1322 sm_levels_input(lin) = NINT(temp)
1323 DO j = jts , MIN(jde-1,jte)
1324 DO i = its , MIN(ide-1,ite)
1325 temp = sm_input(i,lout,j)
1326 sm_input(i,lout,j) = sm_input(i,lin,j)
1327 sm_input(i,lin,j) = temp
1328 END DO
1329 END DO
1330 END IF
1331 END DO innerm
1332 END DO outerm
1333
1334 IF ( flag_soilm000 .NE. 1 ) THEN
1335 DO j = jts , MIN(jde-1,jte)
1336 DO i = its , MIN(ide-1,ite)
1337 sm_input(i,1,j) = sm_input(i,2,j)
1338 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1339 END DO
1340 END DO
1341 END IF
1342
1343 ! Here are the levels that we have from the input for temperature.
1344
1345 IF ( flag_soilt000 .EQ. 1 ) THEN
1346 DO l = 1 , num_st_levels_input
1347 zhave(l) = st_levels_input(l) / 100.
1348 END DO
1349
1350 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1351
1352 z_wantt : DO lwant = 1 , num_soil_layers
1353 z_havet : DO lhave = 1 , num_st_levels_input -1
1354 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1355 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1356 DO j = jts , MIN(jde-1,jte)
1357 DO i = its , MIN(ide-1,ite)
1358 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1359 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1360 ( zhave(lhave+1) - zhave(lhave) )
1361 END DO
1362 END DO
1363 EXIT z_havet
1364 END IF
1365 END DO z_havet
1366 END DO z_wantt
1367
1368 ELSE
1369
1370 zhave(1) = 0.
1371 DO l = 1 , num_st_levels_input
1372 zhave(l+1) = st_levels_input(l) / 100.
1373 END DO
1374 zhave(num_st_levels_input+2) = 300. / 100.
1375
1376 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1377
1378 z_wantt_2 : DO lwant = 1 , num_soil_layers
1379 z_havet_2 : DO lhave = 1 , num_st_levels_input +2
1380 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1381 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1382 DO j = jts , MIN(jde-1,jte)
1383 DO i = its , MIN(ide-1,ite)
1384 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1385 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1386 ( zhave(lhave+1) - zhave(lhave) )
1387 END DO
1388 END DO
1389 EXIT z_havet_2
1390 END IF
1391 END DO z_havet_2
1392 END DO z_wantt_2
1393
1394 END IF
1395
1396 ! Here are the levels that we have from the input for moisture.
1397
1398 IF ( flag_soilm000 .EQ. 1 ) THEN
1399 DO l = 1 , num_sm_levels_input
1400 zhave(l) = sm_levels_input(l) / 100.
1401 END DO
1402
1403 ! Interpolate between the layers we have (zhave) and those that we want (zs).
1404
1405 z_wantm : DO lwant = 1 , num_soil_layers
1406 z_havem : DO lhave = 1 , num_sm_levels_input -1
1407 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1408 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1409 DO j = jts , MIN(jde-1,jte)
1410 DO i = its , MIN(ide-1,ite)
1411 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1412 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1413 ( zhave(lhave+1) - zhave(lhave) )
1414 END DO
1415 END DO
1416 EXIT z_havem
1417 END IF
1418 END DO z_havem
1419 END DO z_wantm
1420
1421 ELSE
1422
1423 zhave(1) = 0.
1424 DO l = 1 , num_sm_levels_input
1425 zhave(l+1) = sm_levels_input(l) / 100.
1426 END DO
1427 zhave(num_sm_levels_input+2) = 300. / 100.
1428
1429 z_wantm_2 : DO lwant = 1 , num_soil_layers
1430 z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
1431 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1432 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1433 DO j = jts , MIN(jde-1,jte)
1434 DO i = its , MIN(ide-1,ite)
1435 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1436 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1437 ( zhave(lhave+1) - zhave(lhave) )
1438 END DO
1439 END DO
1440 EXIT z_havem_2
1441 END IF
1442 END DO z_havem_2
1443 END DO z_wantm_2
1444
1445 END IF
1446 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
1447 ! used, but they will make a more continuous plot.
1448
1449 IF ( flag_sst .EQ. 1 ) THEN
1450 DO j = jts , MIN(jde-1,jte)
1451 DO i = its , MIN(ide-1,ite)
1452 IF ( landmask(i,j) .LT. 0.5 ) THEN
1453 DO l = 1 , num_soil_layers
1454 tslb(i,l,j) = sst(i,j)
1455 tsk(i,j) = sst(i,j)
1456 smois(i,l,j)= 1.0
1457 END DO
1458 END IF
1459 END DO
1460 END DO
1461 ELSE
1462 DO j = jts , MIN(jde-1,jte)
1463 DO i = its , MIN(ide-1,ite)
1464 IF ( landmask(i,j) .LT. 0.5 ) THEN
1465 DO l = 1 , num_soil_layers
1466 tslb(i,l,j)= tsk(i,j)
1467 smois(i,l,j)= 1.0
1468 END DO
1469 END IF
1470 END DO
1471 END DO
1472 END IF
1473
1474 DEALLOCATE (zhave)
1475
1476 END SUBROUTINE init_soil_3_real
1477
1478 END MODULE module_soil_pre
1479
1480 #else
1481
1482 MODULE module_soil_pre
1483
1484 USE module_date_time
1485 USE module_state_description
1486
1487 CONTAINS
1488
1489 SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1490 xland , landusef , isltyp , soilcat , soilctop , &
1491 soilcbot , tmn , &
1492 seaice_threshold , &
1493 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1494 iswater , isice , &
1495 scheme , &
1496 ids , ide , jds , jde , kds , kde , &
1497 ims , ime , jms , jme , kms , kme , &
1498 its , ite , jts , jte , kts , kte )
1499
1500 IMPLICIT NONE
1501
1502 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1503 ims , ime , jms , jme , kms , kme , &
1504 its , ite , jts , jte , kts , kte , &
1505 iswater , isice
1506 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
1507
1508 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
1509 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
1510 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
1511 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1512 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
1513 vegcat, xland , soilcat , tmn
1514 REAL , INTENT(IN) :: seaice_threshold
1515
1516 INTEGER :: i , j , num_seaice_changes , loop
1517 CHARACTER (LEN=132) :: message
1518
1519 fix_seaice : SELECT CASE ( scheme )
1520
1521 CASE ( SLABSCHEME )
1522 DO j = jts , MIN(jde-1,jte)
1523 DO i = its , MIN(ide-1,ite)
1524 IF ( xice(i,j) .GT. 200.0 ) THEN
1525 xice(i,j) = 0.
1526 num_seaice_changes = num_seaice_changes + 1
1527 END IF
1528 END DO
1529 END DO
1530 IF ( num_seaice_changes .GT. 0 ) THEN
1531 WRITE ( message , FMT='(A,I6)' ) &
1532 'Total pre number of sea ice locations removed (due to FLAG values) = ', &
1533 num_seaice_changes
1534 CALL wrf_debug ( 0 , message )
1535 END IF
1536 num_seaice_changes = 0
1537 DO j = jts , MIN(jde-1,jte)
1538 DO i = its , MIN(ide-1,ite)
1539 IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
1540 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
1541 xice(i,j) = 1.
1542 num_seaice_changes = num_seaice_changes + 1
1543 tmn(i,j) = 271.4
1544 vegcat(i,j)=isice
1545 lu_index(i,j)=ivgtyp(i,j)
1546 landmask(i,j)=1.
1547 xland(i,j)=1.
1548 DO loop=1,num_veg_cat
1549 landusef(i,loop,j)=0.
1550 END DO
1551 landusef(i,ivgtyp(i,j),j)=1.
1552
1553 isltyp(i,j) = 16
1554 soilcat(i,j)=isltyp(i,j)
1555 DO loop=1,num_soil_top_cat
1556 soilctop(i,loop,j)=0
1557 END DO
1558 DO loop=1,num_soil_bot_cat
1559 soilcbot(i,loop,j)=0
1560 END DO
1561 soilctop(i,isltyp(i,j),j)=1.
1562 soilcbot(i,isltyp(i,j),j)=1.
1563 END IF
1564 END DO
1565 END DO
1566 IF ( num_seaice_changes .GT. 0 ) THEN
1567 WRITE ( message , FMT='(A,I6)' ) &
1568 'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
1569 CALL wrf_debug ( 0 , message )
1570 END IF
1571
1572 CASE ( LSMSCHEME )
1573 CASE ( RUCLSMSCHEME )
1574
1575 END SELECT fix_seaice
1576
1577 END SUBROUTINE adjust_for_seaice_pre
1578
1579 SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1580 xland , landusef , isltyp , soilcat , soilctop , &
1581 soilcbot , tmn , &
1582 tslb , smois , sh2o , &
1583 seaice_threshold , &
1584 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1585 num_soil_layers , &
1586 iswater , isice , &
1587 scheme , &
1588 ids , ide , jds , jde , kds , kde , &
1589 ims , ime , jms , jme , kms , kme , &
1590 its , ite , jts , jte , kts , kte )
1591
1592 IMPLICIT NONE
1593
1594 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1595 ims , ime , jms , jme , kms , kme , &
1596 its , ite , jts , jte , kts , kte , &
1597 iswater , isice
1598 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
1599 INTEGER , INTENT(IN) :: num_soil_layers
1600
1601 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
1602 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
1603 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
1604 REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
1605 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1606 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
1607 vegcat, xland , soilcat , tmn
1608 REAL , INTENT(IN) :: seaice_threshold
1609 REAL :: total_depth , mid_point_depth
1610
1611 INTEGER :: i , j , num_seaice_changes , loop
1612 CHARACTER (LEN=132) :: message
1613
1614 fix_seaice : SELECT CASE ( scheme )
1615
1616 CASE ( SLABSCHEME )
1617
1618 CASE ( LSMSCHEME )
1619 DO j = jts , MIN(jde-1,jte)
1620 DO i = its , MIN(ide-1,ite)
1621 IF ( xice(i,j) .GT. 200.0 ) THEN
1622 xice(i,j) = 0.
1623 num_seaice_changes = num_seaice_changes + 1
1624 END IF
1625 END DO
1626 END DO
1627 IF ( num_seaice_changes .GT. 0 ) THEN
1628 WRITE ( message , FMT='(A,I6)' ) &
1629 'Total post number of sea ice locations removed (due to FLAG values) = ', &
1630 num_seaice_changes
1631 CALL wrf_debug ( 0 , message )
1632 END IF
1633 num_seaice_changes = 0
1634 DO j = jts , MIN(jde-1,jte)
1635 DO i = its , MIN(ide-1,ite)
1636 IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
1637 ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
1638 xice(i,j) = 1.
1639 num_seaice_changes = num_seaice_changes + 1
1640 tmn(i,j) = 271.16
1641 vegcat(i,j)=isice
1642 lu_index(i,j)=ivgtyp(i,j)
1643 landmask(i,j)=1.
1644 xland(i,j)=1.
1645 DO loop=1,num_veg_cat
1646 landusef(i,loop,j)=0.
1647 END DO
1648 landusef(i,ivgtyp(i,j),j)=1.
1649
1650 isltyp(i,j) = 16
1651 soilcat(i,j)=isltyp(i,j)
1652 DO loop=1,num_soil_top_cat
1653 soilctop(i,loop,j)=0
1654 END DO
1655 DO loop=1,num_soil_bot_cat
1656 soilcbot(i,loop,j)=0
1657 END DO
1658 soilctop(i,isltyp(i,j),j)=1.
1659 soilcbot(i,isltyp(i,j),j)=1.
1660
1661 total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
1662 DO loop = 1,num_soil_layers
1663 mid_point_depth=(total_depth/num_soil_layers)/2. + &
1664 (loop-1)*(total_depth/num_soil_layers)
1665 tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
1666 mid_point_depth*tmn(i,j) ) / total_depth
1667 END DO
1668
1669 DO loop=1,num_soil_layers
1670 smois(i,loop,j) = 1.0
1671 sh2o(i,loop,j) = 0.0
1672 END DO
1673 END IF
1674 END DO
1675 END DO
1676 IF ( num_seaice_changes .GT. 0 ) THEN
1677 WRITE ( message , FMT='(A,I6)' ) &
1678 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
1679 CALL wrf_debug ( 0 , message )
1680 END IF
1681
1682 CASE ( RUCLSMSCHEME )
1683
1684 END SELECT fix_seaice
1685
1686 END SUBROUTINE adjust_for_seaice_post
1687
1688 SUBROUTINE process_percent_cat_new ( landmask , &
1689 landuse_frac , soil_top_cat , soil_bot_cat , &
1690 isltyp , ivgtyp , &
1691 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1692 ids , ide , jds , jde , kds , kde , &
1693 ims , ime , jms , jme , kms , kme , &
1694 its , ite , jts , jte , kts , kte , &
1695 iswater )
1696
1697 IMPLICIT NONE
1698
1699 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1700 ims , ime , jms , jme , kms , kme , &
1701 its , ite , jts , jte , kts , kte , &
1702 iswater
1703 INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
1704 REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
1705 REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
1706 REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
1707 INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1708 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
1709
1710 INTEGER :: i , j , l , ll, dominant_index
1711 REAL :: dominant_value
1712
1713 #ifdef WRF_CHEM
1714 ! REAL :: lwthresh = .99
1715 REAL :: lwthresh = .50
1716 #else
1717 REAL :: lwthresh = .50
1718 #endif
1719
1720 INTEGER , PARAMETER :: iswater_soil = 14
1721 INTEGER :: iforce
1722 CHARACTER (LEN=132) :: message
1723
1724 ! Sanity check on the 50/50 points
1725
1726 DO j = jts , MIN(jde-1,jte)
1727 DO i = its , MIN(ide-1,ite)
1728 dominant_value = landuse_frac(i,iswater,j)
1729 IF ( dominant_value .EQ. lwthresh ) THEN
1730 DO l = 1 , num_veg_cat
1731 IF ( l .EQ. iswater ) CYCLE
1732 IF ( landuse_frac(i,l,j) .EQ. lwthresh ) THEN
1733 write(message, FMT='(I4,I4,A,I4,A,I4)') i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
1734 call wrf_message ( message )
1735 landuse_frac(i,l,j) = lwthresh - .01
1736 landuse_frac(i,iswater,j) = 1. - (lwthresh + 0.01)
1737 END IF
1738 END DO
1739 END IF
1740 END DO
1741 END DO
1742
1743 ! Compute the dominant VEGETATION INDEX.
1744
1745 DO j = jts , MIN(jde-1,jte)
1746 DO i = its , MIN(ide-1,ite)
1747 dominant_value = landuse_frac(i,1,j)
1748 dominant_index = 1
1749 DO l = 2 , num_veg_cat
1750 IF ( l .EQ. iswater ) THEN
1751 ! wait a bit
1752 ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
1753 dominant_value = landuse_frac(i,l,j)
1754 dominant_index = l
1755 END IF
1756 END DO
1757 IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
1758 dominant_value = landuse_frac(i,iswater,j)
1759 dominant_index = iswater
1760 #if 0
1761 si needs to beef up consistency checks before we can use this part
1762 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
1763 ( dominant_value .EQ. lwthresh) ) THEN
1764 ! no op
1765 #else
1766 ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.LT.iswater))THEN
1767 write(message,*)'temporary SI landmask fix'
1768 call wrf_message(trim(message))
1769 ! no op
1770 ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.GT.iswater))THEN
1771 write(message,*)'temporary SI landmask fix'
1772 call wrf_message(trim(message))
1773 dominant_value = landuse_frac(i,iswater,j)
1774 dominant_index = iswater
1775 #endif
1776 ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
1777 ( dominant_value .LT. lwthresh ) ) THEN
1778 dominant_value = landuse_frac(i,iswater,j)
1779 dominant_index = iswater
1780 END IF
1781 IF ( dominant_index .EQ. iswater ) THEN
1782 if(landmask(i,j).gt.lwthresh) then
1783 write(message,*) 'changing to water at point ',i,j
1784 call wrf_message(trim(message))
1785 write(message,*) nint(landuse_frac(i,:,j)*100)
1786 call wrf_message(trim(message))
1787 endif
1788 landmask(i,j) = 0
1789 ELSE IF ( dominant_index .NE. iswater ) THEN
1790 if(landmask(i,j).lt.lwthresh) then
1791 write(message,*) 'changing to land at point ',i,j
1792 call wrf_message(trim(message))
1793 write(message,*) nint(landuse_frac(i,:,j)*100)
1794 call wrf_message(trim(message))
1795 endif
1796 landmask(i,j) = 1
1797 END IF
1798 ivgtyp(i,j) = dominant_index
1799 END DO
1800 END DO
1801
1802 ! Compute the dominant SOIL TEXTURE INDEX, TOP.
1803
1804 iforce = 0
1805 DO i = its , MIN(ide-1,ite)
1806 DO j = jts , MIN(jde-1,jte)
1807 dominant_value = soil_top_cat(i,1,j)
1808 dominant_index = 1
1809 IF ( landmask(i,j) .GT. lwthresh ) THEN
1810 DO l = 2 , num_soil_top_cat
1811 IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
1812 dominant_value = soil_top_cat(i,l,j)
1813 dominant_index = l
1814 END IF
1815 END DO
1816 IF ( dominant_value .LT. 0.01 ) THEN
1817 iforce = iforce + 1
1818 WRITE ( message , FMT = '(A,I4,I4)' ) &
1819 'based on landuse, changing soil to land at point ',i,j
1820 CALL wrf_debug(1,message)
1821 WRITE ( message , FMT = '(16(i3,1x))' ) &
1822 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
1823 CALL wrf_debug(1,message)
1824 WRITE ( message , FMT = '(16(i3,1x))' ) &
1825 nint(soil_top_cat(i,:,j)*100)
1826 CALL wrf_debug(1,message)
1827 dominant_index = 8
1828 END IF
1829 ELSE
1830 dominant_index = iswater_soil
1831 END IF
1832 isltyp(i,j) = dominant_index
1833 END DO
1834 END DO
1835
1836 if(iforce.ne.0)then
1837 WRITE(message,FMT='(A,I4,A,I6)' ) &
1838 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
1839 (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
1840 CALL wrf_debug(0,message)
1841 endif
1842
1843 END SUBROUTINE process_percent_cat_new
1844
1845 SUBROUTINE process_soil_real ( tsk , tmn , &
1846 landmask , sst , &
1847 st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
1848 zs , dzs , tslb , smois , sh2o , &
1849 flag_sst , flag_soilt000, flag_soilm000, &
1850 ids , ide , jds , jde , kds , kde , &
1851 ims , ime , jms , jme , kms , kme , &
1852 its , ite , jts , jte , kts , kte , &
1853 sf_surface_physics , num_soil_layers , real_data_init_type , &
1854 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1855 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
1856
1857 IMPLICIT NONE
1858
1859 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1860 ims , ime , jms , jme , kms , kme , &
1861 its , ite , jts , jte , kts , kte , &
1862 sf_surface_physics , num_soil_layers , real_data_init_type , &
1863 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1864 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
1865
1866 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
1867
1868 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1869
1870 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1871 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1872 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1873 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1874 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1875 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1876
1877 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1878 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1879 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1880 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1881
1882 INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
1883 REAL :: dominant_value
1884
1885 ! Initialize the soil depth, and the soil temperature and moisture.
1886
1887 IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1888 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
1889 CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
1890 landmask , sst , flag_sst , &
1891 ids , ide , jds , jde , kds , kde , &
1892 ims , ime , jms , jme , kms , kme , &
1893 its , ite , jts , jte , kts , kte )
1894
1895 ELSE IF ( ( sf_surface_physics .EQ. 2 .or. sf_surface_physics .EQ. 99 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1896 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
1897 CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
1898 st_input , sm_input , sw_input , landmask , sst , &
1899 zs , dzs , &
1900 st_levels_input , sm_levels_input , sw_levels_input , &
1901 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1902 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1903 flag_sst , flag_soilt000 , flag_soilm000 , &
1904 ids , ide , jds , jde , kds , kde , &
1905 ims , ime , jms , jme , kms , kme , &
1906 its , ite , jts , jte , kts , kte )
1907 ! CALL init_soil_old ( tsk , tmn , &
1908 ! smois , tslb , zs , dzs , num_soil_layers , &
1909 ! st000010_input , st010040_input , st040100_input , st100200_input , &
1910 ! st010200_input , &
1911 ! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
1912 ! sm010200_input , &
1913 ! landmask_input , sst_input , &
1914 ! ids , ide , jds , jde , kds , kde , &
1915 ! ims , ime , jms , jme , kms , kme , &
1916 ! its , ite , jts , jte , kts , kte )
1917 ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1918 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
1919 CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
1920 st_input , sm_input , landmask , sst , &
1921 zs , dzs , &
1922 st_levels_input , sm_levels_input , &
1923 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
1924 num_st_levels_alloc , num_sm_levels_alloc , &
1925 flag_sst , flag_soilt000 , flag_soilm000 , &
1926 ids , ide , jds , jde , kds , kde , &
1927 ims , ime , jms , jme , kms , kme , &
1928 its , ite , jts , jte , kts , kte )
1929 END IF
1930
1931 END SUBROUTINE process_soil_real
1932
1933 SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
1934 ivgtyp,isltyp,tslb,smois, &
1935 tsk,tmn,zs,dzs, &
1936 num_soil_layers, &
1937 sf_surface_physics , &
1938 ids,ide, jds,jde, kds,kde,&
1939 ims,ime, jms,jme, kms,kme,&
1940 its,ite, jts,jte, kts,kte )
1941
1942 IMPLICIT NONE
1943
1944 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
1945 ims,ime, jms,jme, kms,kme, &
1946 its,ite, jts,jte, kts,kte
1947
1948 INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
1949
1950 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
1951
1952 REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
1953
1954 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
1955 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
1956 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1957
1958 ! Local variables.
1959
1960 INTEGER :: itf,jtf
1961
1962 itf=MIN(ite,ide-1)
1963 jtf=MIN(jte,jde-1)
1964
1965 IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1966 CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
1967 CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
1968 ivgtyp,zs,dzs,num_soil_layers, &
1969 ids,ide, jds,jde, kds,kde, &
1970 ims,ime, jms,jme, kms,kme, &
1971 its,ite, jts,jte, kts,kte )
1972 ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1973 CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
1974 CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
1975 ivgtyp,isltyp,tslb,smois,tmn, &
1976 num_soil_layers, &
1977 ids,ide, jds,jde, kds,kde, &
1978 ims,ime, jms,jme, kms,kme, &
1979 its,ite, jts,jte, kts,kte )
1980 ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
1981 CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
1982
1983 END IF
1984
1985 END SUBROUTINE process_soil_ideal
1986
1987 SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
1988 tsk , ter , toposoil , landmask , flag_toposoil , &
1989 st000010 , st010040 , st040100 , st100200 , st010200 , &
1990 flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
1991 soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
1992 flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
1993 ids , ide , jds , jde , kds , kde , &
1994 ims , ime , jms , jme , kms , kme , &
1995 its , ite , jts , jte , kts , kte )
1996
1997 IMPLICIT NONE
1998
1999 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2000 ims , ime , jms , jme , kms , kme , &
2001 its , ite , jts , jte , kts , kte
2002
2003 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
2004 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
2005 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
2006 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
2007
2008 INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
2009 INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
2010 INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
2011
2012 INTEGER :: i , j
2013
2014 REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
2015
2016 ! Do we have a soil field with which to modify soil temperatures?
2017
2018 IF ( flag_toposoil .EQ. 1 ) THEN
2019
2020 DO j = jts , MIN(jde-1,jte)
2021 DO i = its , MIN(ide-1,ite)
2022
2023 ! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
2024 ! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
2025 ! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
2026 ! the difference between the soil elevation and the terrain is greater than 3 km means
2027 ! that the soil data is either all zeros or that the data are inconsistent. Any of these
2028 ! three conditions is grievous enough to induce a WRF fatality. However, if they are at
2029 ! a water point, then we can safely ignore them.
2030
2031 soil_elev_min_val = toposoil(i,j)
2032 soil_elev_max_val = toposoil(i,j)
2033 soil_elev_min_dif = ter(i,j) - toposoil(i,j)
2034 soil_elev_max_dif = ter(i,j) - toposoil(i,j)
2035
2036 IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2037 CYCLE
2038 ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2039 !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
2040 cycle
2041 ! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
2042 ENDIF
2043
2044 IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2045 CYCLE
2046 ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2047 print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
2048 cycle
2049 CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
2050 ENDIF
2051
2052 IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2053 ( landmask(i,j) .LT. 0.5 ) ) THEN
2054 CYCLE
2055 ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2056 ( landmask(i,j) .GT. 0.5 ) ) THEN
2057 print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
2058 cycle
2059 CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
2060 ENDIF
2061
2062 ! For each of the fields that we would like to modify, check to see if it came in from the SI.
2063 ! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
2064 ! are not at a water point.
2065
2066 IF (landmask(i,j) .GT. 0.5 ) THEN
2067
2068 IF ( sf_surface_physics .EQ. 1 ) THEN
2069 tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2070 END IF
2071
2072 tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2073
2074 IF ( flag_st000010 .EQ. 1 ) THEN
2075 st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2076 END IF
2077 IF ( flag_st010040 .EQ. 1 ) THEN
2078 st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2079 END IF
2080 IF ( flag_st040100 .EQ. 1 ) THEN
2081 st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2082 END IF
2083 IF ( flag_st100200 .EQ. 1 ) THEN
2084 st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2085 END IF
2086 IF ( flag_st010200 .EQ. 1 ) THEN
2087 st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2088 END IF
2089
2090 IF ( flag_soilt000 .EQ. 1 ) THEN
2091 soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2092 END IF
2093 IF ( flag_soilt005 .EQ. 1 ) THEN
2094 soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2095 END IF
2096 IF ( flag_soilt020 .EQ. 1 ) THEN
2097 soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2098 END IF
2099 IF ( flag_soilt040 .EQ. 1 ) THEN
2100 soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2101 END IF
2102 IF ( flag_soilt160 .EQ. 1 ) THEN
2103 soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2104 END IF
2105 IF ( flag_soilt300 .EQ. 1 ) THEN
2106 soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2107 END IF
2108
2109 END IF
2110 END DO
2111 END DO
2112
2113 END IF
2114
2115 END SUBROUTINE adjust_soil_temp_new
2116
2117
2118 SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
2119
2120 IMPLICIT NONE
2121
2122 INTEGER, INTENT(IN) :: num_soil_layers
2123
2124 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2125
2126 INTEGER :: l
2127
2128 CHARACTER (LEN=132) :: message
2129
2130 ! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
2131 ! The distance from the ground level to the midpoint of the layer is given by zs.
2132
2133 ! ------- Ground Level ---------- || || || ||
2134 ! || || || || zs(1) = 0.005 m
2135 ! -- -- -- -- -- -- -- -- -- || || || \/
2136 ! || || ||
2137 ! ----------------------------------- || || || \/ dzs(1) = 0.01 m
2138 ! || || ||
2139 ! || || || zs(2) = 0.02
2140 ! -- -- -- -- -- -- -- -- -- || || \/
2141 ! || ||
2142 ! || ||
2143 ! ----------------------------------- || || \/ dzs(2) = 0.02 m
2144 ! || ||
2145 ! || ||
2146 ! || ||
2147 ! || || zs(3) = 0.05
2148 ! -- -- -- -- -- -- -- -- -- || \/
2149 ! ||
2150 ! ||
2151 ! ||
2152 ! ||
2153 ! ----------------------------------- \/ dzs(3) = 0.04 m
2154
2155 IF ( num_soil_layers .NE. 5 ) THEN
2156 write(message,FMT= '(A)') 'Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
2157 CALL wrf_error_fatal ( message )
2158 END IF
2159
2160 dzs(1)=.01
2161 zs(1)=.5*dzs(1)
2162
2163 DO l=2,num_soil_layers
2164 dzs(l)=2*dzs(l-1)
2165 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
2166 ENDDO
2167
2168 END SUBROUTINE init_soil_depth_1
2169
2170 SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
2171
2172 IMPLICIT NONE
2173
2174 INTEGER, INTENT(IN) :: num_soil_layers
2175
2176 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2177
2178 INTEGER :: l
2179
2180 CHARACTER (LEN=132) :: message
2181
2182 dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
2183
2184 IF ( num_soil_layers .NE. 4 ) THEN
2185 write(message,FMT='(A)') 'Usually, the LSM uses 4 layers. Change this in the namelist.'
2186 CALL wrf_error_fatal ( message )
2187 END IF
2188
2189 zs(1)=.5*dzs(1)
2190
2191 DO l=2,num_soil_layers
2192 zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
2193 ENDDO
2194
2195 END SUBROUTINE init_soil_depth_2
2196
2197 SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
2198
2199 IMPLICIT NONE
2200
2201 INTEGER, INTENT(IN) :: num_soil_layers
2202
2203 REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2204
2205 INTEGER :: l
2206
2207 CHARACTER (LEN=132) :: message
2208
2209 ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
2210 ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
2211 ! Other options with number of levels are possible, but
2212 ! WRF users should change consistently the namelist entry with the
2213 ! ZS array in this subroutine.
2214
2215 IF ( num_soil_layers .EQ. 6) THEN
2216 zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
2217 ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2218 ELSEIF ( num_soil_layers .EQ. 9) THEN
2219 zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
2220 ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2221 ENDIF
2222
2223 IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
2224 WRITE(message,FMT= '(A)')'Usually, the RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
2225 CALL wrf_error_fatal ( message )
2226 END IF
2227
2228 END SUBROUTINE init_soil_depth_3
2229
2230 SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
2231 num_soil_layers , real_data_init_type , &
2232 landmask , sst , flag_sst , &
2233 ids , ide , jds , jde , kds , kde , &
2234 ims , ime , jms , jme , kms , kme , &
2235 its , ite , jts , jte , kts , kte )
2236
2237 IMPLICIT NONE
2238
2239 INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
2240 ids , ide , jds , jde , kds , kde , &
2241 ims , ime , jms , jme , kms , kme , &
2242 its , ite , jts , jte , kts , kte
2243
2244 INTEGER , INTENT(IN) :: flag_sst
2245
2246 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2247 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
2248
2249 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2250
2251 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
2252
2253 INTEGER :: i , j , l
2254
2255 ! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
2256 ! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
2257 ! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
2258 ! annual mean temperature.
2259
2260 DO j = jts , MIN(jde-1,jte)
2261 DO i = its , MIN(ide-1,ite)
2262 IF ( landmask(i,j) .GT. 0.5 ) THEN
2263 DO l = 1 , num_soil_layers
2264 tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
2265 tmn(i,j) * ( zs( l) - zs(1) ) ) / &
2266 ( zs(num_soil_layers) - zs(1) )
2267 END DO
2268 ELSE
2269 IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
2270 DO l = 1 , num_soil_layers
2271 tslb(i,l,j)= sst(i,j)
2272 END DO
2273 ELSE
2274 DO l = 1 , num_soil_layers
2275 tslb(i,l,j)= tsk(i,j)
2276 END DO
2277 END IF
2278 END IF
2279 END DO
2280 END DO
2281
2282 END SUBROUTINE init_soil_1_real
2283
2284 SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
2285 ivgtyp,ZS,DZS,num_soil_layers, &
2286 ids,ide, jds,jde, kds,kde, &
2287 ims,ime, jms,jme, kms,kme, &
2288 its,ite, jts,jte, kts,kte )
2289
2290 IMPLICIT NONE
2291
2292 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2293 ims,ime, jms,jme, kms,kme, &
2294 its,ite, jts,jte, kts,kte
2295
2296 INTEGER, INTENT(IN ) :: num_soil_layers
2297
2298 REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
2299 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
2300 INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
2301
2302 REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
2303
2304 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
2305
2306 ! Lcal variables.
2307
2308 INTEGER :: l,j,i,itf,jtf
2309
2310 itf=MIN(ite,ide-1)
2311 jtf=MIN(jte,jde-1)
2312
2313 IF (num_soil_layers.NE.1)THEN
2314 DO j=jts,jtf
2315 DO l=1,num_soil_layers
2316 DO i=its,itf
2317 tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
2318 ( zs(num_soil_layers)-zs(1) )
2319 ENDDO
2320 ENDDO
2321 ENDDO
2322 ENDIF
2323 DO j=jts,jtf
2324 DO i=its,itf
2325 xland(i,j) = 2
2326 ivgtyp(i,j) = 7
2327 ENDDO
2328 ENDDO
2329
2330 END SUBROUTINE init_soil_1_ideal
2331
2332 SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2333 st_input , sm_input , sw_input , landmask , sst , &
2334 zs , dzs , &
2335 st_levels_input , sm_levels_input , sw_levels_input , &
2336 num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2337 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2338 flag_sst , flag_soilt000 , flag_soilmt000 , &
2339 ids , ide , jds , jde , kds , kde , &
2340 ims , ime , jms , jme , kms , kme , &
2341 its , ite , jts , jte , kts , kte )
2342
2343 IMPLICIT NONE
2344
2345 INTEGER , INTENT(IN) :: num_soil_layers , &
2346 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2347 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2348 ids , ide , jds , jde , kds , kde , &
2349 ims , ime , jms , jme , kms , kme , &
2350 its , ite , jts , jte , kts , kte
2351
2352 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilmt000
2353
2354 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2355 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2356 INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2357
2358 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2359 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2360 REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2361 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2362
2363 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2364 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2365 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2366
2367 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2368
2369 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2370
2371 INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2372 REAL :: temp
2373 LOGICAL :: found_levels
2374 CHARACTER(LEN=132) :: message
2375
2376 ! Are there any soil temp and moisture levels - ya know, they are mandatory.
2377
2378 num = num_st_levels_input * num_sm_levels_input
2379
2380 IF ( num .GE. 1 ) THEN
2381
2382 ! Ordered levels that we have data for.
2383
2384 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2385
2386 ! Sort the levels for temperature.
2387
2388 outert : DO lout = 1 , num_st_levels_input-1
2389 innert : DO lin = lout+1 , num_st_levels_input
2390 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2391 temp = st_levels_input(lout)
2392 st_levels_input(lout) = st_levels_input(lin)
2393 st_levels_input(lin) = NINT(temp)
2394 DO j = jts , MIN(jde-1,jte)
2395 DO i = its , MIN(ide-1,ite)
2396 temp = st_input(i,lout+1,j)
2397 st_input(i,lout+1,j) = st_input(i,lin+1,j)
2398 st_input(i,lin+1,j) = temp
2399 END DO
2400 END DO
2401 END IF
2402 END DO innert
2403 END DO outert
2404
2405 IF ( flag_soilt000 .NE. 1 ) THEN
2406 DO j = jts , MIN(jde-1,jte)
2407 DO i = its , MIN(ide-1,ite)
2408 st_input(i,1,j) = tsk(i,j)
2409 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2410 END DO
2411 END DO
2412 ENDIF
2413
2414
2415 #if ( NMM_CORE == 1 )
2416 !new
2417 ! write(message,*) 'st_input(1) in init_soil_2_real'
2418 ! CALL wrf_message(message)
2419 DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2420 ! write(message,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2421 ! CALL wrf_message(message)
2422 ENDDO
2423
2424 ! write(message,*) 'st_input(2) in init_soil_2_real'
2425 ! CALL wrf_message(message)
2426 DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2427 ! write(message,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2428 ! CALL wrf_message(message)
2429 ENDDO
2430
2431 616 format(20(f4.0,1x))
2432 #endif
2433
2434 ! Sort the levels for moisture.
2435
2436 outerm: DO lout = 1 , num_sm_levels_input-1
2437 innerm : DO lin = lout+1 , num_sm_levels_input
2438 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2439 temp = sm_levels_input(lout)
2440 sm_levels_input(lout) = sm_levels_input(lin)
2441 sm_levels_input(lin) = NINT(temp)
2442 DO j = jts , MIN(jde-1,jte)
2443 DO i = its , MIN(ide-1,ite)
2444 temp = sm_input(i,lout+1,j)
2445 sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2446 sm_input(i,lin+1,j) = temp
2447 END DO
2448 END DO
2449 END IF
2450 END DO innerm
2451 END DO outerm
2452 DO j = jts , MIN(jde-1,jte)
2453 DO i = its , MIN(ide-1,ite)
2454 sm_input(i,1,j) = sm_input(i,2,j)
2455 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2456 END DO
2457 END DO
2458
2459 ! Sort the levels for liquid moisture.
2460
2461 outerw: DO lout = 1 , num_sw_levels_input-1
2462 innerw : DO lin = lout+1 , num_sw_levels_input
2463 IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2464 temp = sw_levels_input(lout)
2465 sw_levels_input(lout) = sw_levels_input(lin)
2466 sw_levels_input(lin) = NINT(temp)
2467 DO j = jts , MIN(jde-1,jte)
2468 DO i = its , MIN(ide-1,ite)
2469 temp = sw_input(i,lout+1,j)
2470 sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2471 sw_input(i,lin+1,j) = temp
2472 END DO
2473 END DO
2474 END IF
2475 END DO innerw
2476 END DO outerw
2477 IF ( num_sw_levels_input .GT. 1 ) THEN
2478 DO j = jts , MIN(jde-1,jte)
2479 DO i = its , MIN(ide-1,ite)
2480 sw_input(i,1,j) = sw_input(i,2,j)
2481 sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2482 END DO
2483 END DO
2484 END IF
2485
2486 found_levels = .TRUE.
2487
2488 ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2489
2490 found_levels = .FALSE.
2491
2492 ELSE
2493 CALL wrf_error_fatal ( &
2494 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
2495 END IF
2496
2497 ! Is it OK to continue?
2498
2499 IF ( found_levels ) THEN
2500
2501 ! Here are the levels that we have from the input for temperature. The input levels plus
2502 ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2503
2504 IF (flag_soilt000 .NE. 1) then
2505 zhave(1) = 0.
2506 DO l = 1 , num_st_levels_input
2507 zhave(l+1) = st_levels_input(l) / 100.
2508 END DO
2509 zhave(num_st_levels_input+2) = 300. / 100.
2510 ELSE
2511 DO l = 1 , num_st_levels_input
2512 zhave(l) = st_levels_input(l) / 100.
2513 END DO
2514 ENDIF
2515
2516
2517 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2518
2519 z_wantt : DO lwant = 1 , num_soil_layers
2520 z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2521 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2522 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2523 DO j = jts , MIN(jde-1,jte)
2524 DO i = its , MIN(ide-1,ite)
2525 tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2526 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2527 ( zhave(lhave+1) - zhave(lhave) )
2528 END DO
2529 END DO
2530 EXIT z_havet
2531 END IF
2532 END DO z_havet
2533 END DO z_wantt
2534
2535 ! Here are the levels that we have from the input for moisture. The input levels plus
2536 ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2537 ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2538 ! as the most deep layer's value.
2539
2540 zhave(1) = 0.
2541 DO l = 1 , num_sm_levels_input
2542 zhave(l+1) = sm_levels_input(l) / 100.
2543 END DO
2544 zhave(num_sm_levels_input+2) = 300. / 100.
2545
2546 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2547
2548 z_wantm : DO lwant = 1 , num_soil_layers
2549 z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2550 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2551 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2552 DO j = jts , MIN(jde-1,jte)
2553 DO i = its , MIN(ide-1,ite)
2554 smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2555 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2556 ( zhave(lhave+1) - zhave(lhave) )
2557 END DO
2558 END DO
2559 EXIT z_havem
2560 END IF
2561 END DO z_havem
2562 END DO z_wantm
2563
2564 ! Any liquid soil moisture to worry about?
2565
2566 IF ( num_sw_levels_input .GT. 1 ) THEN
2567
2568 zhave(1) = 0.
2569 DO l = 1 , num_sw_levels_input
2570 zhave(l+1) = sw_levels_input(l) / 100.
2571 END DO
2572 zhave(num_sw_levels_input+2) = 300. / 100.
2573
2574 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2575
2576 z_wantw : DO lwant = 1 , num_soil_layers
2577 z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
2578 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2579 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2580 DO j = jts , MIN(jde-1,jte)
2581 DO i = its , MIN(ide-1,ite)
2582 sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2583 sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2584 ( zhave(lhave+1) - zhave(lhave) )
2585 END DO
2586 END DO
2587 EXIT z_havew
2588 END IF
2589 END DO z_havew
2590 END DO z_wantw
2591
2592 END IF
2593
2594
2595 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2596 ! used, but they will make a more continuous plot.
2597
2598 IF ( flag_sst .EQ. 1 ) THEN
2599 DO j = jts , MIN(jde-1,jte)
2600 DO i = its , MIN(ide-1,ite)
2601 IF ( landmask(i,j) .LT. 0.5 ) THEN
2602 DO l = 1 , num_soil_layers
2603 tslb(i,l,j)= sst(i,j)
2604 smois(i,l,j)= 1.0
2605 sh2o (i,l,j)= 1.0
2606 END DO
2607 END IF
2608 END DO
2609 END DO
2610 ELSE
2611 DO j = jts , MIN(jde-1,jte)
2612 DO i = its , MIN(ide-1,ite)
2613 IF ( landmask(i,j) .LT. 0.5 ) THEN
2614 DO l = 1 , num_soil_layers
2615 tslb(i,l,j)= tsk(i,j)
2616 smois(i,l,j)= 1.0
2617 sh2o (i,l,j)= 1.0
2618 END DO
2619 END IF
2620 END DO
2621 END DO
2622 END IF
2623
2624 DEALLOCATE (zhave)
2625
2626 END IF
2627
2628 END SUBROUTINE init_soil_2_real
2629
2630 SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
2631 ivgtyp,isltyp,tslb,smois,tmn, &
2632 num_soil_layers, &
2633 ids,ide, jds,jde, kds,kde, &
2634 ims,ime, jms,jme, kms,kme, &
2635 its,ite, jts,jte, kts,kte )
2636
2637 IMPLICIT NONE
2638
2639 INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
2640 ims,ime, jms,jme, kms,kme, &
2641 its,ite, jts,jte, kts,kte
2642
2643 INTEGER, INTENT(IN) ::num_soil_layers
2644
2645 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
2646
2647 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
2648
2649 INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
2650
2651 INTEGER :: icm,jcm,itf,jtf
2652 INTEGER :: i,j,l
2653
2654 itf=min0(ite,ide-1)
2655 jtf=min0(jte,jde-1)
2656
2657 icm = ide/2
2658 jcm = jde/2
2659
2660 DO j=jts,jtf
2661 DO l=1,num_soil_layers
2662 DO i=its,itf
2663
2664 smois(i,1,j)=0.10
2665 smois(i,2,j)=0.10
2666 smois(i,3,j)=0.10
2667 smois(i,4,j)=0.10
2668
2669 tslb(i,1,j)=295.
2670 tslb(i,2,j)=297.
2671 tslb(i,3,j)=293.
2672 tslb(i,4,j)=293.
2673
2674 ENDDO
2675 ENDDO
2676 ENDDO
2677
2678 DO j=jts,jtf
2679 DO i=its,itf
2680 xland(i,j) = 2
2681 tmn(i,j) = 294.
2682 xice(i,j) = 0.
2683 vegfra(i,j) = 0.
2684 snow(i,j) = 0.
2685 canwat(i,j) = 0.
2686 ivgtyp(i,j) = 7
2687 isltyp(i,j) = 8
2688 ENDDO
2689 ENDDO
2690
2691 END SUBROUTINE init_soil_2_ideal
2692
2693 SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
2694 st_input , sm_input , landmask, sst, &
2695 zs , dzs , &
2696 st_levels_input , sm_levels_input , &
2697 num_soil_layers , num_st_levels_input , num_sm_levels_input , &
2698 num_st_levels_alloc , num_sm_levels_alloc , &
2699 flag_sst , flag_soilt000 , flag_soilm000 , &
2700 ids , ide , jds , jde , kds , kde , &
2701 ims , ime , jms , jme , kms , kme , &
2702 its , ite , jts , jte , kts , kte )
2703
2704 IMPLICIT NONE
2705
2706 INTEGER , INTENT(IN) :: num_soil_layers , &
2707 num_st_levels_input , num_sm_levels_input , &
2708 num_st_levels_alloc , num_sm_levels_alloc , &
2709 ids , ide , jds , jde , kds , kde , &
2710 ims , ime , jms , jme , kms , kme , &
2711 its , ite , jts , jte , kts , kte
2712
2713 INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2714
2715 INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2716 INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2717
2718 REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2719 REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2720 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2721
2722 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2723 REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2724 REAL , DIMENSION(num_soil_layers) :: zs , dzs
2725
2726 REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
2727
2728 REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2729
2730 INTEGER :: i , j , l , lout , lin , lwant , lhave
2731 REAL :: temp
2732
2733 ! Allocate the soil layer array used for interpolating.
2734
2735 IF ( ( num_st_levels_input .LE. 0 ) .OR. &
2736 ( num_sm_levels_input .LE. 0 ) ) THEN
2737 PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
2738 CALL wrf_error_fatal ( 'no soil data' )
2739 ELSE
2740 IF ( flag_soilt000 .eq. 1 ) THEN
2741 PRINT '(A)',' Assume RUC LSM 6-level input'
2742 ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
2743 ELSE
2744 PRINT '(A)',' Assume non-RUC LSM input'
2745 ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
2746 END IF
2747 END IF
2748
2749 ! Sort the levels for temperature.
2750
2751 outert : DO lout = 1 , num_st_levels_input-1
2752 innert : DO lin = lout+1 , num_st_levels_input
2753 IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2754 temp = st_levels_input(lout)
2755 st_levels_input(lout) = st_levels_input(lin)
2756 st_levels_input(lin) = NINT(temp)
2757 DO j = jts , MIN(jde-1,jte)
2758 DO i = its , MIN(ide-1,ite)
2759 temp = st_input(i,lout,j)
2760 st_input(i,lout,j) = st_input(i,lin,j)
2761 st_input(i,lin,j) = temp
2762 END DO
2763 END DO
2764 END IF
2765 END DO innert
2766 END DO outert
2767
2768 IF ( flag_soilt000 .NE. 1 ) THEN
2769 DO j = jts , MIN(jde-1,jte)
2770 DO i = its , MIN(ide-1,ite)
2771 st_input(i,1,j) = tsk(i,j)
2772 st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2773 END DO
2774 END DO
2775 END IF
2776
2777 ! Sort the levels for moisture.
2778
2779 outerm: DO lout = 1 , num_sm_levels_input-1
2780 innerm : DO lin = lout+1 , num_sm_levels_input
2781 IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2782 temp = sm_levels_input(lout)
2783 sm_levels_input(lout) = sm_levels_input(lin)
2784 sm_levels_input(lin) = NINT(temp)
2785 DO j = jts , MIN(jde-1,jte)
2786 DO i = its , MIN(ide-1,ite)
2787 temp = sm_input(i,lout,j)
2788 sm_input(i,lout,j) = sm_input(i,lin,j)
2789 sm_input(i,lin,j) = temp
2790 END DO
2791 END DO
2792 END IF
2793 END DO innerm
2794 END DO outerm
2795
2796 IF ( flag_soilm000 .NE. 1 ) THEN
2797 DO j = jts , MIN(jde-1,jte)
2798 DO i = its , MIN(ide-1,ite)
2799 sm_input(i,1,j) = sm_input(i,2,j)
2800 sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2801 END DO
2802 END DO
2803 END IF
2804
2805 ! Here are the levels that we have from the input for temperature.
2806
2807 IF ( flag_soilt000 .EQ. 1 ) THEN
2808 DO l = 1 , num_st_levels_input
2809 zhave(l) = st_levels_input(l) / 100.
2810 END DO
2811
2812 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2813
2814 z_wantt : DO lwant = 1 , num_soil_layers
2815 z_havet : DO lhave = 1 , num_st_levels_input -1
2816 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2817 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2818 DO j = jts , MIN(jde-1,jte)
2819 DO i = its , MIN(ide-1,ite)
2820 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2821 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2822 ( zhave(lhave+1) - zhave(lhave) )
2823 END DO
2824 END DO
2825 EXIT z_havet
2826 END IF
2827 END DO z_havet
2828 END DO z_wantt
2829
2830 ELSE
2831
2832 zhave(1) = 0.
2833 DO l = 1 , num_st_levels_input
2834 zhave(l+1) = st_levels_input(l) / 100.
2835 END DO
2836 zhave(num_st_levels_input+2) = 300. / 100.
2837
2838 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2839
2840 z_wantt_2 : DO lwant = 1 , num_soil_layers
2841 z_havet_2 : DO lhave = 1 , num_st_levels_input +2
2842 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2843 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2844 DO j = jts , MIN(jde-1,jte)
2845 DO i = its , MIN(ide-1,ite)
2846 tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2847 st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2848 ( zhave(lhave+1) - zhave(lhave) )
2849 END DO
2850 END DO
2851 EXIT z_havet_2
2852 END IF
2853 END DO z_havet_2
2854 END DO z_wantt_2
2855
2856 END IF
2857
2858 ! Here are the levels that we have from the input for moisture.
2859
2860 IF ( flag_soilm000 .EQ. 1 ) THEN
2861 DO l = 1 , num_sm_levels_input
2862 zhave(l) = sm_levels_input(l) / 100.
2863 END DO
2864
2865 ! Interpolate between the layers we have (zhave) and those that we want (zs).
2866
2867 z_wantm : DO lwant = 1 , num_soil_layers
2868 z_havem : DO lhave = 1 , num_sm_levels_input -1
2869 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2870 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2871 DO j = jts , MIN(jde-1,jte)
2872 DO i = its , MIN(ide-1,ite)
2873 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2874 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2875 ( zhave(lhave+1) - zhave(lhave) )
2876 END DO
2877 END DO
2878 EXIT z_havem
2879 END IF
2880 END DO z_havem
2881 END DO z_wantm
2882
2883 ELSE
2884
2885 zhave(1) = 0.
2886 DO l = 1 , num_sm_levels_input
2887 zhave(l+1) = sm_levels_input(l) / 100.
2888 END DO
2889 zhave(num_sm_levels_input+2) = 300. / 100.
2890
2891 z_wantm_2 : DO lwant = 1 , num_soil_layers
2892 z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
2893 IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2894 ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2895 DO j = jts , MIN(jde-1,jte)
2896 DO i = its , MIN(ide-1,ite)
2897 smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2898 sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2899 ( zhave(lhave+1) - zhave(lhave) )
2900 END DO
2901 END DO
2902 EXIT z_havem_2
2903 END IF
2904 END DO z_havem_2
2905 END DO z_wantm_2
2906
2907 END IF
2908 ! Over water, put in reasonable values for soil temperature and moisture. These won't be
2909 ! used, but they will make a more continuous plot.
2910
2911 IF ( flag_sst .EQ. 1 ) THEN
2912 DO j = jts , MIN(jde-1,jte)
2913 DO i = its , MIN(ide-1,ite)
2914 IF ( landmask(i,j) .LT. 0.5 ) THEN
2915 DO l = 1 , num_soil_layers
2916 tslb(i,l,j) = sst(i,j)
2917 tsk(i,j) = sst(i,j)
2918 smois(i,l,j)= 1.0
2919 END DO
2920 END IF
2921 END DO
2922 END DO
2923 ELSE
2924 DO j = jts , MIN(jde-1,jte)
2925 DO i = its , MIN(ide-1,ite)
2926 IF ( landmask(i,j) .LT. 0.5 ) THEN
2927 DO l = 1 , num_soil_layers
2928 tslb(i,l,j)= tsk(i,j)
2929 smois(i,l,j)= 1.0
2930 END DO
2931 END IF
2932 END DO
2933 END DO
2934 END IF
2935
2936 DEALLOCATE (zhave)
2937
2938 END SUBROUTINE init_soil_3_real
2939
2940 END MODULE module_soil_pre
2941
2942 #endif
2943