Context Navigation

module_soil_pre.F

Last change on this file was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 144.4 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	ELSE IF ( ( sf_surface_physics .EQ. 7 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
481	CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
482	CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
483	st_input , sm_input , sw_input, landmask , sst , &
484	zs , dzs , &
485	st_levels_input , sm_levels_input , sw_levels_input, &
486	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
487	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
488	flag_sst , flag_soilt000 , flag_soilm000 , &
489	ids , ide , jds , jde , kds , kde , &
490	ims , ime , jms , jme , kms , kme , &
491	its , ite , jts , jte , kts , kte )
492	END IF
493
494	END SUBROUTINE process_soil_real
495
496	SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
497	ivgtyp,isltyp,tslb,smois, &
498	tsk,tmn,zs,dzs, &
499	num_soil_layers, &
500	sf_surface_physics , &
501	ids,ide, jds,jde, kds,kde,&
502	ims,ime, jms,jme, kms,kme,&
503	its,ite, jts,jte, kts,kte )
504
505	IMPLICIT NONE
506
507	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
508	ims,ime, jms,jme, kms,kme, &
509	its,ite, jts,jte, kts,kte
510
511	INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
512
513	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
514
515	REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
516
517	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
518	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
519	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
520
521	! Local variables.
522
523	INTEGER :: itf,jtf
524
525	itf=MIN(ite,ide-1)
526	jtf=MIN(jte,jde-1)
527
528	IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
529	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
530	CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
531	ivgtyp,zs,dzs,num_soil_layers, &
532	ids,ide, jds,jde, kds,kde, &
533	ims,ime, jms,jme, kms,kme, &
534	its,ite, jts,jte, kts,kte )
535	ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
536	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
537	CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
538	ivgtyp,isltyp,tslb,smois,tmn, &
539	num_soil_layers, &
540	ids,ide, jds,jde, kds,kde, &
541	ims,ime, jms,jme, kms,kme, &
542	its,ite, jts,jte, kts,kte )
543	ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
544	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
545
546	END IF
547
548	END SUBROUTINE process_soil_ideal
549
550	SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
551	tsk , ter , toposoil , landmask , flag_toposoil , flag_tavgsfc , &
552	st000010 , st010040 , st040100 , st100200 , st010200 , &
553	flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
554	st000007 , st007028 , st028100 , st100255 , &
555	flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
556	soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
557	flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
558	ids , ide , jds , jde , kds , kde , &
559	ims , ime , jms , jme , kms , kme , &
560	its , ite , jts , jte , kts , kte )
561
562	IMPLICIT NONE
563
564	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
565	ims , ime , jms , jme , kms , kme , &
566	its , ite , jts , jte , kts , kte
567
568	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
569	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
570	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
571	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000007 , st007028 , st028100 , st100255
572	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
573
574	INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
575	INTEGER , INTENT(IN) :: flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255
576	INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
577	INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
578
579	INTEGER :: i , j
580
581	REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
582
583	! Adjust the annual mean temperature as if it is based on from a sea-level elevation
584	! if the value used is from the actual annula mean data set. If the input field to
585	! be used for tmn is one of the first-guess input temp fields, need to do an adjustment
586	! only on the diff in topo from the model terrain and the first-guess terrain.
587
588	SELECT CASE ( sf_surface_physics )
589
590	CASE (LSMSCHEME)
591	DO j = jts , MIN(jde-1,jte)
592	DO i = its , MIN(ide-1,ite)
593	IF (landmask(i,j) .GT. 0.5 ) THEN
594	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
595	END IF
596	END DO
597	END DO
598
599	CASE (RUCLSMSCHEME)
600	DO j = jts , MIN(jde-1,jte)
601	DO i = its , MIN(ide-1,ite)
602	IF (landmask(i,j) .GT. 0.5 ) THEN
603	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
604	END IF
605	END DO
606	END DO
607
608	END SELECT
609
610
611	! Do we have a soil field with which to modify soil temperatures?
612
613	IF ( flag_toposoil .EQ. 1 ) THEN
614
615	DO j = jts , MIN(jde-1,jte)
616	DO i = its , MIN(ide-1,ite)
617
618	! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
619	! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
620	! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
621	! the difference between the soil elevation and the terrain is greater than 3 km means
622	! that the soil data is either all zeros or that the data are inconsistent. Any of these
623	! three conditions is grievous enough to induce a WRF fatality. However, if they are at
624	! a water point, then we can safely ignore them.
625
626	soil_elev_min_val = toposoil(i,j)
627	soil_elev_max_val = toposoil(i,j)
628	soil_elev_min_dif = ter(i,j) - toposoil(i,j)
629	soil_elev_max_dif = ter(i,j) - toposoil(i,j)
630
631	IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
632	CYCLE
633	ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
634	!print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
635	cycle
636	! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
637	ENDIF
638
639	IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
640	CYCLE
641	ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
642	print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
643	cycle
644	CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
645	ENDIF
646
647	IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
648	( landmask(i,j) .LT. 0.5 ) ) THEN
649	CYCLE
650	ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
651	( landmask(i,j) .GT. 0.5 ) ) THEN
652	print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
653	cycle
654	CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
655	ENDIF
656
657	! For each of the fields that we would like to modify, check to see if it came in from the SI.
658	! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
659	! are not at a water point.
660
661	IF (landmask(i,j) .GT. 0.5 ) THEN
662
663	IF ( sf_surface_physics .EQ. SLABSCHEME ) THEN
664	IF ( ( flag_tavgsfc .EQ. 1 ) .OR. &
665	( flag_st010040 .EQ. 1 ) .OR. &
666	( flag_st000010 .EQ. 1 ) .OR. &
667	( flag_soilt020 .EQ. 1 ) .OR. &
668	( flag_st007028 .EQ. 1 ) ) THEN
669	tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
670	ELSE
671	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
672	END IF
673	END IF
674
675	tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
676
677	IF ( flag_st000010 .EQ. 1 ) THEN
678	st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
679	END IF
680	IF ( flag_st010040 .EQ. 1 ) THEN
681	st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
682	END IF
683	IF ( flag_st040100 .EQ. 1 ) THEN
684	st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
685	END IF
686	IF ( flag_st100200 .EQ. 1 ) THEN
687	st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
688	END IF
689	IF ( flag_st010200 .EQ. 1 ) THEN
690	st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
691	END IF
692
693	IF ( flag_st000007 .EQ. 1 ) THEN
694	st000007(i,j) = st000007(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
695	END IF
696	IF ( flag_st007028 .EQ. 1 ) THEN
697	st007028(i,j) = st007028(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
698	END IF
699	IF ( flag_st028100 .EQ. 1 ) THEN
700	st028100(i,j) = st028100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
701	END IF
702	IF ( flag_st100255 .EQ. 1 ) THEN
703	st100255(i,j) = st100255(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
704	END IF
705
706	IF ( flag_soilt000 .EQ. 1 ) THEN
707	soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
708	END IF
709	IF ( flag_soilt005 .EQ. 1 ) THEN
710	soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
711	END IF
712	IF ( flag_soilt020 .EQ. 1 ) THEN
713	soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
714	END IF
715	IF ( flag_soilt040 .EQ. 1 ) THEN
716	soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
717	END IF
718	IF ( flag_soilt160 .EQ. 1 ) THEN
719	soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
720	END IF
721	IF ( flag_soilt300 .EQ. 1 ) THEN
722	soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
723	END IF
724
725	END IF
726	END DO
727	END DO
728
729	END IF
730
731	END SUBROUTINE adjust_soil_temp_new
732
733
734	SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
735
736	IMPLICIT NONE
737
738	INTEGER, INTENT(IN) :: num_soil_layers
739
740	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
741
742	INTEGER :: l
743
744	! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
745	! The distance from the ground level to the midpoint of the layer is given by zs.
746
747	! ------- Ground Level ---------- \|\| \|\| \|\| \|\|
748	! \|\| \|\| \|\| \|\| zs(1) = 0.005 m
749	! -- -- -- -- -- -- -- -- -- \|\| \|\| \|\| \/
750	! \|\| \|\| \|\|
751	! ----------------------------------- \|\| \|\| \|\| \/ dzs(1) = 0.01 m
752	! \|\| \|\| \|\|
753	! \|\| \|\| \|\| zs(2) = 0.02
754	! -- -- -- -- -- -- -- -- -- \|\| \|\| \/
755	! \|\| \|\|
756	! \|\| \|\|
757	! ----------------------------------- \|\| \|\| \/ dzs(2) = 0.02 m
758	! \|\| \|\|
759	! \|\| \|\|
760	! \|\| \|\|
761	! \|\| \|\| zs(3) = 0.05
762	! -- -- -- -- -- -- -- -- -- \|\| \/
763	! \|\|
764	! \|\|
765	! \|\|
766	! \|\|
767	! ----------------------------------- \/ dzs(3) = 0.04 m
768
769	IF ( num_soil_layers .NE. 5 ) THEN
770	PRINT '(A)','Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
771	CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
772	END IF
773
774	dzs(1)=.01
775	zs(1)=.5*dzs(1)
776
777	DO l=2,num_soil_layers
778	dzs(l)=2*dzs(l-1)
779	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
780	ENDDO
781
782	END SUBROUTINE init_soil_depth_1
783
784	SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
785
786	IMPLICIT NONE
787
788	INTEGER, INTENT(IN) :: num_soil_layers
789
790	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
791
792	INTEGER :: l
793
794	dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
795
796	IF ( num_soil_layers .NE. 4 ) THEN
797	PRINT '(A)','Usually, the LSM uses 4 layers. Change this in the namelist.'
798	CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
799	END IF
800
801	zs(1)=.5*dzs(1)
802
803	DO l=2,num_soil_layers
804	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
805	ENDDO
806
807	END SUBROUTINE init_soil_depth_2
808
809	SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
810
811	IMPLICIT NONE
812
813	INTEGER, INTENT(IN) :: num_soil_layers
814
815	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
816
817	INTEGER :: l
818
819	CHARACTER (LEN=132) :: message
820
821	! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
822	! ZS is specified in the namelist: num_soil_layers = 6 or 9.
823	! Other options with number of levels are possible, but
824	! WRF users should change consistently the namelist entry with the
825	! ZS array in this subroutine.
826
827	IF ( num_soil_layers .EQ. 6) THEN
828	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
829	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
830	ELSEIF ( num_soil_layers .EQ. 9) THEN
831	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
832	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
833	ENDIF
834
835	IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
836	write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
837	CALL wrf_error_fatal ( message )
838	END IF
839
840	END SUBROUTINE init_soil_depth_3
841
842	SUBROUTINE init_soil_depth_7 ( zs , dzs , num_soil_layers )
843
844	IMPLICIT NONE
845
846	INTEGER, INTENT(IN) :: num_soil_layers
847
848	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
849
850	INTEGER :: l
851
852	dzs = (/ 0.01 , 0.99 /)
853
854	IF ( num_soil_layers .NE. 2 ) THEN
855	PRINT '(A)','Usually, the PX LSM uses 2 layers. Change this in the namelist.'
856	CALL wrf_error_fatal ( 'PXLSM_uses_2_layers' )
857	END IF
858
859	zs(1) = 0.5 * dzs(1)
860	zs(2) = dzs(1) + 0.5 * dzs(2)
861
862	END SUBROUTINE init_soil_depth_7
863
864	SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
865	num_soil_layers , real_data_init_type , &
866	landmask , sst , flag_sst , &
867	ids , ide , jds , jde , kds , kde , &
868	ims , ime , jms , jme , kms , kme , &
869	its , ite , jts , jte , kts , kte )
870
871	IMPLICIT NONE
872
873	INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
874	ids , ide , jds , jde , kds , kde , &
875	ims , ime , jms , jme , kms , kme , &
876	its , ite , jts , jte , kts , kte
877
878	INTEGER , INTENT(IN) :: flag_sst
879
880	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
881	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
882
883	REAL , DIMENSION(num_soil_layers) :: zs , dzs
884
885	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
886
887	INTEGER :: i , j , l
888
889	! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
890	! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
891	! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
892	! annual mean temperature.
893
894	DO j = jts , MIN(jde-1,jte)
895	DO i = its , MIN(ide-1,ite)
896	IF ( landmask(i,j) .GT. 0.5 ) THEN
897	DO l = 1 , num_soil_layers
898	tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
899	tmn(i,j) * ( zs( l) - zs(1) ) ) / &
900	( zs(num_soil_layers) - zs(1) )
901	END DO
902	ELSE
903	IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
904	DO l = 1 , num_soil_layers
905	tslb(i,l,j)= sst(i,j)
906	END DO
907	ELSE
908	DO l = 1 , num_soil_layers
909	tslb(i,l,j)= tsk(i,j)
910	END DO
911	END IF
912	END IF
913	END DO
914	END DO
915
916	END SUBROUTINE init_soil_1_real
917
918	SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
919	ivgtyp,ZS,DZS,num_soil_layers, &
920	ids,ide, jds,jde, kds,kde, &
921	ims,ime, jms,jme, kms,kme, &
922	its,ite, jts,jte, kts,kte )
923
924	IMPLICIT NONE
925
926	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
927	ims,ime, jms,jme, kms,kme, &
928	its,ite, jts,jte, kts,kte
929
930	INTEGER, INTENT(IN ) :: num_soil_layers
931
932	REAL, DIMENSION( ims:ime , num_soil_layers , jms:jme ), INTENT(OUT) :: tslb
933	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
934	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
935
936	REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
937
938	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
939
940	! Lcal variables.
941
942	INTEGER :: l,j,i,itf,jtf
943
944	itf=MIN(ite,ide-1)
945	jtf=MIN(jte,jde-1)
946
947	IF (num_soil_layers.NE.1)THEN
948	DO j=jts,jtf
949	DO l=1,num_soil_layers
950	DO i=its,itf
951	tslb(i,l,j)=( tsk(i,j)(zs(num_soil_layers)-zs(l)) + tmn(i,j)(zs(l)-zs(1)) ) / &
952	( zs(num_soil_layers)-zs(1) )
953	ENDDO
954	ENDDO
955	ENDDO
956	ENDIF
957
958	! DO j=jts,jtf
959	! DO i=its,itf
960	! xland(i,j) = 2
961	! ivgtyp(i,j) = 7
962	! ENDDO
963	! ENDDO
964
965	END SUBROUTINE init_soil_1_ideal
966
967	SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
968	st_input , sm_input , sw_input , landmask , sst , &
969	zs , dzs , &
970	st_levels_input , sm_levels_input , sw_levels_input , &
971	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
972	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
973	flag_sst , flag_soilt000 , flag_soilm000 , &
974	ids , ide , jds , jde , kds , kde , &
975	ims , ime , jms , jme , kms , kme , &
976	its , ite , jts , jte , kts , kte )
977
978	IMPLICIT NONE
979
980	INTEGER , INTENT(IN) :: num_soil_layers , &
981	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
982	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
983	ids , ide , jds , jde , kds , kde , &
984	ims , ime , jms , jme , kms , kme , &
985	its , ite , jts , jte , kts , kte
986
987	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
988
989	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
990	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
991	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
992
993	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
994	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
995	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
996	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
997
998	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
999	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1000	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1001
1002	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1003
1004	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1005
1006	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1007	REAL :: temp
1008	LOGICAL :: found_levels
1009
1010	CHARACTER (LEN=132) :: message
1011
1012	! Are there any soil temp and moisture levels - ya know, they are mandatory.
1013
1014	num = num_st_levels_input * num_sm_levels_input
1015
1016	IF ( num .GE. 1 ) THEN
1017
1018	! Ordered levels that we have data for.
1019
1020	!tgs add option to initialize from RUCLSM
1021	IF ( flag_soilt000 .eq. 1 ) THEN
1022	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1023	CALL wrf_message ( message )
1024	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
1025	ELSE
1026	write(message, FMT='(A)') ' Assume Noah LSM input'
1027	CALL wrf_message ( message )
1028	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1029	END IF
1030
1031
1032	! Sort the levels for temperature.
1033
1034	outert : DO lout = 1 , num_st_levels_input-1
1035	innert : DO lin = lout+1 , num_st_levels_input
1036	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1037	temp = st_levels_input(lout)
1038	st_levels_input(lout) = st_levels_input(lin)
1039	st_levels_input(lin) = NINT(temp)
1040	DO j = jts , MIN(jde-1,jte)
1041	DO i = its , MIN(ide-1,ite)
1042	temp = st_input(i,lout+1,j)
1043	st_input(i,lout+1,j) = st_input(i,lin+1,j)
1044	st_input(i,lin+1,j) = temp
1045	END DO
1046	END DO
1047	END IF
1048	END DO innert
1049	END DO outert
1050	!tgs add IF
1051	IF ( flag_soilt000 .NE. 1 ) THEN
1052	DO j = jts , MIN(jde-1,jte)
1053	DO i = its , MIN(ide-1,ite)
1054	st_input(i,1,j) = tsk(i,j)
1055	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1056	END DO
1057	END DO
1058	ENDIF
1059
1060	! Sort the levels for moisture.
1061
1062	outerm: DO lout = 1 , num_sm_levels_input-1
1063	innerm : DO lin = lout+1 , num_sm_levels_input
1064	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1065	temp = sm_levels_input(lout)
1066	sm_levels_input(lout) = sm_levels_input(lin)
1067	sm_levels_input(lin) = NINT(temp)
1068	DO j = jts , MIN(jde-1,jte)
1069	DO i = its , MIN(ide-1,ite)
1070	temp = sm_input(i,lout+1,j)
1071	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1072	sm_input(i,lin+1,j) = temp
1073	END DO
1074	END DO
1075	END IF
1076	END DO innerm
1077	END DO outerm
1078	!tgs add IF
1079	IF ( flag_soilm000 .NE. 1 ) THEN
1080	DO j = jts , MIN(jde-1,jte)
1081	DO i = its , MIN(ide-1,ite)
1082	sm_input(i,1,j) = sm_input(i,2,j)
1083	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1084	END DO
1085	END DO
1086	ENDIF
1087
1088	! Sort the levels for liquid moisture.
1089
1090	outerw: DO lout = 1 , num_sw_levels_input-1
1091	innerw : DO lin = lout+1 , num_sw_levels_input
1092	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1093	temp = sw_levels_input(lout)
1094	sw_levels_input(lout) = sw_levels_input(lin)
1095	sw_levels_input(lin) = NINT(temp)
1096	DO j = jts , MIN(jde-1,jte)
1097	DO i = its , MIN(ide-1,ite)
1098	temp = sw_input(i,lout+1,j)
1099	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1100	sw_input(i,lin+1,j) = temp
1101	END DO
1102	END DO
1103	END IF
1104	END DO innerw
1105	END DO outerw
1106	IF ( num_sw_levels_input .GT. 1 ) THEN
1107	DO j = jts , MIN(jde-1,jte)
1108	DO i = its , MIN(ide-1,ite)
1109	sw_input(i,1,j) = sw_input(i,2,j)
1110	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1111	END DO
1112	END DO
1113	END IF
1114
1115	found_levels = .TRUE.
1116
1117	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
1118
1119	found_levels = .FALSE.
1120
1121	ELSE
1122	CALL wrf_error_fatal ( &
1123	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
1124	END IF
1125
1126	! Is it OK to continue?
1127
1128	IF ( found_levels ) THEN
1129
1130	!tgs: Here are the levels that we have from the input for temperature.
1131
1132	IF ( flag_soilt000 .EQ. 1 ) THEN
1133	DO l = 1 , num_st_levels_input
1134	zhave(l) = st_levels_input(l) / 100.
1135	END DO
1136
1137	! Interpolate between the layers we have (zhave) and those that we want (zs).
1138
1139	z_wantt : DO lwant = 1 , num_soil_layers
1140	z_havet : DO lhave = 1 , num_st_levels_input -1
1141	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1142	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1143	DO j = jts , MIN(jde-1,jte)
1144	DO i = its , MIN(ide-1,ite)
1145	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1146	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1147	( zhave(lhave+1) - zhave(lhave) )
1148	END DO
1149	END DO
1150	EXIT z_havet
1151	END IF
1152	END DO z_havet
1153	END DO z_wantt
1154
1155	ELSE
1156
1157	! Here are the levels that we have from the input for temperature. The input levels plus
1158	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1159
1160	zhave(1) = 0.
1161	DO l = 1 , num_st_levels_input
1162	zhave(l+1) = st_levels_input(l) / 100.
1163	END DO
1164	zhave(num_st_levels_input+2) = 300. / 100.
1165
1166	! Interpolate between the layers we have (zhave) and those that we want (zs).
1167
1168	z_wantt_2: DO lwant = 1 , num_soil_layers
1169	z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
1170	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1171	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1172	DO j = jts , MIN(jde-1,jte)
1173	DO i = its , MIN(ide-1,ite)
1174	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1175	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1176	( zhave(lhave+1) - zhave(lhave) )
1177	END DO
1178	END DO
1179	EXIT z_havet_2
1180	END IF
1181	END DO z_havet_2
1182	END DO z_wantt_2
1183
1184	END IF
1185
1186	!tgs:
1187	IF ( flag_soilm000 .EQ. 1 ) THEN
1188	DO l = 1 , num_sm_levels_input
1189	zhave(l) = sm_levels_input(l) / 100.
1190	END DO
1191
1192	! Interpolate between the layers we have (zhave) and those that we want (zs).
1193
1194	z_wantm : DO lwant = 1 , num_soil_layers
1195	z_havem : DO lhave = 1 , num_sm_levels_input -1
1196	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1197	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1198	DO j = jts , MIN(jde-1,jte)
1199	DO i = its , MIN(ide-1,ite)
1200	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1201	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1202	( zhave(lhave+1) - zhave(lhave) )
1203	END DO
1204	END DO
1205	EXIT z_havem
1206	END IF
1207	END DO z_havem
1208	END DO z_wantm
1209
1210	ELSE
1211	! Here are the levels that we have from the input for moisture. The input levels plus
1212	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
1213	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
1214	! as the most deep layer's value.
1215
1216	zhave(1) = 0.
1217	DO l = 1 , num_sm_levels_input
1218	zhave(l+1) = sm_levels_input(l) / 100.
1219	END DO
1220	zhave(num_sm_levels_input+2) = 300. / 100.
1221
1222	! Interpolate between the layers we have (zhave) and those that we want (zs).
1223
1224	z_wantm_2 : DO lwant = 1 , num_soil_layers
1225	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
1226	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1227	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1228	DO j = jts , MIN(jde-1,jte)
1229	DO i = its , MIN(ide-1,ite)
1230	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1231	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1232	( zhave(lhave+1) - zhave(lhave) )
1233	END DO
1234	END DO
1235	EXIT z_havem_2
1236	END IF
1237	END DO z_havem_2
1238	END DO z_wantm_2
1239	ENDIF
1240
1241	! Any liquid soil moisture to worry about?
1242
1243	IF ( num_sw_levels_input .GT. 1 ) THEN
1244
1245	zhave(1) = 0.
1246	DO l = 1 , num_sw_levels_input
1247	zhave(l+1) = sw_levels_input(l) / 100.
1248	END DO
1249	zhave(num_sw_levels_input+2) = 300. / 100.
1250
1251	! Interpolate between the layers we have (zhave) and those that we want (zs).
1252
1253	z_wantw : DO lwant = 1 , num_soil_layers
1254	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1255	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1256	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1257	DO j = jts , MIN(jde-1,jte)
1258	DO i = its , MIN(ide-1,ite)
1259	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1260	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1261	( zhave(lhave+1) - zhave(lhave) )
1262	END DO
1263	END DO
1264	EXIT z_havew
1265	END IF
1266	END DO z_havew
1267	END DO z_wantw
1268
1269	END IF
1270
1271
1272	! Over water, put in reasonable values for soil temperature and moisture. These won't be
1273	! used, but they will make a more continuous plot.
1274
1275	IF ( flag_sst .EQ. 1 ) THEN
1276	DO j = jts , MIN(jde-1,jte)
1277	DO i = its , MIN(ide-1,ite)
1278	IF ( landmask(i,j) .LT. 0.5 ) THEN
1279	DO l = 1 , num_soil_layers
1280	tslb(i,l,j)= sst(i,j)
1281	smois(i,l,j)= 1.0
1282	sh2o (i,l,j)= 1.0
1283	END DO
1284	END IF
1285	END DO
1286	END DO
1287	ELSE
1288	DO j = jts , MIN(jde-1,jte)
1289	DO i = its , MIN(ide-1,ite)
1290	IF ( landmask(i,j) .LT. 0.5 ) THEN
1291	DO l = 1 , num_soil_layers
1292	tslb(i,l,j)= tsk(i,j)
1293	smois(i,l,j)= 1.0
1294	sh2o (i,l,j)= 1.0
1295	END DO
1296	END IF
1297	END DO
1298	END DO
1299	END IF
1300
1301	DEALLOCATE (zhave)
1302
1303	END IF
1304
1305	END SUBROUTINE init_soil_2_real
1306
1307	SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
1308	ivgtyp,isltyp,tslb,smois,tmn, &
1309	num_soil_layers, &
1310	ids,ide, jds,jde, kds,kde, &
1311	ims,ime, jms,jme, kms,kme, &
1312	its,ite, jts,jte, kts,kte )
1313
1314	IMPLICIT NONE
1315
1316	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
1317	ims,ime, jms,jme, kms,kme, &
1318	its,ite, jts,jte, kts,kte
1319
1320	INTEGER, INTENT(IN) ::num_soil_layers
1321
1322	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
1323
1324	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
1325
1326	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1327
1328	INTEGER :: icm,jcm,itf,jtf
1329	INTEGER :: i,j,l
1330
1331	itf=min0(ite,ide-1)
1332	jtf=min0(jte,jde-1)
1333
1334	icm = ide/2
1335	jcm = jde/2
1336
1337	DO j=jts,jtf
1338	DO l=1,num_soil_layers
1339	DO i=its,itf
1340	if (xland(i,j) .lt. 1.5) then
1341	smois(i,1,j)=0.30
1342	smois(i,2,j)=0.30
1343	smois(i,3,j)=0.30
1344	smois(i,4,j)=0.30
1345
1346	tslb(i,1,j)=290.
1347	tslb(i,2,j)=290.
1348	tslb(i,3,j)=290.
1349	tslb(i,4,j)=290.
1350	endif
1351	ENDDO
1352	ENDDO
1353	ENDDO
1354
1355	END SUBROUTINE init_soil_2_ideal
1356
1357	SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
1358	st_input , sm_input , landmask, sst, &
1359	zs , dzs , &
1360	st_levels_input , sm_levels_input , &
1361	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
1362	num_st_levels_alloc , num_sm_levels_alloc , &
1363	flag_sst , flag_soilt000 , flag_soilm000 , &
1364	ids , ide , jds , jde , kds , kde , &
1365	ims , ime , jms , jme , kms , kme , &
1366	its , ite , jts , jte , kts , kte )
1367
1368	IMPLICIT NONE
1369
1370	INTEGER , INTENT(IN) :: num_soil_layers , &
1371	num_st_levels_input , num_sm_levels_input , &
1372	num_st_levels_alloc , num_sm_levels_alloc , &
1373	ids , ide , jds , jde , kds , kde , &
1374	ims , ime , jms , jme , kms , kme , &
1375	its , ite , jts , jte , kts , kte
1376
1377	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
1378
1379	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1380	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1381
1382	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1383	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1384	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1385
1386	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1387	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1388	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1389
1390	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
1391
1392	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1393
1394	INTEGER :: i , j , l , lout , lin , lwant , lhave, k
1395	REAL :: temp
1396
1397	CHARACTER (LEN=132) :: message
1398
1399	! Allocate the soil layer array used for interpolating.
1400
1401	IF ( ( num_st_levels_input .LE. 0 ) .OR. &
1402	( num_sm_levels_input .LE. 0 ) ) THEN
1403	write (message, FMT='(A)')&
1404	'No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
1405	CALL wrf_error_fatal ( message )
1406	ELSE
1407	IF ( flag_soilt000 .eq. 1 ) THEN
1408	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1409	CALL wrf_message ( message )
1410	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
1411	ELSE
1412	write(message, FMT='(A)') ' Assume non-RUC LSM input'
1413	CALL wrf_message ( message )
1414	ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
1415	END IF
1416	END IF
1417
1418	! Sort the levels for temperature.
1419
1420	outert : DO lout = 1 , num_st_levels_input-1
1421	innert : DO lin = lout+1 , num_st_levels_input
1422	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1423	temp = st_levels_input(lout)
1424	st_levels_input(lout) = st_levels_input(lin)
1425	st_levels_input(lin) = NINT(temp)
1426	DO j = jts , MIN(jde-1,jte)
1427	DO i = its , MIN(ide-1,ite)
1428	temp = st_input(i,lout,j)
1429	st_input(i,lout,j) = st_input(i,lin,j)
1430	st_input(i,lin,j) = temp
1431	END DO
1432	END DO
1433	END IF
1434	END DO innert
1435	END DO outert
1436
1437	IF ( flag_soilt000 .NE. 1 ) THEN
1438	DO j = jts , MIN(jde-1,jte)
1439	DO i = its , MIN(ide-1,ite)
1440	st_input(i,1,j) = tsk(i,j)
1441	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1442	END DO
1443	END DO
1444	END IF
1445
1446	! Sort the levels for moisture.
1447
1448	outerm: DO lout = 1 , num_sm_levels_input-1
1449	innerm : DO lin = lout+1 , num_sm_levels_input
1450	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1451	temp = sm_levels_input(lout)
1452	sm_levels_input(lout) = sm_levels_input(lin)
1453	sm_levels_input(lin) = NINT(temp)
1454	DO j = jts , MIN(jde-1,jte)
1455	DO i = its , MIN(ide-1,ite)
1456	temp = sm_input(i,lout,j)
1457	sm_input(i,lout,j) = sm_input(i,lin,j)
1458	sm_input(i,lin,j) = temp
1459	END DO
1460	END DO
1461	END IF
1462	END DO innerm
1463	END DO outerm
1464
1465	IF ( flag_soilm000 .NE. 1 ) THEN
1466	DO j = jts , MIN(jde-1,jte)
1467	DO i = its , MIN(ide-1,ite)
1468	sm_input(i,1,j) = (sm_input(i,2,j)-sm_input(i,3,j))/ &
1469	(st_levels_input(2)-st_levels_input(1))*st_levels_input(1)+ &
1470	sm_input(i,2,j)
1471	! sm_input(i,1,j) = sm_input(i,2,j)
1472	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1473	END DO
1474	END DO
1475	END IF
1476
1477	! Here are the levels that we have from the input for temperature.
1478
1479	IF ( flag_soilt000 .EQ. 1 ) THEN
1480	DO l = 1 , num_st_levels_input
1481	zhave(l) = st_levels_input(l) / 100.
1482	END DO
1483
1484	! Interpolate between the layers we have (zhave) and those that we want (zs).
1485
1486
1487	z_wantt : DO lwant = 1 , num_soil_layers
1488	z_havet : DO lhave = 1 , num_st_levels_input -1
1489	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1490	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1491	DO j = jts , MIN(jde-1,jte)
1492	DO i = its , MIN(ide-1,ite)
1493	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1494	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1495	( zhave(lhave+1) - zhave(lhave) )
1496	END DO
1497	END DO
1498	EXIT z_havet
1499	END IF
1500	END DO z_havet
1501	END DO z_wantt
1502
1503	ELSE
1504
1505	zhave(1) = 0.
1506	DO l = 1 , num_st_levels_input
1507	zhave(l+1) = st_levels_input(l) / 100.
1508	END DO
1509	zhave(num_st_levels_input+2) = 300. / 100.
1510
1511	! Interpolate between the layers we have (zhave) and those that we want (zs).
1512
1513	z_wantt_2 : DO lwant = 1 , num_soil_layers
1514	z_havet_2 : DO lhave = 1 , num_st_levels_input +2
1515	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1516	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1517	DO j = jts , MIN(jde-1,jte)
1518	DO i = its , MIN(ide-1,ite)
1519	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1520	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1521	( zhave(lhave+1) - zhave(lhave) )
1522	END DO
1523	END DO
1524	EXIT z_havet_2
1525	END IF
1526	END DO z_havet_2
1527	END DO z_wantt_2
1528
1529	END IF
1530
1531	! Here are the levels that we have from the input for moisture.
1532
1533	IF ( flag_soilm000 .EQ. 1 ) THEN
1534	DO l = 1 , num_sm_levels_input
1535	zhave(l) = sm_levels_input(l) / 100.
1536	END DO
1537
1538	! Interpolate between the layers we have (zhave) and those that we want (zs).
1539
1540	z_wantm : DO lwant = 1 , num_soil_layers
1541	z_havem : DO lhave = 1 , num_sm_levels_input -1
1542	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1543	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1544	DO j = jts , MIN(jde-1,jte)
1545	DO i = its , MIN(ide-1,ite)
1546	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1547	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1548	( zhave(lhave+1) - zhave(lhave) )
1549	END DO
1550	END DO
1551	EXIT z_havem
1552	END IF
1553	END DO z_havem
1554	END DO z_wantm
1555
1556	ELSE
1557
1558	zhave(1) = 0.
1559	DO l = 1 , num_sm_levels_input
1560	zhave(l+1) = sm_levels_input(l) / 100.
1561	END DO
1562	zhave(num_sm_levels_input+2) = 300. / 100.
1563
1564	z_wantm_2 : DO lwant = 1 , num_soil_layers
1565	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
1566	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1567	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1568	DO j = jts , MIN(jde-1,jte)
1569	DO i = its , MIN(ide-1,ite)
1570	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1571	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1572	( zhave(lhave+1) - zhave(lhave) )
1573	END DO
1574	END DO
1575	EXIT z_havem_2
1576	END IF
1577	END DO z_havem_2
1578	END DO z_wantm_2
1579
1580	END IF
1581	! Over water, put in reasonable values for soil temperature and moisture. These won't be
1582	! used, but they will make a more continuous plot.
1583
1584	IF ( flag_sst .EQ. 1 ) THEN
1585	DO j = jts , MIN(jde-1,jte)
1586	DO i = its , MIN(ide-1,ite)
1587	IF ( landmask(i,j) .LT. 0.5 ) THEN
1588	DO l = 1 , num_soil_layers
1589	tslb(i,l,j) = sst(i,j)
1590	tsk(i,j) = sst(i,j)
1591	smois(i,l,j)= 1.0
1592	END DO
1593	END IF
1594	END DO
1595	END DO
1596	ELSE
1597	DO j = jts , MIN(jde-1,jte)
1598	DO i = its , MIN(ide-1,ite)
1599	IF ( landmask(i,j) .LT. 0.5 ) THEN
1600	DO l = 1 , num_soil_layers
1601	tslb(i,l,j)= tsk(i,j)
1602	smois(i,l,j)= 1.0
1603	END DO
1604	END IF
1605	END DO
1606	END DO
1607	END IF
1608
1609	DEALLOCATE (zhave)
1610
1611	END SUBROUTINE init_soil_3_real
1612	SUBROUTINE init_soil_7_real ( tsk , tmn , smois , sh2o , tslb , &
1613	st_input , sm_input , sw_input , landmask , sst , &
1614	zs , dzs , &
1615	st_levels_input , sm_levels_input , sw_levels_input , &
1616	num_soil_layers , num_st_levels_input , &
1617	num_sm_levels_input , num_sw_levels_input , &
1618	num_st_levels_alloc , num_sm_levels_alloc , &
1619	num_sw_levels_alloc , &
1620	flag_sst , flag_soilt000 , flag_soilmt000 , &
1621	ids , ide , jds , jde , kds , kde , &
1622	ims , ime , jms , jme , kms , kme , &
1623	its , ite , jts , jte , kts , kte )
1624
1625	! for soil temperature and moisture initialization for the PX LSM
1626
1627	IMPLICIT NONE
1628
1629	INTEGER , INTENT(IN) :: num_soil_layers , &
1630	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1631	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1632	ids , ide , jds , jde , kds , kde , &
1633	ims , ime , jms , jme , kms , kme , &
1634	its , ite , jts , jte , kts , kte
1635
1636	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilmt000
1637
1638	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1639	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1640	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1641
1642	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1643	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1644	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1645	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1646
1647	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1648	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1649	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1650
1651	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1652
1653	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1654
1655	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1656	REAL :: temp
1657	LOGICAL :: found_levels
1658
1659	! Are there any soil temp and moisture levels - ya know, they are mandatory.
1660
1661	num = num_st_levels_input * num_sm_levels_input
1662
1663	IF ( num .GE. 1 ) THEN
1664
1665	! Ordered levels that we have data for.
1666
1667	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1668
1669	! Sort the levels for temperature.
1670	outert : DO lout = 1 , num_st_levels_input-1
1671	innert : DO lin = lout+1 , num_st_levels_input
1672	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1673	temp = st_levels_input(lout)
1674	st_levels_input(lout) = st_levels_input(lin)
1675	st_levels_input(lin) = NINT(temp)
1676	DO j = jts , MIN(jde-1,jte)
1677	DO i = its , MIN(ide-1,ite)
1678	temp = st_input(i,lout+1,j)
1679	st_input(i,lout+1,j) = st_input(i,lin+1,j)
1680	st_input(i,lin+1,j) = temp
1681	END DO
1682	END DO
1683	END IF
1684	END DO innert
1685	END DO outert
1686	DO j = jts , MIN(jde-1,jte)
1687	DO i = its , MIN(ide-1,ite)
1688	st_input(i,1,j) = tsk(i,j)
1689	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1690	END DO
1691	END DO
1692
1693	! Sort the levels for moisture.
1694
1695	outerm: DO lout = 1 , num_sm_levels_input-1
1696	innerm : DO lin = lout+1 , num_sm_levels_input
1697	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1698	temp = sm_levels_input(lout)
1699	sm_levels_input(lout) = sm_levels_input(lin)
1700	sm_levels_input(lin) = NINT(temp)
1701	DO j = jts , MIN(jde-1,jte)
1702	DO i = its , MIN(ide-1,ite)
1703	temp = sm_input(i,lout+1,j)
1704	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1705	sm_input(i,lin+1,j) = temp
1706	END DO
1707	END DO
1708	END IF
1709	END DO innerm
1710	END DO outerm
1711	DO j = jts , MIN(jde-1,jte)
1712	DO i = its , MIN(ide-1,ite)
1713	sm_input(i,1,j) = sm_input(i,2,j)
1714	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1715	END DO
1716	END DO
1717
1718	! Sort the levels for liquid moisture.
1719
1720	outerw: DO lout = 1 , num_sw_levels_input-1
1721	innerw : DO lin = lout+1 , num_sw_levels_input
1722	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1723	temp = sw_levels_input(lout)
1724	sw_levels_input(lout) = sw_levels_input(lin)
1725	sw_levels_input(lin) = NINT(temp)
1726	DO j = jts , MIN(jde-1,jte)
1727	DO i = its , MIN(ide-1,ite)
1728	temp = sw_input(i,lout+1,j)
1729	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1730	sw_input(i,lin+1,j) = temp
1731	END DO
1732	END DO
1733	END IF
1734	END DO innerw
1735	END DO outerw
1736	IF ( num_sw_levels_input .GT. 1 ) THEN
1737	DO j = jts , MIN(jde-1,jte)
1738	DO i = its , MIN(ide-1,ite)
1739	sw_input(i,1,j) = sw_input(i,2,j)
1740	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1741	END DO
1742	END DO
1743	END IF
1744
1745	found_levels = .TRUE.
1746
1747	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
1748
1749	found_levels = .FALSE.
1750
1751	ELSE
1752	CALL wrf_error_fatal ( &
1753	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for PX LSM.' )
1754	END IF
1755
1756	! Is it OK to continue?
1757
1758	IF ( found_levels ) THEN
1759
1760	! Here are the levels that we have from the input for temperature. The input levels plus
1761	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1762
1763	zhave(1) = 0.
1764	DO l = 1 , num_st_levels_input
1765	zhave(l+1) = st_levels_input(l) / 100.
1766	END DO
1767	zhave(num_st_levels_input+2) = 300. / 100.
1768
1769	! Interpolate between the layers we have (zhave) and those that we want (zs).
1770
1771	z_wantt : DO lwant = 1 , num_soil_layers
1772	z_havet : DO lhave = 1 , num_st_levels_input +2 -1
1773	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1774	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1775	DO j = jts , MIN(jde-1,jte)
1776	DO i = its , MIN(ide-1,ite)
1777	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1778	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1779	( zhave(lhave+1) - zhave(lhave) )
1780	END DO
1781	END DO
1782	EXIT z_havet
1783	END IF
1784	END DO z_havet
1785	END DO z_wantt
1786
1787	! Here are the levels that we have from the input for moisture. The input levels plus
1788	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
1789	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
1790	! as the most deep layer's value.
1791
1792	zhave(1) = 0.
1793	DO l = 1 , num_sm_levels_input
1794	zhave(l+1) = sm_levels_input(l) / 100.
1795	END DO
1796	zhave(num_sm_levels_input+2) = 300. / 100.
1797
1798	! Interpolate between the layers we have (zhave) and those that we want (zs).
1799
1800	z_wantm : DO lwant = 1 , num_soil_layers
1801	z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
1802	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1803	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1804	DO j = jts , MIN(jde-1,jte)
1805	DO i = its , MIN(ide-1,ite)
1806	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1807	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1808	( zhave(lhave+1) - zhave(lhave) )
1809	END DO
1810	END DO
1811	EXIT z_havem
1812	END IF
1813	END DO z_havem
1814	END DO z_wantm
1815
1816	! Any liquid soil moisture to worry about?
1817
1818	IF ( num_sw_levels_input .GT. 1 ) THEN
1819
1820	zhave(1) = 0.
1821	DO l = 1 , num_sw_levels_input
1822	zhave(l+1) = sw_levels_input(l) / 100.
1823	END DO
1824	zhave(num_sw_levels_input+2) = 300. / 100.
1825
1826	! Interpolate between the layers we have (zhave) and those that we want (zs).
1827
1828	z_wantw : DO lwant = 1 , num_soil_layers
1829	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1830	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1831	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1832	DO j = jts , MIN(jde-1,jte)
1833	DO i = its , MIN(ide-1,ite)
1834	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1835	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1836	( zhave(lhave+1) - zhave(lhave) )
1837	END DO
1838	END DO
1839	EXIT z_havew
1840	END IF
1841	END DO z_havew
1842	END DO z_wantw
1843
1844	END IF
1845
1846
1847	! Over water, put in reasonable values for soil temperature and moisture. These won't be
1848	! used, but they will make a more continuous plot.
1849
1850	DO j = jts , MIN(jde-1,jte)
1851	DO i = its , MIN(ide-1,ite)
1852	tslb(i,1,j)= tsk(i,j)
1853	tslb(i,2,j)= tmn(i,j)
1854	END DO
1855	END DO
1856
1857	IF ( flag_sst .EQ. 1 ) THEN
1858	DO j = jts , MIN(jde-1,jte)
1859	DO i = its , MIN(ide-1,ite)
1860	IF ( landmask(i,j) .LT. 0.5 ) THEN
1861	DO l = 1 , num_soil_layers
1862	tslb(i,l,j)= sst(i,j)
1863	smois(i,l,j)= 1.0
1864	sh2o (i,l,j)= 1.0
1865	END DO
1866	END IF
1867	END DO
1868	END DO
1869	ELSE
1870	DO j = jts , MIN(jde-1,jte)
1871	DO i = its , MIN(ide-1,ite)
1872	IF ( landmask(i,j) .LT. 0.5 ) THEN
1873	DO l = 1 , num_soil_layers
1874	tslb(i,l,j)= tsk(i,j)
1875	smois(i,l,j)= 1.0
1876	sh2o (i,l,j)= 1.0
1877	END DO
1878	END IF
1879	END DO
1880	END DO
1881	END IF
1882
1883	DEALLOCATE (zhave)
1884
1885	END IF
1886
1887	END SUBROUTINE init_soil_7_real
1888
1889
1890	END MODULE module_soil_pre
1891
1892	#else
1893
1894	MODULE module_soil_pre
1895
1896	USE module_date_time
1897	USE module_state_description
1898
1899	CONTAINS
1900
1901	SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1902	xland , landusef , isltyp , soilcat , soilctop , &
1903	soilcbot , tmn , &
1904	seaice_threshold , &
1905	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1906	iswater , isice , &
1907	scheme , &
1908	ids , ide , jds , jde , kds , kde , &
1909	ims , ime , jms , jme , kms , kme , &
1910	its , ite , jts , jte , kts , kte )
1911
1912	IMPLICIT NONE
1913
1914	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1915	ims , ime , jms , jme , kms , kme , &
1916	its , ite , jts , jte , kts , kte , &
1917	iswater , isice
1918	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
1919
1920	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
1921	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
1922	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
1923	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1924	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
1925	vegcat, xland , soilcat , tmn
1926	REAL , INTENT(IN) :: seaice_threshold
1927
1928	INTEGER :: i , j , num_seaice_changes , loop
1929	CHARACTER (LEN=132) :: message
1930
1931	fix_seaice : SELECT CASE ( scheme )
1932
1933	CASE ( SLABSCHEME )
1934	DO j = jts , MIN(jde-1,jte)
1935	DO i = its , MIN(ide-1,ite)
1936	IF ( xice(i,j) .GT. 200.0 ) THEN
1937	xice(i,j) = 0.
1938	num_seaice_changes = num_seaice_changes + 1
1939	END IF
1940	END DO
1941	END DO
1942	IF ( num_seaice_changes .GT. 0 ) THEN
1943	WRITE ( message , FMT='(A,I6)' ) &
1944	'Total pre number of sea ice locations removed (due to FLAG values) = ', &
1945	num_seaice_changes
1946	CALL wrf_debug ( 0 , message )
1947	END IF
1948	num_seaice_changes = 0
1949	DO j = jts , MIN(jde-1,jte)
1950	DO i = its , MIN(ide-1,ite)
1951	IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
1952	( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
1953	xice(i,j) = 1.
1954	num_seaice_changes = num_seaice_changes + 1
1955	tmn(i,j) = 271.4
1956	vegcat(i,j)=isice
1957	lu_index(i,j)=ivgtyp(i,j)
1958	landmask(i,j)=1.
1959	xland(i,j)=1.
1960	DO loop=1,num_veg_cat
1961	landusef(i,loop,j)=0.
1962	END DO
1963	landusef(i,ivgtyp(i,j),j)=1.
1964
1965	isltyp(i,j) = 16
1966	soilcat(i,j)=isltyp(i,j)
1967	DO loop=1,num_soil_top_cat
1968	soilctop(i,loop,j)=0
1969	END DO
1970	DO loop=1,num_soil_bot_cat
1971	soilcbot(i,loop,j)=0
1972	END DO
1973	soilctop(i,isltyp(i,j),j)=1.
1974	soilcbot(i,isltyp(i,j),j)=1.
1975	END IF
1976	END DO
1977	END DO
1978	IF ( num_seaice_changes .GT. 0 ) THEN
1979	WRITE ( message , FMT='(A,I6)' ) &
1980	'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
1981	CALL wrf_debug ( 0 , message )
1982	END IF
1983
1984	CASE ( LSMSCHEME )
1985	CASE ( RUCLSMSCHEME )
1986
1987	END SELECT fix_seaice
1988
1989	END SUBROUTINE adjust_for_seaice_pre
1990
1991	SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1992	xland , landusef , isltyp , soilcat , soilctop , &
1993	soilcbot , tmn , &
1994	tslb , smois , sh2o , &
1995	seaice_threshold , &
1996	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1997	num_soil_layers , &
1998	iswater , isice , &
1999	scheme , &
2000	ids , ide , jds , jde , kds , kde , &
2001	ims , ime , jms , jme , kms , kme , &
2002	its , ite , jts , jte , kts , kte )
2003
2004	IMPLICIT NONE
2005
2006	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2007	ims , ime , jms , jme , kms , kme , &
2008	its , ite , jts , jte , kts , kte , &
2009	iswater , isice
2010	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
2011	INTEGER , INTENT(IN) :: num_soil_layers
2012
2013	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
2014	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
2015	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
2016	REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
2017	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
2018	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
2019	vegcat, xland , soilcat , tmn
2020	REAL , INTENT(IN) :: seaice_threshold
2021	REAL :: total_depth , mid_point_depth
2022
2023	INTEGER :: i , j , num_seaice_changes , loop
2024	CHARACTER (LEN=132) :: message
2025
2026	fix_seaice : SELECT CASE ( scheme )
2027
2028	CASE ( SLABSCHEME )
2029
2030	CASE ( LSMSCHEME )
2031	DO j = jts , MIN(jde-1,jte)
2032	DO i = its , MIN(ide-1,ite)
2033	IF ( xice(i,j) .GT. 200.0 ) THEN
2034	xice(i,j) = 0.
2035	num_seaice_changes = num_seaice_changes + 1
2036	END IF
2037	END DO
2038	END DO
2039	IF ( num_seaice_changes .GT. 0 ) THEN
2040	WRITE ( message , FMT='(A,I6)' ) &
2041	'Total post number of sea ice locations removed (due to FLAG values) = ', &
2042	num_seaice_changes
2043	CALL wrf_debug ( 0 , message )
2044	END IF
2045	num_seaice_changes = 0
2046	DO j = jts , MIN(jde-1,jte)
2047	DO i = its , MIN(ide-1,ite)
2048	IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
2049	( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
2050	xice(i,j) = 1.
2051	num_seaice_changes = num_seaice_changes + 1
2052	tmn(i,j) = 271.16
2053	vegcat(i,j)=isice
2054	lu_index(i,j)=ivgtyp(i,j)
2055	landmask(i,j)=1.
2056	xland(i,j)=1.
2057	DO loop=1,num_veg_cat
2058	landusef(i,loop,j)=0.
2059	END DO
2060	landusef(i,ivgtyp(i,j),j)=1.
2061
2062	isltyp(i,j) = 16
2063	soilcat(i,j)=isltyp(i,j)
2064	DO loop=1,num_soil_top_cat
2065	soilctop(i,loop,j)=0
2066	END DO
2067	DO loop=1,num_soil_bot_cat
2068	soilcbot(i,loop,j)=0
2069	END DO
2070	soilctop(i,isltyp(i,j),j)=1.
2071	soilcbot(i,isltyp(i,j),j)=1.
2072
2073	total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
2074	DO loop = 1,num_soil_layers
2075	mid_point_depth=(total_depth/num_soil_layers)/2. + &
2076	(loop-1)*(total_depth/num_soil_layers)
2077	tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
2078	mid_point_depth*tmn(i,j) ) / total_depth
2079	END DO
2080
2081	DO loop=1,num_soil_layers
2082	smois(i,loop,j) = 1.0
2083	sh2o(i,loop,j) = 0.0
2084	END DO
2085	END IF
2086	END DO
2087	END DO
2088	IF ( num_seaice_changes .GT. 0 ) THEN
2089	WRITE ( message , FMT='(A,I6)' ) &
2090	'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
2091	CALL wrf_debug ( 0 , message )
2092	END IF
2093
2094	CASE ( RUCLSMSCHEME )
2095
2096	END SELECT fix_seaice
2097
2098	END SUBROUTINE adjust_for_seaice_post
2099
2100	SUBROUTINE process_percent_cat_new ( landmask , &
2101	landuse_frac , soil_top_cat , soil_bot_cat , &
2102	isltyp , ivgtyp , &
2103	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
2104	ids , ide , jds , jde , kds , kde , &
2105	ims , ime , jms , jme , kms , kme , &
2106	its , ite , jts , jte , kts , kte , &
2107	iswater )
2108
2109	IMPLICIT NONE
2110
2111	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2112	ims , ime , jms , jme , kms , kme , &
2113	its , ite , jts , jte , kts , kte , &
2114	iswater
2115	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
2116	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
2117	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
2118	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
2119	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
2120	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
2121
2122	INTEGER :: i , j , l , ll, dominant_index
2123	REAL :: dominant_value
2124
2125	#ifdef WRF_CHEM
2126	! REAL :: lwthresh = .99
2127	REAL :: lwthresh = .50
2128	#else
2129	REAL :: lwthresh = .50
2130	#endif
2131
2132	INTEGER , PARAMETER :: iswater_soil = 14
2133	INTEGER :: iforce
2134	CHARACTER (LEN=1024) :: message
2135
2136	! Sanity check on the 50/50 points
2137
2138	DO j = jts , MIN(jde-1,jte)
2139	DO i = its , MIN(ide-1,ite)
2140	dominant_value = landuse_frac(i,iswater,j)
2141	IF ( dominant_value .EQ. lwthresh ) THEN
2142	DO l = 1 , num_veg_cat
2143	IF ( l .EQ. iswater ) CYCLE
2144	IF ( landuse_frac(i,l,j) .EQ. lwthresh ) THEN
2145	write(message, FMT='(I4,I4,A,I4,A,F12.4)') i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
2146	call wrf_message ( message )
2147	landuse_frac(i,l,j) = lwthresh - .01
2148	landuse_frac(i,iswater,j) = 1. - (lwthresh + 0.01)
2149	END IF
2150	END DO
2151	END IF
2152	END DO
2153	END DO
2154
2155	! Compute the dominant VEGETATION INDEX.
2156
2157	DO j = jts , MIN(jde-1,jte)
2158	DO i = its , MIN(ide-1,ite)
2159	dominant_value = landuse_frac(i,1,j)
2160	dominant_index = 1
2161	DO l = 2 , num_veg_cat
2162	IF ( l .EQ. iswater ) THEN
2163	! wait a bit
2164	ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
2165	dominant_value = landuse_frac(i,l,j)
2166	dominant_index = l
2167	END IF
2168	END DO
2169	IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
2170	dominant_value = landuse_frac(i,iswater,j)
2171	dominant_index = iswater
2172	#if 0
2173	si needs to beef up consistency checks before we can use this part
2174	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
2175	( dominant_value .EQ. lwthresh) ) THEN
2176	! no op
2177	#else
2178	ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.LT.iswater))THEN
2179	write(message,*)'temporary SI landmask fix'
2180	call wrf_message(trim(message))
2181	! no op
2182	ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.GT.iswater))THEN
2183	write(message,*)'temporary SI landmask fix'
2184	call wrf_message(trim(message))
2185	dominant_value = landuse_frac(i,iswater,j)
2186	dominant_index = iswater
2187	#endif
2188	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
2189	( dominant_value .LT. lwthresh ) ) THEN
2190	dominant_value = landuse_frac(i,iswater,j)
2191	dominant_index = iswater
2192	END IF
2193	IF ( dominant_index .EQ. iswater ) THEN
2194	if(landmask(i,j).gt.lwthresh) then
2195	write(message,*) 'changing to water at point ',i,j
2196	call wrf_message(trim(message))
2197	write(message,) nint(landuse_frac(i,:,j)100)
2198	call wrf_message(trim(message))
2199	endif
2200	landmask(i,j) = 0
2201	ELSE IF ( dominant_index .NE. iswater ) THEN
2202	if(landmask(i,j).lt.lwthresh) then
2203	write(message,*) 'changing to land at point ',i,j
2204	call wrf_message(trim(message))
2205	write(message,) nint(landuse_frac(i,:,j)100)
2206	call wrf_message(trim(message))
2207	endif
2208	landmask(i,j) = 1
2209	END IF
2210	ivgtyp(i,j) = dominant_index
2211	END DO
2212	END DO
2213
2214	! Compute the dominant SOIL TEXTURE INDEX, TOP.
2215
2216	iforce = 0
2217	DO i = its , MIN(ide-1,ite)
2218	DO j = jts , MIN(jde-1,jte)
2219	dominant_value = soil_top_cat(i,1,j)
2220	dominant_index = 1
2221	IF ( landmask(i,j) .GT. lwthresh ) THEN
2222	DO l = 2 , num_soil_top_cat
2223	IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
2224	dominant_value = soil_top_cat(i,l,j)
2225	dominant_index = l
2226	END IF
2227	END DO
2228	IF ( dominant_value .LT. 0.01 ) THEN
2229	iforce = iforce + 1
2230	WRITE ( message , FMT = '(A,I4,I4)' ) &
2231	'based on landuse, changing soil to land at point ',i,j
2232	CALL wrf_debug(1,message)
2233	WRITE ( message , FMT = '(16(i3,1x))' ) &
2234	1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
2235	CALL wrf_debug(1,message)
2236	WRITE ( message , FMT = '(16(i3,1x))' ) &
2237	nint(soil_top_cat(i,:,j)*100)
2238	CALL wrf_debug(1,message)
2239	dominant_index = 8
2240	END IF
2241	ELSE
2242	dominant_index = iswater_soil
2243	END IF
2244	isltyp(i,j) = dominant_index
2245	END DO
2246	END DO
2247
2248	if(iforce.ne.0)then
2249	WRITE(message,FMT='(A,I4,A,I6)' ) &
2250	'forcing artificial silty clay loam at ',iforce,' points, out of ',&
2251	(MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
2252	CALL wrf_debug(0,message)
2253	endif
2254
2255	END SUBROUTINE process_percent_cat_new
2256
2257	SUBROUTINE process_soil_real ( tsk , tmn , &
2258	landmask , sst , &
2259	st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
2260	zs , dzs , tslb , smois , sh2o , &
2261	flag_sst , flag_soilt000, flag_soilm000, &
2262	ids , ide , jds , jde , kds , kde , &
2263	ims , ime , jms , jme , kms , kme , &
2264	its , ite , jts , jte , kts , kte , &
2265	sf_surface_physics , num_soil_layers , real_data_init_type , &
2266	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2267	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
2268
2269	IMPLICIT NONE
2270
2271	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2272	ims , ime , jms , jme , kms , kme , &
2273	its , ite , jts , jte , kts , kte , &
2274	sf_surface_physics , num_soil_layers , real_data_init_type , &
2275	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2276	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
2277
2278	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2279
2280	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2281
2282	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2283	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2284	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2285	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2286	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2287	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2288
2289	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2290	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2291	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2292	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2293
2294	INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
2295	REAL :: dominant_value
2296
2297	! Initialize the soil depth, and the soil temperature and moisture.
2298
2299	IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2300	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2301	CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
2302	landmask , sst , flag_sst , &
2303	ids , ide , jds , jde , kds , kde , &
2304	ims , ime , jms , jme , kms , kme , &
2305	its , ite , jts , jte , kts , kte )
2306
2307	ELSE IF ( ( sf_surface_physics .EQ. 2 .or. sf_surface_physics .EQ. 99 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2308	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2309	CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2310	st_input , sm_input , sw_input , landmask , sst , &
2311	zs , dzs , &
2312	st_levels_input , sm_levels_input , sw_levels_input , &
2313	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2314	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2315	flag_sst , flag_soilt000 , flag_soilm000 , &
2316	ids , ide , jds , jde , kds , kde , &
2317	ims , ime , jms , jme , kms , kme , &
2318	its , ite , jts , jte , kts , kte )
2319	! CALL init_soil_old ( tsk , tmn , &
2320	! smois , tslb , zs , dzs , num_soil_layers , &
2321	! st000010_input , st010040_input , st040100_input , st100200_input , &
2322	! st010200_input , &
2323	! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
2324	! sm010200_input , &
2325	! landmask_input , sst_input , &
2326	! ids , ide , jds , jde , kds , kde , &
2327	! ims , ime , jms , jme , kms , kme , &
2328	! its , ite , jts , jte , kts , kte )
2329	ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2330	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2331	CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
2332	st_input , sm_input , landmask , sst , &
2333	zs , dzs , &
2334	st_levels_input , sm_levels_input , &
2335	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
2336	num_st_levels_alloc , num_sm_levels_alloc , &
2337	flag_sst , flag_soilt000 , flag_soilm000 , &
2338	ids , ide , jds , jde , kds , kde , &
2339	ims , ime , jms , jme , kms , kme , &
2340	its , ite , jts , jte , kts , kte )
2341	END IF
2342
2343	END SUBROUTINE process_soil_real
2344
2345	SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
2346	ivgtyp,isltyp,tslb,smois, &
2347	tsk,tmn,zs,dzs, &
2348	num_soil_layers, &
2349	sf_surface_physics , &
2350	ids,ide, jds,jde, kds,kde,&
2351	ims,ime, jms,jme, kms,kme,&
2352	its,ite, jts,jte, kts,kte )
2353
2354	IMPLICIT NONE
2355
2356	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
2357	ims,ime, jms,jme, kms,kme, &
2358	its,ite, jts,jte, kts,kte
2359
2360	INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
2361
2362	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
2363
2364	REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
2365
2366	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
2367	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
2368	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
2369
2370	! Local variables.
2371
2372	INTEGER :: itf,jtf
2373
2374	itf=MIN(ite,ide-1)
2375	jtf=MIN(jte,jde-1)
2376
2377	IF ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2378	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2379	CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
2380	ivgtyp,zs,dzs,num_soil_layers, &
2381	ids,ide, jds,jde, kds,kde, &
2382	ims,ime, jms,jme, kms,kme, &
2383	its,ite, jts,jte, kts,kte )
2384	ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2385	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2386	CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
2387	ivgtyp,isltyp,tslb,smois,tmn, &
2388	num_soil_layers, &
2389	ids,ide, jds,jde, kds,kde, &
2390	ims,ime, jms,jme, kms,kme, &
2391	its,ite, jts,jte, kts,kte )
2392	ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2393	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2394
2395	END IF
2396
2397	END SUBROUTINE process_soil_ideal
2398
2399	SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
2400	tsk , ter , toposoil , landmask , flag_toposoil , &
2401	st000010 , st010040 , st040100 , st100200 , st010200 , &
2402	flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
2403	soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
2404	flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
2405	ids , ide , jds , jde , kds , kde , &
2406	ims , ime , jms , jme , kms , kme , &
2407	its , ite , jts , jte , kts , kte )
2408
2409	IMPLICIT NONE
2410
2411	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2412	ims , ime , jms , jme , kms , kme , &
2413	its , ite , jts , jte , kts , kte
2414
2415	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
2416	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
2417	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
2418	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
2419
2420	INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
2421	INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
2422	INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
2423
2424	INTEGER :: i , j
2425
2426	REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
2427
2428	! Do we have a soil field with which to modify soil temperatures?
2429
2430	IF ( flag_toposoil .EQ. 1 ) THEN
2431
2432	DO j = jts , MIN(jde-1,jte)
2433	DO i = its , MIN(ide-1,ite)
2434
2435	! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
2436	! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
2437	! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
2438	! the difference between the soil elevation and the terrain is greater than 3 km means
2439	! that the soil data is either all zeros or that the data are inconsistent. Any of these
2440	! three conditions is grievous enough to induce a WRF fatality. However, if they are at
2441	! a water point, then we can safely ignore them.
2442
2443	soil_elev_min_val = toposoil(i,j)
2444	soil_elev_max_val = toposoil(i,j)
2445	soil_elev_min_dif = ter(i,j) - toposoil(i,j)
2446	soil_elev_max_dif = ter(i,j) - toposoil(i,j)
2447
2448	IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2449	CYCLE
2450	ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2451	!print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
2452	cycle
2453	! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
2454	ENDIF
2455
2456	IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2457	CYCLE
2458	ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2459	print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
2460	cycle
2461	CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
2462	ENDIF
2463
2464	IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2465	( landmask(i,j) .LT. 0.5 ) ) THEN
2466	CYCLE
2467	ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2468	( landmask(i,j) .GT. 0.5 ) ) THEN
2469	print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
2470	cycle
2471	CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
2472	ENDIF
2473
2474	! For each of the fields that we would like to modify, check to see if it came in from the SI.
2475	! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
2476	! are not at a water point.
2477
2478	IF (landmask(i,j) .GT. 0.5 ) THEN
2479
2480	IF ( sf_surface_physics .EQ. 1 .OR. sf_surface_physics .EQ. 7) THEN
2481	tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2482	END IF
2483
2484	tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2485
2486	IF ( flag_st000010 .EQ. 1 ) THEN
2487	st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2488	END IF
2489	IF ( flag_st010040 .EQ. 1 ) THEN
2490	st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2491	END IF
2492	IF ( flag_st040100 .EQ. 1 ) THEN
2493	st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2494	END IF
2495	IF ( flag_st100200 .EQ. 1 ) THEN
2496	st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2497	END IF
2498	IF ( flag_st010200 .EQ. 1 ) THEN
2499	st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2500	END IF
2501
2502	IF ( flag_soilt000 .EQ. 1 ) THEN
2503	soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2504	END IF
2505	IF ( flag_soilt005 .EQ. 1 ) THEN
2506	soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2507	END IF
2508	IF ( flag_soilt020 .EQ. 1 ) THEN
2509	soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2510	END IF
2511	IF ( flag_soilt040 .EQ. 1 ) THEN
2512	soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2513	END IF
2514	IF ( flag_soilt160 .EQ. 1 ) THEN
2515	soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2516	END IF
2517	IF ( flag_soilt300 .EQ. 1 ) THEN
2518	soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2519	END IF
2520
2521	END IF
2522	END DO
2523	END DO
2524
2525	END IF
2526
2527	END SUBROUTINE adjust_soil_temp_new
2528
2529
2530	SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
2531
2532	IMPLICIT NONE
2533
2534	INTEGER, INTENT(IN) :: num_soil_layers
2535
2536	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2537
2538	INTEGER :: l
2539
2540	CHARACTER (LEN=132) :: message
2541
2542	! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
2543	! The distance from the ground level to the midpoint of the layer is given by zs.
2544
2545	! ------- Ground Level ---------- \|\| \|\| \|\| \|\|
2546	! \|\| \|\| \|\| \|\| zs(1) = 0.005 m
2547	! -- -- -- -- -- -- -- -- -- \|\| \|\| \|\| \/
2548	! \|\| \|\| \|\|
2549	! ----------------------------------- \|\| \|\| \|\| \/ dzs(1) = 0.01 m
2550	! \|\| \|\| \|\|
2551	! \|\| \|\| \|\| zs(2) = 0.02
2552	! -- -- -- -- -- -- -- -- -- \|\| \|\| \/
2553	! \|\| \|\|
2554	! \|\| \|\|
2555	! ----------------------------------- \|\| \|\| \/ dzs(2) = 0.02 m
2556	! \|\| \|\|
2557	! \|\| \|\|
2558	! \|\| \|\|
2559	! \|\| \|\| zs(3) = 0.05
2560	! -- -- -- -- -- -- -- -- -- \|\| \/
2561	! \|\|
2562	! \|\|
2563	! \|\|
2564	! \|\|
2565	! ----------------------------------- \/ dzs(3) = 0.04 m
2566
2567	IF ( num_soil_layers .NE. 5 ) THEN
2568	write(message,FMT= '(A)') 'Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
2569	CALL wrf_error_fatal ( message )
2570	END IF
2571
2572	dzs(1)=.01
2573	zs(1)=.5*dzs(1)
2574
2575	DO l=2,num_soil_layers
2576	dzs(l)=2*dzs(l-1)
2577	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
2578	ENDDO
2579
2580	END SUBROUTINE init_soil_depth_1
2581
2582	SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
2583
2584	IMPLICIT NONE
2585
2586	INTEGER, INTENT(IN) :: num_soil_layers
2587
2588	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2589
2590	INTEGER :: l
2591
2592	CHARACTER (LEN=132) :: message
2593
2594	dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
2595
2596	IF ( num_soil_layers .NE. 4 ) THEN
2597	write(message,FMT='(A)') 'Usually, the LSM uses 4 layers. Change this in the namelist.'
2598	CALL wrf_error_fatal ( message )
2599	END IF
2600
2601	zs(1)=.5*dzs(1)
2602
2603	DO l=2,num_soil_layers
2604	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
2605	ENDDO
2606
2607	END SUBROUTINE init_soil_depth_2
2608
2609	SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
2610
2611	IMPLICIT NONE
2612
2613	INTEGER, INTENT(IN) :: num_soil_layers
2614
2615	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2616
2617	INTEGER :: l
2618
2619	CHARACTER (LEN=132) :: message
2620
2621	! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
2622	! ZS is specified in the namelist: num_soil_layers = 6 or 9.
2623	! Other options with number of levels are possible, but
2624	! WRF users should change consistently the namelist entry with the
2625	! ZS array in this subroutine.
2626
2627	IF ( num_soil_layers .EQ. 6) THEN
2628	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
2629	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2630	ELSEIF ( num_soil_layers .EQ. 9) THEN
2631	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
2632	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2633	ENDIF
2634
2635	IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
2636	WRITE(message,FMT= '(A)')'Usually, the RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
2637	CALL wrf_error_fatal ( message )
2638	END IF
2639
2640	END SUBROUTINE init_soil_depth_3
2641
2642	SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
2643	num_soil_layers , real_data_init_type , &
2644	landmask , sst , flag_sst , &
2645	ids , ide , jds , jde , kds , kde , &
2646	ims , ime , jms , jme , kms , kme , &
2647	its , ite , jts , jte , kts , kte )
2648
2649	IMPLICIT NONE
2650
2651	INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
2652	ids , ide , jds , jde , kds , kde , &
2653	ims , ime , jms , jme , kms , kme , &
2654	its , ite , jts , jte , kts , kte
2655
2656	INTEGER , INTENT(IN) :: flag_sst
2657
2658	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2659	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
2660
2661	REAL , DIMENSION(num_soil_layers) :: zs , dzs
2662
2663	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
2664
2665	INTEGER :: i , j , l
2666
2667	! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
2668	! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
2669	! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
2670	! annual mean temperature.
2671
2672	DO j = jts , MIN(jde-1,jte)
2673	DO i = its , MIN(ide-1,ite)
2674	IF ( landmask(i,j) .GT. 0.5 ) THEN
2675	DO l = 1 , num_soil_layers
2676	tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
2677	tmn(i,j) * ( zs( l) - zs(1) ) ) / &
2678	( zs(num_soil_layers) - zs(1) )
2679	END DO
2680	ELSE
2681	IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
2682	DO l = 1 , num_soil_layers
2683	tslb(i,l,j)= sst(i,j)
2684	END DO
2685	ELSE
2686	DO l = 1 , num_soil_layers
2687	tslb(i,l,j)= tsk(i,j)
2688	END DO
2689	END IF
2690	END IF
2691	END DO
2692	END DO
2693
2694	END SUBROUTINE init_soil_1_real
2695
2696	SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
2697	ivgtyp,ZS,DZS,num_soil_layers, &
2698	ids,ide, jds,jde, kds,kde, &
2699	ims,ime, jms,jme, kms,kme, &
2700	its,ite, jts,jte, kts,kte )
2701
2702	IMPLICIT NONE
2703
2704	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2705	ims,ime, jms,jme, kms,kme, &
2706	its,ite, jts,jte, kts,kte
2707
2708	INTEGER, INTENT(IN ) :: num_soil_layers
2709
2710	REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
2711	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
2712	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
2713
2714	REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
2715
2716	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
2717
2718	! Lcal variables.
2719
2720	INTEGER :: l,j,i,itf,jtf
2721
2722	itf=MIN(ite,ide-1)
2723	jtf=MIN(jte,jde-1)
2724
2725	IF (num_soil_layers.NE.1)THEN
2726	DO j=jts,jtf
2727	DO l=1,num_soil_layers
2728	DO i=its,itf
2729	tslb(i,l,j)=( tsk(i,j)(zs(num_soil_layers)-zs(l)) + tmn(i,j)(zs(l)-zs(1)) ) / &
2730	( zs(num_soil_layers)-zs(1) )
2731	ENDDO
2732	ENDDO
2733	ENDDO
2734	ENDIF
2735	DO j=jts,jtf
2736	DO i=its,itf
2737	xland(i,j) = 2
2738	ivgtyp(i,j) = 7
2739	ENDDO
2740	ENDDO
2741
2742	END SUBROUTINE init_soil_1_ideal
2743
2744	SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2745	st_input , sm_input , sw_input , landmask , sst , &
2746	zs , dzs , &
2747	st_levels_input , sm_levels_input , sw_levels_input , &
2748	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2749	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2750	flag_sst , flag_soilt000 , flag_soilm000 , &
2751	ids , ide , jds , jde , kds , kde , &
2752	ims , ime , jms , jme , kms , kme , &
2753	its , ite , jts , jte , kts , kte )
2754
2755	IMPLICIT NONE
2756
2757	INTEGER , INTENT(IN) :: num_soil_layers , &
2758	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2759	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2760	ids , ide , jds , jde , kds , kde , &
2761	ims , ime , jms , jme , kms , kme , &
2762	its , ite , jts , jte , kts , kte
2763
2764	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2765
2766	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2767	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2768	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2769
2770	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2771	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2772	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2773	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2774
2775	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2776	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2777	REAL , DIMENSION(num_soil_layers) :: zs , dzs
2778
2779	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2780
2781	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2782
2783	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2784	REAL :: temp
2785	LOGICAL :: found_levels
2786
2787	CHARACTER (LEN=132) :: message
2788
2789	! Are there any soil temp and moisture levels - ya know, they are mandatory.
2790
2791	num = num_st_levels_input * num_sm_levels_input
2792
2793	IF ( num .GE. 1 ) THEN
2794
2795	! Ordered levels that we have data for.
2796	IF ( flag_soilt000 .eq. 1 ) THEN
2797	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
2798	CALL wrf_message ( message )
2799	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
2800	ELSE
2801	write(message, FMT='(A)') ' Assume Noah LSM input'
2802	CALL wrf_message ( message )
2803	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2804	END IF
2805
2806
2807	! Sort the levels for temperature.
2808	!print*,'num_st_levels_input, st_levels_input', num_st_levels_input, st_levels_input
2809	!print*,'num_sm_levels_input,num_sw_levels_input',num_sm_levels_input,num_sw_levels_input
2810
2811	outert : DO lout = 1 , num_st_levels_input-1
2812	innert : DO lin = lout+1 , num_st_levels_input
2813	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2814	temp = st_levels_input(lout)
2815	st_levels_input(lout) = st_levels_input(lin)
2816	st_levels_input(lin) = NINT(temp)
2817	DO j = jts , MIN(jde-1,jte)
2818	DO i = its , MIN(ide-1,ite)
2819	temp = st_input(i,lout+1,j)
2820	st_input(i,lout+1,j) = st_input(i,lin+1,j)
2821	st_input(i,lin+1,j) = temp
2822	END DO
2823	END DO
2824	END IF
2825	END DO innert
2826	END DO outert
2827	!tgs add IF
2828	IF ( flag_soilt000 .NE. 1 ) THEN
2829	DO j = jts , MIN(jde-1,jte)
2830	DO i = its , MIN(ide-1,ite)
2831	st_input(i,1,j) = tsk(i,j)
2832	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2833	END DO
2834	END DO
2835	ENDIF
2836
2837
2838	#if ( NMM_CORE == 1 )
2839	!new
2840	! write(0,*) 'st_input(1) in init_soil_2_real'
2841	! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2842	! write(0,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2843	! ENDDO
2844
2845	! write(0,*) 'st_input(2) in init_soil_2_real'
2846	! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2847	! write(0,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2848	! ENDDO
2849
2850	616 format(20(f4.0,1x))
2851	#endif
2852
2853	! Sort the levels for moisture.
2854
2855	outerm: DO lout = 1 , num_sm_levels_input-1
2856	innerm : DO lin = lout+1 , num_sm_levels_input
2857	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2858	temp = sm_levels_input(lout)
2859	sm_levels_input(lout) = sm_levels_input(lin)
2860	sm_levels_input(lin) = NINT(temp)
2861	DO j = jts , MIN(jde-1,jte)
2862	DO i = its , MIN(ide-1,ite)
2863	temp = sm_input(i,lout+1,j)
2864	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2865	sm_input(i,lin+1,j) = temp
2866	END DO
2867	END DO
2868	END IF
2869	END DO innerm
2870	END DO outerm
2871	!tgs add IF
2872	IF ( flag_soilm000 .NE. 1 ) THEN
2873	DO j = jts , MIN(jde-1,jte)
2874	DO i = its , MIN(ide-1,ite)
2875	sm_input(i,1,j) = sm_input(i,2,j)
2876	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2877	END DO
2878	END DO
2879	ENDIF
2880
2881	! Sort the levels for liquid moisture.
2882
2883	outerw: DO lout = 1 , num_sw_levels_input-1
2884	innerw : DO lin = lout+1 , num_sw_levels_input
2885	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2886	temp = sw_levels_input(lout)
2887	sw_levels_input(lout) = sw_levels_input(lin)
2888	sw_levels_input(lin) = NINT(temp)
2889	DO j = jts , MIN(jde-1,jte)
2890	DO i = its , MIN(ide-1,ite)
2891	temp = sw_input(i,lout+1,j)
2892	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2893	sw_input(i,lin+1,j) = temp
2894	END DO
2895	END DO
2896	END IF
2897	END DO innerw
2898	END DO outerw
2899	IF ( num_sw_levels_input .GT. 1 ) THEN
2900	DO j = jts , MIN(jde-1,jte)
2901	DO i = its , MIN(ide-1,ite)
2902	sw_input(i,1,j) = sw_input(i,2,j)
2903	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2904	END DO
2905	END DO
2906	END IF
2907
2908	found_levels = .TRUE.
2909
2910	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2911
2912	found_levels = .FALSE.
2913
2914	ELSE
2915	CALL wrf_error_fatal ( &
2916	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
2917	END IF
2918
2919	! Is it OK to continue?
2920
2921	IF ( found_levels ) THEN
2922
2923	!tgs add IF
2924	IF ( flag_soilt000 .NE.1) THEN
2925	!tgs initialize from Noah data
2926	! Here are the levels that we have from the input for temperature. The input levels plus
2927	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2928
2929	zhave(1) = 0.
2930	DO l = 1 , num_st_levels_input
2931	zhave(l+1) = st_levels_input(l) / 100.
2932	END DO
2933	zhave(num_st_levels_input+2) = 300. / 100.
2934
2935	! Interpolate between the layers we have (zhave) and those that we want (zs).
2936
2937	z_wantt : DO lwant = 1 , num_soil_layers
2938	z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2939	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2940	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2941	DO j = jts , MIN(jde-1,jte)
2942	DO i = its , MIN(ide-1,ite)
2943	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2944	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2945	( zhave(lhave+1) - zhave(lhave) )
2946	END DO
2947	END DO
2948	EXIT z_havet
2949	END IF
2950	END DO z_havet
2951	END DO z_wantt
2952
2953	ELSE
2954
2955	!tgs initialize from RUCLSM data
2956	DO l = 1 , num_st_levels_input
2957	zhave(l) = st_levels_input(l) / 100.
2958	END DO
2959
2960	! Interpolate between the layers we have (zhave) and those that we want (zs).
2961
2962
2963	z_wantt_2 : DO lwant = 1 , num_soil_layers
2964	z_havet_2 : DO lhave = 1 , num_st_levels_input -1
2965	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2966	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2967	DO j = jts , MIN(jde-1,jte)
2968	DO i = its , MIN(ide-1,ite)
2969	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
2970	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2971	( zhave(lhave+1) - zhave(lhave) )
2972	END DO
2973	END DO
2974	EXIT z_havet_2
2975	END IF
2976	END DO z_havet_2
2977	END DO z_wantt_2
2978
2979	ENDIF
2980
2981
2982	IF ( flag_soilm000 .NE. 1 ) THEN
2983	!tgs initialize from Noah
2984	! Here are the levels that we have from the input for moisture. The input levels plus
2985	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2986	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2987	! as the most deep layer's value.
2988
2989	zhave(1) = 0.
2990	DO l = 1 , num_sm_levels_input
2991	zhave(l+1) = sm_levels_input(l) / 100.
2992	END DO
2993	zhave(num_sm_levels_input+2) = 300. / 100.
2994
2995	! Interpolate between the layers we have (zhave) and those that we want (zs).
2996
2997	z_wantm : DO lwant = 1 , num_soil_layers
2998	z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2999	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3000	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3001	DO j = jts , MIN(jde-1,jte)
3002	DO i = its , MIN(ide-1,ite)
3003	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
3004	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3005	( zhave(lhave+1) - zhave(lhave) )
3006	END DO
3007	END DO
3008	EXIT z_havem
3009	END IF
3010	END DO z_havem
3011	END DO z_wantm
3012
3013	ELSE
3014
3015	!tgs initialize from RUCLSM data
3016	DO l = 1 , num_sm_levels_input
3017	zhave(l) = sm_levels_input(l) / 100.
3018	END DO
3019
3020	! Interpolate between the layers we have (zhave) and those that we want (zs).
3021
3022	z_wantm_2 : DO lwant = 1 , num_soil_layers
3023	z_havem_2 : DO lhave = 1 , num_sm_levels_input -1
3024	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3025	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3026	DO j = jts , MIN(jde-1,jte)
3027	DO i = its , MIN(ide-1,ite)
3028	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3029	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3030	( zhave(lhave+1) - zhave(lhave) )
3031	END DO
3032	END DO
3033	EXIT z_havem_2
3034	END IF
3035	END DO z_havem_2
3036	END DO z_wantm_2
3037
3038	ENDIF
3039	! Any liquid soil moisture to worry about?
3040
3041	IF ( num_sw_levels_input .GT. 1 ) THEN
3042
3043	zhave(1) = 0.
3044	DO l = 1 , num_sw_levels_input
3045	zhave(l+1) = sw_levels_input(l) / 100.
3046	END DO
3047	zhave(num_sw_levels_input+2) = 300. / 100.
3048
3049	! Interpolate between the layers we have (zhave) and those that we want (zs).
3050
3051	z_wantw : DO lwant = 1 , num_soil_layers
3052	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
3053	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3054	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3055	DO j = jts , MIN(jde-1,jte)
3056	DO i = its , MIN(ide-1,ite)
3057	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
3058	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3059	( zhave(lhave+1) - zhave(lhave) )
3060	END DO
3061	END DO
3062	EXIT z_havew
3063	END IF
3064	END DO z_havew
3065	END DO z_wantw
3066
3067	END IF
3068
3069
3070	! Over water, put in reasonable values for soil temperature and moisture. These won't be
3071	! used, but they will make a more continuous plot.
3072
3073	IF ( flag_sst .EQ. 1 ) THEN
3074	DO j = jts , MIN(jde-1,jte)
3075	DO i = its , MIN(ide-1,ite)
3076	IF ( landmask(i,j) .LT. 0.5 ) THEN
3077	DO l = 1 , num_soil_layers
3078	tslb(i,l,j)= sst(i,j)
3079	!tgs add line for tsk
3080	tsk(i,j) = sst(i,j)
3081	smois(i,l,j)= 1.0
3082	sh2o (i,l,j)= 1.0
3083	END DO
3084	END IF
3085	END DO
3086	END DO
3087	ELSE
3088	DO j = jts , MIN(jde-1,jte)
3089	DO i = its , MIN(ide-1,ite)
3090	IF ( landmask(i,j) .LT. 0.5 ) THEN
3091	DO l = 1 , num_soil_layers
3092	tslb(i,l,j)= tsk(i,j)
3093	smois(i,l,j)= 1.0
3094	sh2o (i,l,j)= 1.0
3095	END DO
3096	END IF
3097	END DO
3098	END DO
3099	END IF
3100
3101	DEALLOCATE (zhave)
3102
3103	END IF
3104
3105	END SUBROUTINE init_soil_2_real
3106
3107	SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
3108	ivgtyp,isltyp,tslb,smois,tmn, &
3109	num_soil_layers, &
3110	ids,ide, jds,jde, kds,kde, &
3111	ims,ime, jms,jme, kms,kme, &
3112	its,ite, jts,jte, kts,kte )
3113
3114	IMPLICIT NONE
3115
3116	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
3117	ims,ime, jms,jme, kms,kme, &
3118	its,ite, jts,jte, kts,kte
3119
3120	INTEGER, INTENT(IN) ::num_soil_layers
3121
3122	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
3123
3124	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
3125
3126	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
3127
3128	INTEGER :: icm,jcm,itf,jtf
3129	INTEGER :: i,j,l
3130
3131	itf=min0(ite,ide-1)
3132	jtf=min0(jte,jde-1)
3133
3134	icm = ide/2
3135	jcm = jde/2
3136
3137	DO j=jts,jtf
3138	DO l=1,num_soil_layers
3139	DO i=its,itf
3140
3141	smois(i,1,j)=0.10
3142	smois(i,2,j)=0.10
3143	smois(i,3,j)=0.10
3144	smois(i,4,j)=0.10
3145
3146	tslb(i,1,j)=295.
3147	tslb(i,2,j)=297.
3148	tslb(i,3,j)=293.
3149	tslb(i,4,j)=293.
3150
3151	ENDDO
3152	ENDDO
3153	ENDDO
3154
3155	DO j=jts,jtf
3156	DO i=its,itf
3157	xland(i,j) = 2
3158	tmn(i,j) = 294.
3159	xice(i,j) = 0.
3160	vegfra(i,j) = 0.
3161	snow(i,j) = 0.
3162	canwat(i,j) = 0.
3163	ivgtyp(i,j) = 7
3164	isltyp(i,j) = 8
3165	ENDDO
3166	ENDDO
3167
3168	END SUBROUTINE init_soil_2_ideal
3169
3170	SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
3171	st_input , sm_input , landmask, sst, &
3172	zs , dzs , &
3173	st_levels_input , sm_levels_input , &
3174	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
3175	num_st_levels_alloc , num_sm_levels_alloc , &
3176	flag_sst , flag_soilt000 , flag_soilm000 , &
3177	ids , ide , jds , jde , kds , kde , &
3178	ims , ime , jms , jme , kms , kme , &
3179	its , ite , jts , jte , kts , kte )
3180
3181	IMPLICIT NONE
3182
3183	INTEGER , INTENT(IN) :: num_soil_layers , &
3184	num_st_levels_input , num_sm_levels_input , &
3185	num_st_levels_alloc , num_sm_levels_alloc , &
3186	ids , ide , jds , jde , kds , kde , &
3187	ims , ime , jms , jme , kms , kme , &
3188	its , ite , jts , jte , kts , kte
3189
3190	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
3191
3192	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
3193	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
3194
3195	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
3196	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
3197	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
3198
3199	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
3200	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
3201	REAL , DIMENSION(num_soil_layers) :: zs , dzs
3202
3203	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
3204
3205	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
3206
3207	INTEGER :: i , j , l , lout , lin , lwant , lhave
3208	REAL :: temp
3209
3210	! Allocate the soil layer array used for interpolating.
3211
3212	IF ( ( num_st_levels_input .LE. 0 ) .OR. &
3213	( num_sm_levels_input .LE. 0 ) ) THEN
3214	PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
3215	CALL wrf_error_fatal ( 'no soil data' )
3216	ELSE
3217	IF ( flag_soilt000 .eq. 1 ) THEN
3218	PRINT '(A)',' Assume RUC LSM 6-level input'
3219	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
3220	ELSE
3221	PRINT '(A)',' Assume non-RUC LSM input'
3222	ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
3223	END IF
3224	END IF
3225
3226	! Sort the levels for temperature.
3227
3228	outert : DO lout = 1 , num_st_levels_input-1
3229	innert : DO lin = lout+1 , num_st_levels_input
3230	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
3231	temp = st_levels_input(lout)
3232	st_levels_input(lout) = st_levels_input(lin)
3233	st_levels_input(lin) = NINT(temp)
3234	DO j = jts , MIN(jde-1,jte)
3235	DO i = its , MIN(ide-1,ite)
3236	temp = st_input(i,lout,j)
3237	st_input(i,lout,j) = st_input(i,lin,j)
3238	st_input(i,lin,j) = temp
3239	END DO
3240	END DO
3241	END IF
3242	END DO innert
3243	END DO outert
3244
3245	IF ( flag_soilt000 .NE. 1 ) THEN
3246	DO j = jts , MIN(jde-1,jte)
3247	DO i = its , MIN(ide-1,ite)
3248	st_input(i,1,j) = tsk(i,j)
3249	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
3250	END DO
3251	END DO
3252	END IF
3253
3254	! Sort the levels for moisture.
3255
3256	outerm: DO lout = 1 , num_sm_levels_input-1
3257	innerm : DO lin = lout+1 , num_sm_levels_input
3258	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
3259	temp = sm_levels_input(lout)
3260	sm_levels_input(lout) = sm_levels_input(lin)
3261	sm_levels_input(lin) = NINT(temp)
3262	DO j = jts , MIN(jde-1,jte)
3263	DO i = its , MIN(ide-1,ite)
3264	temp = sm_input(i,lout,j)
3265	sm_input(i,lout,j) = sm_input(i,lin,j)
3266	sm_input(i,lin,j) = temp
3267	END DO
3268	END DO
3269	END IF
3270	END DO innerm
3271	END DO outerm
3272
3273	IF ( flag_soilm000 .NE. 1 ) THEN
3274	DO j = jts , MIN(jde-1,jte)
3275	DO i = its , MIN(ide-1,ite)
3276	sm_input(i,1,j) = sm_input(i,2,j)
3277	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
3278	END DO
3279	END DO
3280	END IF
3281
3282	! Here are the levels that we have from the input for temperature.
3283
3284	IF ( flag_soilt000 .EQ. 1 ) THEN
3285	DO l = 1 , num_st_levels_input
3286	zhave(l) = st_levels_input(l) / 100.
3287	END DO
3288
3289	! Interpolate between the layers we have (zhave) and those that we want (zs).
3290
3291	z_wantt : DO lwant = 1 , num_soil_layers
3292	z_havet : DO lhave = 1 , num_st_levels_input -1
3293	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3294	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3295	DO j = jts , MIN(jde-1,jte)
3296	DO i = its , MIN(ide-1,ite)
3297	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3298	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3299	( zhave(lhave+1) - zhave(lhave) )
3300	END DO
3301	END DO
3302	EXIT z_havet
3303	END IF
3304	END DO z_havet
3305	END DO z_wantt
3306
3307	ELSE
3308
3309	zhave(1) = 0.
3310	DO l = 1 , num_st_levels_input
3311	zhave(l+1) = st_levels_input(l) / 100.
3312	END DO
3313	zhave(num_st_levels_input+2) = 300. / 100.
3314
3315	! Interpolate between the layers we have (zhave) and those that we want (zs).
3316
3317	z_wantt_2 : DO lwant = 1 , num_soil_layers
3318	z_havet_2 : DO lhave = 1 , num_st_levels_input +2
3319	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3320	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3321	DO j = jts , MIN(jde-1,jte)
3322	DO i = its , MIN(ide-1,ite)
3323	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3324	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3325	( zhave(lhave+1) - zhave(lhave) )
3326	END DO
3327	END DO
3328	EXIT z_havet_2
3329	END IF
3330	END DO z_havet_2
3331	END DO z_wantt_2
3332
3333	END IF
3334
3335	! Here are the levels that we have from the input for moisture.
3336
3337	IF ( flag_soilm000 .EQ. 1 ) THEN
3338	DO l = 1 , num_sm_levels_input
3339	zhave(l) = sm_levels_input(l) / 100.
3340	END DO
3341
3342	! Interpolate between the layers we have (zhave) and those that we want (zs).
3343
3344	z_wantm : DO lwant = 1 , num_soil_layers
3345	z_havem : DO lhave = 1 , num_sm_levels_input -1
3346	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3347	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3348	DO j = jts , MIN(jde-1,jte)
3349	DO i = its , MIN(ide-1,ite)
3350	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3351	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3352	( zhave(lhave+1) - zhave(lhave) )
3353	END DO
3354	END DO
3355	EXIT z_havem
3356	END IF
3357	END DO z_havem
3358	END DO z_wantm
3359
3360	ELSE
3361
3362	zhave(1) = 0.
3363	DO l = 1 , num_sm_levels_input
3364	zhave(l+1) = sm_levels_input(l) / 100.
3365	END DO
3366	zhave(num_sm_levels_input+2) = 300. / 100.
3367
3368	z_wantm_2 : DO lwant = 1 , num_soil_layers
3369	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
3370	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3371	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3372	DO j = jts , MIN(jde-1,jte)
3373	DO i = its , MIN(ide-1,ite)
3374	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3375	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3376	( zhave(lhave+1) - zhave(lhave) )
3377	END DO
3378	END DO
3379	EXIT z_havem_2
3380	END IF
3381	END DO z_havem_2
3382	END DO z_wantm_2
3383
3384	END IF
3385	! Over water, put in reasonable values for soil temperature and moisture. These won't be
3386	! used, but they will make a more continuous plot.
3387
3388	IF ( flag_sst .EQ. 1 ) THEN
3389	DO j = jts , MIN(jde-1,jte)
3390	DO i = its , MIN(ide-1,ite)
3391	IF ( landmask(i,j) .LT. 0.5 ) THEN
3392	DO l = 1 , num_soil_layers
3393	tslb(i,l,j) = sst(i,j)
3394	tsk(i,j) = sst(i,j)
3395	smois(i,l,j)= 1.0
3396	END DO
3397	END IF
3398	END DO
3399	END DO
3400	ELSE
3401	DO j = jts , MIN(jde-1,jte)
3402	DO i = its , MIN(ide-1,ite)
3403	IF ( landmask(i,j) .LT. 0.5 ) THEN
3404	DO l = 1 , num_soil_layers
3405	tslb(i,l,j)= tsk(i,j)
3406	smois(i,l,j)= 1.0
3407	END DO
3408	END IF
3409	END DO
3410	END DO
3411	END IF
3412
3413	DEALLOCATE (zhave)
3414
3415	END SUBROUTINE init_soil_3_real
3416
3417	END MODULE module_soil_pre
3418
3419	#endif
3420

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

Download in other formats:

Original Format