Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_soil_pre.F @ 3709

Last change on this file since 3709 was 11, checked in by aslmd, 15 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 124.3 KB

Line
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)+.5dzs(l-1)+.5dzs(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)+.5dzs(l-1)+.5dzs(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)+.5dzs(l-1)+.5dzs(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)+.5dzs(l-1)+.5dzs(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
2375	! Are there any soil temp and moisture levels - ya know, they are mandatory.
2376
2377	num = num_st_levels_input * num_sm_levels_input
2378
2379	IF ( num .GE. 1 ) THEN
2380
2381	! Ordered levels that we have data for.
2382
2383	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2384
2385	! Sort the levels for temperature.
2386
2387	outert : DO lout = 1 , num_st_levels_input-1
2388	innert : DO lin = lout+1 , num_st_levels_input
2389	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2390	temp = st_levels_input(lout)
2391	st_levels_input(lout) = st_levels_input(lin)
2392	st_levels_input(lin) = NINT(temp)
2393	DO j = jts , MIN(jde-1,jte)
2394	DO i = its , MIN(ide-1,ite)
2395	temp = st_input(i,lout+1,j)
2396	st_input(i,lout+1,j) = st_input(i,lin+1,j)
2397	st_input(i,lin+1,j) = temp
2398	END DO
2399	END DO
2400	END IF
2401	END DO innert
2402	END DO outert
2403	DO j = jts , MIN(jde-1,jte)
2404	DO i = its , MIN(ide-1,ite)
2405	st_input(i,1,j) = tsk(i,j)
2406	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2407	END DO
2408	END DO
2409
2410
2411	#if ( NMM_CORE == 1 )
2412	!new
2413	! write(0,*) 'st_input(1) in init_soil_2_real'
2414	DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2415	! write(0,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2416	ENDDO
2417
2418	! write(0,*) 'st_input(2) in init_soil_2_real'
2419	DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2420	! write(0,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2421	ENDDO
2422
2423	616 format(20(f4.0,1x))
2424	#endif
2425
2426	! Sort the levels for moisture.
2427
2428	outerm: DO lout = 1 , num_sm_levels_input-1
2429	innerm : DO lin = lout+1 , num_sm_levels_input
2430	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2431	temp = sm_levels_input(lout)
2432	sm_levels_input(lout) = sm_levels_input(lin)
2433	sm_levels_input(lin) = NINT(temp)
2434	DO j = jts , MIN(jde-1,jte)
2435	DO i = its , MIN(ide-1,ite)
2436	temp = sm_input(i,lout+1,j)
2437	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2438	sm_input(i,lin+1,j) = temp
2439	END DO
2440	END DO
2441	END IF
2442	END DO innerm
2443	END DO outerm
2444	DO j = jts , MIN(jde-1,jte)
2445	DO i = its , MIN(ide-1,ite)
2446	sm_input(i,1,j) = sm_input(i,2,j)
2447	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2448	END DO
2449	END DO
2450
2451	! Sort the levels for liquid moisture.
2452
2453	outerw: DO lout = 1 , num_sw_levels_input-1
2454	innerw : DO lin = lout+1 , num_sw_levels_input
2455	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2456	temp = sw_levels_input(lout)
2457	sw_levels_input(lout) = sw_levels_input(lin)
2458	sw_levels_input(lin) = NINT(temp)
2459	DO j = jts , MIN(jde-1,jte)
2460	DO i = its , MIN(ide-1,ite)
2461	temp = sw_input(i,lout+1,j)
2462	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2463	sw_input(i,lin+1,j) = temp
2464	END DO
2465	END DO
2466	END IF
2467	END DO innerw
2468	END DO outerw
2469	IF ( num_sw_levels_input .GT. 1 ) THEN
2470	DO j = jts , MIN(jde-1,jte)
2471	DO i = its , MIN(ide-1,ite)
2472	sw_input(i,1,j) = sw_input(i,2,j)
2473	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2474	END DO
2475	END DO
2476	END IF
2477
2478	found_levels = .TRUE.
2479
2480	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2481
2482	found_levels = .FALSE.
2483
2484	ELSE
2485	CALL wrf_error_fatal ( &
2486	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
2487	END IF
2488
2489	! Is it OK to continue?
2490
2491	IF ( found_levels ) THEN
2492
2493	! Here are the levels that we have from the input for temperature. The input levels plus
2494	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2495
2496	zhave(1) = 0.
2497	DO l = 1 , num_st_levels_input
2498	zhave(l+1) = st_levels_input(l) / 100.
2499	END DO
2500	zhave(num_st_levels_input+2) = 300. / 100.
2501
2502	! Interpolate between the layers we have (zhave) and those that we want (zs).
2503
2504	z_wantt : DO lwant = 1 , num_soil_layers
2505	z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2506	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2507	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2508	DO j = jts , MIN(jde-1,jte)
2509	DO i = its , MIN(ide-1,ite)
2510	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2511	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2512	( zhave(lhave+1) - zhave(lhave) )
2513	END DO
2514	END DO
2515	EXIT z_havet
2516	END IF
2517	END DO z_havet
2518	END DO z_wantt
2519
2520	! Here are the levels that we have from the input for moisture. The input levels plus
2521	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2522	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2523	! as the most deep layer's value.
2524
2525	zhave(1) = 0.
2526	DO l = 1 , num_sm_levels_input
2527	zhave(l+1) = sm_levels_input(l) / 100.
2528	END DO
2529	zhave(num_sm_levels_input+2) = 300. / 100.
2530
2531	! Interpolate between the layers we have (zhave) and those that we want (zs).
2532
2533	z_wantm : DO lwant = 1 , num_soil_layers
2534	z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2535	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2536	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2537	DO j = jts , MIN(jde-1,jte)
2538	DO i = its , MIN(ide-1,ite)
2539	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2540	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2541	( zhave(lhave+1) - zhave(lhave) )
2542	END DO
2543	END DO
2544	EXIT z_havem
2545	END IF
2546	END DO z_havem
2547	END DO z_wantm
2548
2549	! Any liquid soil moisture to worry about?
2550
2551	IF ( num_sw_levels_input .GT. 1 ) THEN
2552
2553	zhave(1) = 0.
2554	DO l = 1 , num_sw_levels_input
2555	zhave(l+1) = sw_levels_input(l) / 100.
2556	END DO
2557	zhave(num_sw_levels_input+2) = 300. / 100.
2558
2559	! Interpolate between the layers we have (zhave) and those that we want (zs).
2560
2561	z_wantw : DO lwant = 1 , num_soil_layers
2562	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
2563	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2564	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2565	DO j = jts , MIN(jde-1,jte)
2566	DO i = its , MIN(ide-1,ite)
2567	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2568	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2569	( zhave(lhave+1) - zhave(lhave) )
2570	END DO
2571	END DO
2572	EXIT z_havew
2573	END IF
2574	END DO z_havew
2575	END DO z_wantw
2576
2577	END IF
2578
2579
2580	! Over water, put in reasonable values for soil temperature and moisture. These won't be
2581	! used, but they will make a more continuous plot.
2582
2583	IF ( flag_sst .EQ. 1 ) THEN
2584	DO j = jts , MIN(jde-1,jte)
2585	DO i = its , MIN(ide-1,ite)
2586	IF ( landmask(i,j) .LT. 0.5 ) THEN
2587	DO l = 1 , num_soil_layers
2588	tslb(i,l,j)= sst(i,j)
2589	smois(i,l,j)= 1.0
2590	sh2o (i,l,j)= 1.0
2591	END DO
2592	END IF
2593	END DO
2594	END DO
2595	ELSE
2596	DO j = jts , MIN(jde-1,jte)
2597	DO i = its , MIN(ide-1,ite)
2598	IF ( landmask(i,j) .LT. 0.5 ) THEN
2599	DO l = 1 , num_soil_layers
2600	tslb(i,l,j)= tsk(i,j)
2601	smois(i,l,j)= 1.0
2602	sh2o (i,l,j)= 1.0
2603	END DO
2604	END IF
2605	END DO
2606	END DO
2607	END IF
2608
2609	DEALLOCATE (zhave)
2610
2611	END IF
2612
2613	END SUBROUTINE init_soil_2_real
2614
2615	SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
2616	ivgtyp,isltyp,tslb,smois,tmn, &
2617	num_soil_layers, &
2618	ids,ide, jds,jde, kds,kde, &
2619	ims,ime, jms,jme, kms,kme, &
2620	its,ite, jts,jte, kts,kte )
2621
2622	IMPLICIT NONE
2623
2624	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
2625	ims,ime, jms,jme, kms,kme, &
2626	its,ite, jts,jte, kts,kte
2627
2628	INTEGER, INTENT(IN) ::num_soil_layers
2629
2630	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
2631
2632	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
2633
2634	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
2635
2636	INTEGER :: icm,jcm,itf,jtf
2637	INTEGER :: i,j,l
2638
2639	itf=min0(ite,ide-1)
2640	jtf=min0(jte,jde-1)
2641
2642	icm = ide/2
2643	jcm = jde/2
2644
2645	DO j=jts,jtf
2646	DO l=1,num_soil_layers
2647	DO i=its,itf
2648
2649	smois(i,1,j)=0.10
2650	smois(i,2,j)=0.10
2651	smois(i,3,j)=0.10
2652	smois(i,4,j)=0.10
2653
2654	tslb(i,1,j)=295.
2655	tslb(i,2,j)=297.
2656	tslb(i,3,j)=293.
2657	tslb(i,4,j)=293.
2658
2659	ENDDO
2660	ENDDO
2661	ENDDO
2662
2663	DO j=jts,jtf
2664	DO i=its,itf
2665	xland(i,j) = 2
2666	tmn(i,j) = 294.
2667	xice(i,j) = 0.
2668	vegfra(i,j) = 0.
2669	snow(i,j) = 0.
2670	canwat(i,j) = 0.
2671	ivgtyp(i,j) = 7
2672	isltyp(i,j) = 8
2673	ENDDO
2674	ENDDO
2675
2676	END SUBROUTINE init_soil_2_ideal
2677
2678	SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
2679	st_input , sm_input , landmask, sst, &
2680	zs , dzs , &
2681	st_levels_input , sm_levels_input , &
2682	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
2683	num_st_levels_alloc , num_sm_levels_alloc , &
2684	flag_sst , flag_soilt000 , flag_soilm000 , &
2685	ids , ide , jds , jde , kds , kde , &
2686	ims , ime , jms , jme , kms , kme , &
2687	its , ite , jts , jte , kts , kte )
2688
2689	IMPLICIT NONE
2690
2691	INTEGER , INTENT(IN) :: num_soil_layers , &
2692	num_st_levels_input , num_sm_levels_input , &
2693	num_st_levels_alloc , num_sm_levels_alloc , &
2694	ids , ide , jds , jde , kds , kde , &
2695	ims , ime , jms , jme , kms , kme , &
2696	its , ite , jts , jte , kts , kte
2697
2698	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2699
2700	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2701	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2702
2703	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2704	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2705	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2706
2707	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2708	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2709	REAL , DIMENSION(num_soil_layers) :: zs , dzs
2710
2711	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
2712
2713	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2714
2715	INTEGER :: i , j , l , lout , lin , lwant , lhave
2716	REAL :: temp
2717
2718	! Allocate the soil layer array used for interpolating.
2719
2720	IF ( ( num_st_levels_input .LE. 0 ) .OR. &
2721	( num_sm_levels_input .LE. 0 ) ) THEN
2722	PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
2723	CALL wrf_error_fatal ( 'no soil data' )
2724	ELSE
2725	IF ( flag_soilt000 .eq. 1 ) THEN
2726	PRINT '(A)',' Assume RUC LSM 6-level input'
2727	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
2728	ELSE
2729	PRINT '(A)',' Assume non-RUC LSM input'
2730	ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
2731	END IF
2732	END IF
2733
2734	! Sort the levels for temperature.
2735
2736	outert : DO lout = 1 , num_st_levels_input-1
2737	innert : DO lin = lout+1 , num_st_levels_input
2738	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2739	temp = st_levels_input(lout)
2740	st_levels_input(lout) = st_levels_input(lin)
2741	st_levels_input(lin) = NINT(temp)
2742	DO j = jts , MIN(jde-1,jte)
2743	DO i = its , MIN(ide-1,ite)
2744	temp = st_input(i,lout,j)
2745	st_input(i,lout,j) = st_input(i,lin,j)
2746	st_input(i,lin,j) = temp
2747	END DO
2748	END DO
2749	END IF
2750	END DO innert
2751	END DO outert
2752
2753	IF ( flag_soilt000 .NE. 1 ) THEN
2754	DO j = jts , MIN(jde-1,jte)
2755	DO i = its , MIN(ide-1,ite)
2756	st_input(i,1,j) = tsk(i,j)
2757	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2758	END DO
2759	END DO
2760	END IF
2761
2762	! Sort the levels for moisture.
2763
2764	outerm: DO lout = 1 , num_sm_levels_input-1
2765	innerm : DO lin = lout+1 , num_sm_levels_input
2766	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2767	temp = sm_levels_input(lout)
2768	sm_levels_input(lout) = sm_levels_input(lin)
2769	sm_levels_input(lin) = NINT(temp)
2770	DO j = jts , MIN(jde-1,jte)
2771	DO i = its , MIN(ide-1,ite)
2772	temp = sm_input(i,lout,j)
2773	sm_input(i,lout,j) = sm_input(i,lin,j)
2774	sm_input(i,lin,j) = temp
2775	END DO
2776	END DO
2777	END IF
2778	END DO innerm
2779	END DO outerm
2780
2781	IF ( flag_soilm000 .NE. 1 ) THEN
2782	DO j = jts , MIN(jde-1,jte)
2783	DO i = its , MIN(ide-1,ite)
2784	sm_input(i,1,j) = sm_input(i,2,j)
2785	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2786	END DO
2787	END DO
2788	END IF
2789
2790	! Here are the levels that we have from the input for temperature.
2791
2792	IF ( flag_soilt000 .EQ. 1 ) THEN
2793	DO l = 1 , num_st_levels_input
2794	zhave(l) = st_levels_input(l) / 100.
2795	END DO
2796
2797	! Interpolate between the layers we have (zhave) and those that we want (zs).
2798
2799	z_wantt : DO lwant = 1 , num_soil_layers
2800	z_havet : DO lhave = 1 , num_st_levels_input -1
2801	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2802	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2803	DO j = jts , MIN(jde-1,jte)
2804	DO i = its , MIN(ide-1,ite)
2805	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2806	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2807	( zhave(lhave+1) - zhave(lhave) )
2808	END DO
2809	END DO
2810	EXIT z_havet
2811	END IF
2812	END DO z_havet
2813	END DO z_wantt
2814
2815	ELSE
2816
2817	zhave(1) = 0.
2818	DO l = 1 , num_st_levels_input
2819	zhave(l+1) = st_levels_input(l) / 100.
2820	END DO
2821	zhave(num_st_levels_input+2) = 300. / 100.
2822
2823	! Interpolate between the layers we have (zhave) and those that we want (zs).
2824
2825	z_wantt_2 : DO lwant = 1 , num_soil_layers
2826	z_havet_2 : DO lhave = 1 , num_st_levels_input +2
2827	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2828	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2829	DO j = jts , MIN(jde-1,jte)
2830	DO i = its , MIN(ide-1,ite)
2831	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2832	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2833	( zhave(lhave+1) - zhave(lhave) )
2834	END DO
2835	END DO
2836	EXIT z_havet_2
2837	END IF
2838	END DO z_havet_2
2839	END DO z_wantt_2
2840
2841	END IF
2842
2843	! Here are the levels that we have from the input for moisture.
2844
2845	IF ( flag_soilm000 .EQ. 1 ) THEN
2846	DO l = 1 , num_sm_levels_input
2847	zhave(l) = sm_levels_input(l) / 100.
2848	END DO
2849
2850	! Interpolate between the layers we have (zhave) and those that we want (zs).
2851
2852	z_wantm : DO lwant = 1 , num_soil_layers
2853	z_havem : DO lhave = 1 , num_sm_levels_input -1
2854	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2855	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2856	DO j = jts , MIN(jde-1,jte)
2857	DO i = its , MIN(ide-1,ite)
2858	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2859	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2860	( zhave(lhave+1) - zhave(lhave) )
2861	END DO
2862	END DO
2863	EXIT z_havem
2864	END IF
2865	END DO z_havem
2866	END DO z_wantm
2867
2868	ELSE
2869
2870	zhave(1) = 0.
2871	DO l = 1 , num_sm_levels_input
2872	zhave(l+1) = sm_levels_input(l) / 100.
2873	END DO
2874	zhave(num_sm_levels_input+2) = 300. / 100.
2875
2876	z_wantm_2 : DO lwant = 1 , num_soil_layers
2877	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
2878	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2879	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2880	DO j = jts , MIN(jde-1,jte)
2881	DO i = its , MIN(ide-1,ite)
2882	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2883	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2884	( zhave(lhave+1) - zhave(lhave) )
2885	END DO
2886	END DO
2887	EXIT z_havem_2
2888	END IF
2889	END DO z_havem_2
2890	END DO z_wantm_2
2891
2892	END IF
2893	! Over water, put in reasonable values for soil temperature and moisture. These won't be
2894	! used, but they will make a more continuous plot.
2895
2896	IF ( flag_sst .EQ. 1 ) THEN
2897	DO j = jts , MIN(jde-1,jte)
2898	DO i = its , MIN(ide-1,ite)
2899	IF ( landmask(i,j) .LT. 0.5 ) THEN
2900	DO l = 1 , num_soil_layers
2901	tslb(i,l,j) = sst(i,j)
2902	tsk(i,j) = sst(i,j)
2903	smois(i,l,j)= 1.0
2904	END DO
2905	END IF
2906	END DO
2907	END DO
2908	ELSE
2909	DO j = jts , MIN(jde-1,jte)
2910	DO i = its , MIN(ide-1,ite)
2911	IF ( landmask(i,j) .LT. 0.5 ) THEN
2912	DO l = 1 , num_soil_layers
2913	tslb(i,l,j)= tsk(i,j)
2914	smois(i,l,j)= 1.0
2915	END DO
2916	END IF
2917	END DO
2918	END DO
2919	END IF
2920
2921	DEALLOCATE (zhave)
2922
2923	END SUBROUTINE init_soil_3_real
2924
2925	END MODULE module_soil_pre
2926
2927	#endif
2928

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format