Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_soil_pre.F @ 1531

Last change on this file since 1531 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 161.5 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	CHARACTER (LEN=3) :: num_cat_count
9	INTEGER , PARAMETER , DIMENSION(0:300) :: ints = &
10	(/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, &
11	10, 11, 12, 13, 14, 15, 16, 17, 18, 19, &
12	20, 21, 22, 23, 24, 25, 26, 27, 28, 29, &
13	30, 31, 32, 33, 34, 35, 36, 37, 38, 39, &
14	40, 41, 42, 43, 44, 45, 46, 47, 48, 49, &
15	50, 51, 52, 53, 54, 55, 56, 57, 58, 59, &
16	60, 61, 62, 63, 64, 65, 66, 67, 68, 69, &
17	70, 71, 72, 73, 74, 75, 76, 77, 78, 79, &
18	80, 81, 82, 83, 84, 85, 86, 87, 88, 89, &
19	90, 91, 92, 93, 94, 95, 96, 97, 98, 99, &
20	100, 101, 102, 103, 104, 105, 106, 107, 108, 109, &
21	110, 111, 112, 113, 114, 115, 116, 117, 118, 119, &
22	120, 121, 122, 123, 124, 125, 126, 127, 128, 129, &
23	130, 131, 132, 133, 134, 135, 136, 137, 138, 139, &
24	140, 141, 142, 143, 144, 145, 146, 147, 148, 149, &
25	150, 151, 152, 153, 154, 155, 156, 157, 158, 159, &
26	160, 161, 162, 163, 164, 165, 166, 167, 168, 169, &
27	170, 171, 172, 173, 174, 175, 176, 177, 178, 179, &
28	180, 181, 182, 183, 184, 185, 186, 187, 188, 189, &
29	190, 191, 192, 193, 194, 195, 196, 197, 198, 199, &
30	200, 201, 202, 203, 204, 205, 206, 207, 208, 209, &
31	210, 211, 212, 213, 214, 215, 216, 217, 218, 219, &
32	220, 221, 222, 223, 224, 225, 226, 227, 228, 229, &
33	230, 231, 232, 233, 234, 235, 236, 237, 238, 239, &
34	240, 241, 242, 243, 244, 245, 246, 247, 248, 249, &
35	250, 251, 252, 253, 254, 255, 256, 257, 258, 259, &
36	260, 261, 262, 263, 264, 265, 266, 267, 268, 269, &
37	270, 271, 272, 273, 274, 275, 276, 277, 278, 279, &
38	280, 281, 282, 283, 284, 285, 286, 287, 288, 289, &
39	290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300 /)
40
41	! Excluded middle processing
42
43	LOGICAL , SAVE :: hold_ups
44	INTEGER , SAVE :: em_width
45
46	LOGICAL , EXTERNAL :: skip_middle_points_t
47
48	CONTAINS
49
50	SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
51	xland , landusef , isltyp , soilcat , soilctop , &
52	soilcbot , tmn , &
53	seaice_threshold , &
54	fractional_seaice, &
55	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
56	iswater , isice , &
57	scheme , &
58	ids , ide , jds , jde , kds , kde , &
59	ims , ime , jms , jme , kms , kme , &
60	its , ite , jts , jte , kts , kte )
61
62	IMPLICIT NONE
63
64	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
65	ims , ime , jms , jme , kms , kme , &
66	its , ite , jts , jte , kts , kte , &
67	iswater , isice
68	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
69
70	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
71	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
72	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
73	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
74	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
75	vegcat, xland , soilcat , tmn
76	REAL , INTENT(IN) :: seaice_threshold
77
78	INTEGER :: i , j , num_seaice_changes , loop
79	CHARACTER (LEN=132) :: message
80
81	INTEGER, INTENT(IN) :: fractional_seaice
82	REAL :: XICE_THRESHOLD
83
84	IF ( FRACTIONAL_SEAICE == 0 ) THEN
85	xice_threshold = 0.5
86	ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
87	xice_threshold = 0.02
88	ENDIF
89
90	num_seaice_changes = 0
91	fix_seaice : SELECT CASE ( scheme )
92
93	CASE ( SLABSCHEME )
94	DO j = jts , MIN(jde-1,jte)
95	DO i = its , MIN(ide-1,ite)
96	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
97	IF ( xice(i,j) .GT. 200.0 ) THEN
98	xice(i,j) = 0.
99	num_seaice_changes = num_seaice_changes + 1
100	END IF
101	END DO
102	END DO
103	IF ( num_seaice_changes .GT. 0 ) THEN
104	WRITE ( message , FMT='(A,I6)' ) &
105	'Total pre number of sea ice locations removed (due to FLAG values) = ', &
106	num_seaice_changes
107	CALL wrf_debug ( 0 , message )
108	END IF
109	num_seaice_changes = 0
110	DO j = jts , MIN(jde-1,jte)
111	DO i = its , MIN(ide-1,ite)
112	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
113	IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
114	( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
115	IF ( FRACTIONAL_SEAICE == 0 ) THEN
116	xice(i,j) = 1.0
117	ENDIF
118	num_seaice_changes = num_seaice_changes + 1
119	if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
120	vegcat(i,j)=isice
121	ivgtyp(i,j)=isice
122	lu_index(i,j)=isice
123	landmask(i,j)=1.
124	xland(i,j)=1.
125	DO loop=1,num_veg_cat
126	landusef(i,loop,j)=0.
127	END DO
128	landusef(i,ivgtyp(i,j),j)=1.
129
130	isltyp(i,j) = 16
131	soilcat(i,j)=isltyp(i,j)
132	DO loop=1,num_soil_top_cat
133	soilctop(i,loop,j)=0
134	END DO
135	DO loop=1,num_soil_bot_cat
136	soilcbot(i,loop,j)=0
137	END DO
138	soilctop(i,isltyp(i,j),j)=1.
139	soilcbot(i,isltyp(i,j),j)=1.
140	ELSE
141	xice(i,j) = 0.0
142	END IF
143	END DO
144	END DO
145	IF ( num_seaice_changes .GT. 0 ) THEN
146	WRITE ( message , FMT='(A,I6)' ) &
147	'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
148	CALL wrf_debug ( 0 , message )
149	END IF
150
151	CASE ( LSMSCHEME , RUCLSMSCHEME )
152	num_seaice_changes = 0
153	DO j = jts , MIN(jde-1,jte)
154	DO i = its , MIN(ide-1,ite)
155	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
156	IF ( landmask(i,j) .GT. 0.5 ) THEN
157	if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
158	xice(i,j) = 0.
159	END IF
160	END DO
161	END DO
162	IF ( num_seaice_changes .GT. 0 ) THEN
163	WRITE ( message , FMT='(A,I6)' ) &
164	'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
165	CALL wrf_debug ( 0 , message )
166	END IF
167
168	END SELECT fix_seaice
169
170	END SUBROUTINE adjust_for_seaice_pre
171
172	SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
173	xland , landusef , isltyp , soilcat , soilctop , &
174	soilcbot , tmn , vegfra , &
175	tslb , smois , sh2o , &
176	seaice_threshold , &
177	sst , flag_sst , &
178	fractional_seaice, &
179	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
180	num_soil_layers , &
181	iswater , isice , &
182	scheme , &
183	ids , ide , jds , jde , kds , kde , &
184	ims , ime , jms , jme , kms , kme , &
185	its , ite , jts , jte , kts , kte )
186
187	IMPLICIT NONE
188
189	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
190	ims , ime , jms , jme , kms , kme , &
191	its , ite , jts , jte , kts , kte , &
192	iswater , isice
193	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
194	INTEGER , INTENT(IN) :: num_soil_layers
195
196	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
197	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
198	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
199	REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
200	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN):: sst
201	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
202	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
203	vegcat, xland , soilcat , tmn , &
204	tsk_old , vegfra
205	INTEGER , INTENT(IN) :: flag_sst
206	REAL , INTENT(IN) :: seaice_threshold
207	REAL :: total_depth , mid_point_depth
208
209	INTEGER :: i , j , num_seaice_changes , loop
210	CHARACTER (LEN=132) :: message
211
212
213	INTEGER, INTENT(IN) :: fractional_seaice
214	real :: xice_threshold
215
216	IF ( FRACTIONAL_SEAICE == 0 ) THEN
217	xice_threshold = 0.5
218	ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
219	xice_threshold = 0.02
220	ENDIF
221	num_seaice_changes = 0
222	fix_seaice : SELECT CASE ( scheme )
223
224	CASE ( SLABSCHEME )
225
226	CASE ( LSMSCHEME )
227	DO j = jts , MIN(jde-1,jte)
228	DO i = its , MIN(ide-1,ite)
229	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
230	IF ( xice(i,j) .GT. 200.0 ) THEN
231	xice(i,j) = 0.
232	num_seaice_changes = num_seaice_changes + 1
233	END IF
234	END DO
235	END DO
236	IF ( num_seaice_changes .GT. 0 ) THEN
237	WRITE ( message , FMT='(A,I6)' ) &
238	'Total post number of sea ice locations removed (due to FLAG values) = ', &
239	num_seaice_changes
240	CALL wrf_debug ( 0 , message )
241	END IF
242	num_seaice_changes = 0
243	DO j = jts , MIN(jde-1,jte)
244	DO i = its , MIN(ide-1,ite)
245	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
246	IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
247	( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
248	tsk(i,j) = tsk_old(i,j)
249	END IF
250	IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
251	( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
252	print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
253	CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
254	ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
255	( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
256	IF ( FRACTIONAL_SEAICE == 0 ) THEN
257	xice(i,j) = 1.0
258	ENDIF
259	num_seaice_changes = num_seaice_changes + 1
260	if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
261	vegcat(i,j)=isice
262	ivgtyp(i,j)=isice
263	lu_index(i,j)=isice
264	landmask(i,j)=1.
265	xland(i,j)=1.
266	vegfra(i,j)=0.
267	DO loop=1,num_veg_cat
268	landusef(i,loop,j)=0.
269	END DO
270	landusef(i,ivgtyp(i,j),j)=1.
271
272	tsk_old(i,j) = tsk(i,j)
273
274	isltyp(i,j) = 16
275	soilcat(i,j)=isltyp(i,j)
276	DO loop=1,num_soil_top_cat
277	soilctop(i,loop,j)=0
278	END DO
279	DO loop=1,num_soil_bot_cat
280	soilcbot(i,loop,j)=0
281	END DO
282	soilctop(i,isltyp(i,j),j)=1.
283	soilcbot(i,isltyp(i,j),j)=1.
284
285	total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
286	DO loop = 1,num_soil_layers
287	mid_point_depth=(total_depth/num_soil_layers)/2. + &
288	(loop-1)*(total_depth/num_soil_layers)
289	tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
290	mid_point_depth*tmn(i,j) ) / total_depth
291	END DO
292
293	DO loop=1,num_soil_layers
294	smois(i,loop,j) = 1.0
295	sh2o(i,loop,j) = 0.0
296	END DO
297	ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
298	xice(i,j) = 0.
299	END IF
300	END DO
301	END DO
302	IF ( num_seaice_changes .GT. 0 ) THEN
303	WRITE ( message , FMT='(A,I6)' ) &
304	'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
305	CALL wrf_debug ( 0 , message )
306	END IF
307
308	CASE ( RUCLSMSCHEME )
309	DO j = jts , MIN(jde-1,jte)
310	DO i = its , MIN(ide-1,ite)
311	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
312	IF ( xice(i,j) .GT. 200.0 ) THEN
313	xice(i,j) = 0.
314	num_seaice_changes = num_seaice_changes + 1
315	END IF
316	END DO
317	END DO
318	IF ( num_seaice_changes .GT. 0 ) THEN
319	WRITE ( message , FMT='(A,I6)' ) &
320	'Total post number of sea ice locations removed (due to FLAG values) = ', &
321	num_seaice_changes
322	CALL wrf_debug ( 0 , message )
323	END IF
324	num_seaice_changes = 0
325	DO j = jts , MIN(jde-1,jte)
326	DO i = its , MIN(ide-1,ite)
327	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
328	IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
329	( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
330	tsk(i,j) = tsk_old(i,j)
331	END IF
332	IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
333	( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
334	print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
335	CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
336	ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
337	( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
338	IF ( FRACTIONAL_SEAICE == 0 ) THEN
339	xice(i,j) = 1.0
340	ELSE
341	xice(i,j)=max(0.25,xice(i,j))
342	ENDIF
343	num_seaice_changes = num_seaice_changes + 1
344	if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
345	vegcat(i,j)=isice
346	ivgtyp(i,j)=isice
347	lu_index(i,j)=isice
348	landmask(i,j)=1.
349	xland(i,j)=1.
350	vegfra(i,j)=0.
351	DO loop=1,num_veg_cat
352	landusef(i,loop,j)=0.
353	END DO
354	landusef(i,ivgtyp(i,j),j)=1.
355
356	!tgs - compute blended sea ice/water skin temperature
357	if(flag_sst.eq.1) then
358	tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j))) &
359	+(1-xice(i,j))*sst(i,j)
360	else
361	tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j))) &
362	+(1-xice(i,j))*tsk(i,j)
363	endif
364	tsk_old(i,j) = tsk(i,j)
365
366	isltyp(i,j) = 16
367	soilcat(i,j)=isltyp(i,j)
368	DO loop=1,num_soil_top_cat
369	soilctop(i,loop,j)=0
370	END DO
371	DO loop=1,num_soil_bot_cat
372	soilcbot(i,loop,j)=0
373	END DO
374	soilctop(i,isltyp(i,j),j)=1.
375	soilcbot(i,isltyp(i,j),j)=1.
376
377	total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
378	tslb(i,1,j) = tsk(i,j)
379	tslb(i,num_soil_layers,j) = tmn(i,j)
380	DO loop = 2,num_soil_layers-1
381	mid_point_depth=(total_depth/num_soil_layers)/4. + &
382	(loop-2)*(total_depth/num_soil_layers)
383	tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
384	mid_point_depth*tmn(i,j) ) / total_depth
385	END DO
386
387	DO loop=1,num_soil_layers
388	smois(i,loop,j) = 1.0
389	sh2o(i,loop,j) = 0.0
390	END DO
391	ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
392	xice(i,j) = 0.
393	END IF
394	END DO
395	END DO
396	IF ( num_seaice_changes .GT. 0 ) THEN
397	WRITE ( message , FMT='(A,I6)' ) &
398	'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
399	CALL wrf_debug ( 0 , message )
400	END IF
401
402
403	END SELECT fix_seaice
404
405	END SUBROUTINE adjust_for_seaice_post
406
407	SUBROUTINE process_percent_cat_new ( landmask , &
408	landuse_frac , soil_top_cat , soil_bot_cat , &
409	isltyp , ivgtyp , &
410	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
411	ids , ide , jds , jde , kds , kde , &
412	ims , ime , jms , jme , kms , kme , &
413	its , ite , jts , jte , kts , kte , &
414	iswater )
415
416	IMPLICIT NONE
417
418	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
419	ims , ime , jms , jme , kms , kme , &
420	its , ite , jts , jte , kts , kte , &
421	iswater
422	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
423	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
424	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
425	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
426	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
427	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
428
429	INTEGER :: i , j , l , ll, dominant_index
430	REAL :: dominant_value
431
432	#ifdef WRF_CHEM
433	! REAL :: lwthresh = .99
434	REAL :: lwthresh = .50
435	#else
436	REAL :: lwthresh = .50
437	#endif
438
439	INTEGER , PARAMETER :: iswater_soil = 14
440	INTEGER :: iforce
441	CHARACTER (LEN=132) :: message
442	CHARACTER(LEN=256) :: mminlu
443	LOGICAL :: aggregate_lu
444	integer :: change_water , change_land
445	change_water = 0
446	change_land = 0
447
448	! Sanity check on the 50/50 points
449
450	DO j = jts , MIN(jde-1,jte)
451	DO i = its , MIN(ide-1,ite)
452	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
453	dominant_value = landuse_frac(i,iswater,j)
454	IF ( dominant_value .EQ. lwthresh ) THEN
455	DO l = 1 , num_veg_cat
456	IF ( l .EQ. iswater ) CYCLE
457	IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
458	PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
459	landuse_frac(i,l,j) = lwthresh - .01
460	landuse_frac(i,iswater,j) = lwthresh + 0.01
461	ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
462	PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
463	landuse_frac(i,l,j) = lwthresh + .01
464	landuse_frac(i,iswater,j) = lwthresh - 0.01
465	END IF
466	END DO
467	END IF
468	END DO
469	END DO
470
471	! Compute the aggregate of the vegetation/land use categories. Lump all of the grasses together,
472	! all of the shrubs, all of the trees, etc. Choose the correct set of available land use
473	! categories. Also, make sure that the user wants to actually avail themselves of aforementioned
474	! opportunity, as mayhaps they don't.
475
476	CALL nl_get_mminlu ( 1 , mminlu )
477	CALL nl_get_aggregate_lu ( 1 , aggregate_lu )
478	IF ( aggregate_lu ) THEN
479	DO j = jts , MIN(jde-1,jte)
480	DO i = its , MIN(ide-1,ite)
481	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
482	CALL aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu(1:4) )
483	END DO
484	END DO
485	END IF
486
487	! Compute the dominant VEGETATION INDEX.
488
489	DO j = jts , MIN(jde-1,jte)
490	DO i = its , MIN(ide-1,ite)
491	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
492	dominant_value = landuse_frac(i,1,j)
493	dominant_index = 1
494	DO l = 2 , num_veg_cat
495	IF ( l .EQ. iswater ) THEN
496	! wait a bit
497	ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
498	dominant_value = landuse_frac(i,l,j)
499	dominant_index = l
500	END IF
501	END DO
502	IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
503	dominant_value = landuse_frac(i,iswater,j)
504	dominant_index = iswater
505	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
506	( landmask(i,j) .LT. 0.5) .AND. &
507	( dominant_value .EQ. lwthresh) ) THEN
508	dominant_value = landuse_frac(i,iswater,j)
509	dominant_index = iswater
510	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
511	( landmask(i,j) .GT. 0.5) .AND. &
512	( dominant_value .EQ. lwthresh) ) THEN
513	!no op
514	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
515	( dominant_value .LT. lwthresh ) ) THEN
516	dominant_value = landuse_frac(i,iswater,j)
517	dominant_index = iswater
518	END IF
519	IF ( dominant_index .EQ. iswater ) THEN
520	if(landmask(i,j).gt.lwthresh) then
521	!print *,'changing to water at point ',i,j
522	!WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
523	!WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
524	!CALL wrf_debug(1,message)
525	!WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
526	!CALL wrf_debug(1,message)
527	change_water=change_water+1
528	endif
529	landmask(i,j) = 0
530	ELSE IF ( dominant_index .NE. iswater ) THEN
531	if(landmask(i,j).lt.lwthresh) then
532	!print *,'changing to land at point ',i,j
533	!WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
534	!WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
535	!CALL wrf_debug(1,message)
536	!WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
537	!CALL wrf_debug(1,message)
538	change_land=change_land+1
539	endif
540	landmask(i,j) = 1
541	END IF
542	ivgtyp(i,j) = dominant_index
543	END DO
544	END DO
545
546	! Compute the dominant SOIL TEXTURE INDEX, TOP.
547
548	iforce = 0
549	DO i = its , MIN(ide-1,ite)
550	DO j = jts , MIN(jde-1,jte)
551	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
552	dominant_value = soil_top_cat(i,1,j)
553	dominant_index = 1
554	IF ( landmask(i,j) .GT. lwthresh ) THEN
555	DO l = 2 , num_soil_top_cat
556	IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
557	dominant_value = soil_top_cat(i,l,j)
558	dominant_index = l
559	END IF
560	END DO
561	IF ( dominant_value .LT. 0.01 ) THEN
562	iforce = iforce + 1
563	WRITE ( message , FMT = '(A,I4,I4)' ) &
564	'based on landuse, changing soil to land at point ',i,j
565	CALL wrf_debug(1,message)
566	WRITE ( num_cat_count , FMT = '(I3)' ) num_soil_top_cat
567	WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) (ints(l),l=1,num_soil_top_cat)
568	CALL wrf_debug(1,message)
569	WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) &
570	((nint(soil_top_cat(i,ints(l),j)*100)), l=1,num_soil_top_cat)
571	CALL wrf_debug(1,message)
572	dominant_index = 8
573	END IF
574	ELSE
575	dominant_index = iswater_soil
576	END IF
577	isltyp(i,j) = dominant_index
578	END DO
579	END DO
580
581	if(iforce.ne.0)then
582	WRITE(message,FMT='(A,I4,A,I6)' ) &
583	'forcing artificial silty clay loam at ',iforce,' points, out of ',&
584	(MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
585	CALL wrf_debug(0,message)
586	endif
587	print *,'LAND CHANGE = ',change_land
588	print *,'WATER CHANGE = ',change_water
589
590	END SUBROUTINE process_percent_cat_new
591
592	SUBROUTINE process_soil_real ( tsk , tmn , tavgsfc, &
593	landmask , sst , ht, toposoil, &
594	st_input , sm_input , sw_input , &
595	st_levels_input , sm_levels_input , sw_levels_input , &
596	zs , dzs , tslb , smois , sh2o , &
597	flag_sst , flag_tavgsfc, flag_soilhgt, &
598	flag_soil_layers, flag_soil_levels, &
599	ids , ide , jds , jde , kds , kde , &
600	ims , ime , jms , jme , kms , kme , &
601	its , ite , jts , jte , kts , kte , &
602	sf_surface_physics , num_soil_layers , real_data_init_type , &
603	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
604	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
605
606	IMPLICIT NONE
607
608	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
609	ims , ime , jms , jme , kms , kme , &
610	its , ite , jts , jte , kts , kte , &
611	sf_surface_physics , num_soil_layers , real_data_init_type , &
612	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
613	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
614
615	INTEGER , INTENT(IN) :: flag_sst, flag_tavgsfc
616	INTEGER , INTENT(IN) :: flag_soil_layers, flag_soil_levels, flag_soilhgt
617
618	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
619
620	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
621	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
622	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
623	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
624	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
625	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
626
627	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
628	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
629	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tavgsfc, ht, toposoil
630	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk, tmn
631
632	INTEGER :: i , j , k, l , dominant_index , num_soil_cat , num_veg_cat, closest_layer
633	REAL :: dominant_value, closest_depth, diff_cm
634	REAL , ALLOCATABLE , DIMENSION(:) :: depth ! Soil layer thicknesses (cm)
635	REAL, PARAMETER :: get_temp_closest_to = 30. ! use soil temperature closest to this depth (cm)
636	REAL, PARAMETER :: something_big = 1.e6 ! Initialize closest depth as something big (cm)
637	INTEGER :: something_far = 1000 ! Soil array index far away
638	CHARACTER (LEN=132) :: message
639
640	! Case statement for tmn initialization
641	! Need to have a reasonable default value for annual mean deeeeep soil temperature
642	! For sf_surface_physics = 1, we want to use close to a 30 cm value
643	! for the bottom level of the soil temps.
644	! NOTE: We are assuming that soil_layers are the same for each grid point
645
646	fix_bottom_level_for_temp : SELECT CASE ( sf_surface_physics )
647	CASE (SLABSCHEME)
648	IF ( flag_tavgsfc .EQ. 1 ) THEN
649	CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
650	DO j = jts , MIN(jde-1,jte)
651	DO i = its , MIN(ide-1,ite)
652	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
653	tmn(i,j) = tavgsfc(i,j)
654	END DO
655	END DO
656	ELSE
657	! Look for soil temp close to 30 cm
658	closest_layer = something_far
659	closest_depth = something_big
660	DO k = 1, num_st_levels_input
661	diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
662	IF ( diff_cm < closest_depth ) THEN
663	closest_depth = diff_cm
664	closest_layer = k
665	END IF
666	END DO
667	IF ( closest_layer == something_far ) THEN
668	CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
669	CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
670	ELSE
671	write(message, FMT='(A,F7.2,A,I3)')&
672	'Soil temperature closest to ',get_temp_closest_to, &
673	' at level ',st_levels_input(closest_layer)
674	CALL wrf_debug ( 0 , message )
675	DO j = jts , MIN(jde-1,jte)
676	DO i = its , MIN(ide-1,ite)
677	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
678	tmn(i,j) = st_input(i,closest_layer+1,j)
679	END DO
680	END DO
681	END IF
682	END IF
683
684	CASE (LSMSCHEME)
685
686	CASE (RUCLSMSCHEME)
687
688	CASE (PXLSMSCHEME)
689	! When the input data from met_em is in layers, there is an additional level added to the beginning
690	! of the array to define the surface, which is why we add the extra value (closest_layer+1)
691	IF ( flag_tavgsfc .EQ. 1 ) THEN
692	CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
693	DO j = jts , MIN(jde-1,jte)
694	DO i = its , MIN(ide-1,ite)
695	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
696	tmn(i,j) = tavgsfc(i,j)
697	END DO
698	END DO
699	ELSE
700	! Look for soil temp close to 30 cm
701	closest_layer = num_st_levels_input+1
702	closest_depth = something_big
703	DO k = 1, num_st_levels_input
704	diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
705	IF ( diff_cm < closest_depth ) THEN
706	closest_depth = diff_cm
707	closest_layer = k
708	END IF
709	END DO
710	IF ( closest_layer == num_st_levels_input + 1 ) THEN
711	CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
712	CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
713	ELSE
714	write(message, FMT='(A,F7.2,A,I3)')&
715	'Soil temperature closest to ',get_temp_closest_to, &
716	' at level ',st_levels_input(closest_layer)
717	CALL wrf_debug ( 0 , message )
718	DO j = jts , MIN(jde-1,jte)
719	DO i = its , MIN(ide-1,ite)
720	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
721	tmn(i,j) = st_input(i,closest_layer+1,j)
722	END DO
723	END DO
724	END IF
725	END IF
726
727	#if 0
728	! Loop over layers and do a weighted mean
729	IF ( ALLOCATED ( depth ) ) DEALLOCATE ( depth )
730	ALLOCATE ( depth(num_st_levels_input) )
731	IF ( flag_soil_layers == 1 ) THEN
732	DO k = num_st_levels_input, 2, -1
733	depth(k) = st_levels_input(k) - st_levels_input(k-1)
734	END DO
735	depth(1) = st_levels_input(1)
736	ELSE IF ( flag_soil_levels == 1 ) THEN
737	DO k = 2, num_st_levels_input
738	depth(k) = st_levels_input(k) - st_levels_input(k-1)
739	END DO
740	depth(1) = 0.
741	END IF
742
743	DO j = jts , MIN(jde-1,jte)
744	DO i = its , MIN(ide-1,ite)
745	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
746	tmn(i,j) = 0.
747	DO k = 1, num_st_levels_input
748	tmn(i,j) = tmn(i,j) + depth(k) * st_input(i,k,j)
749	END DO
750	END DO
751	END DO
752	DEALLOCATE ( depth )
753
754	#endif
755
756	END SELECT fix_bottom_level_for_temp
757
758	! Adjust the various soil temperature values depending on the difference in
759	! elevation between the current model's elevation and the incoming data's
760	! orography.
761
762	adjust_soil : SELECT CASE ( sf_surface_physics )
763
764	CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME, PXLSMSCHEME )
765	CALL adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ht , &
766	toposoil , landmask , st_input, st_levels_input, &
767	flag_soilhgt , flag_tavgsfc , &
768	flag_soil_layers , flag_soil_levels, &
769	num_st_levels_input, num_st_levels_alloc, &
770	ids , ide , jds , jde , kds , kde , &
771	ims , ime , jms , jme , kms , kme , &
772	its , ite , jts , jte , kts , kte )
773
774	END SELECT adjust_soil
775
776	! Initialize the soil depth, and the soil temperature and moisture.
777
778	IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
779	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
780	CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
781	landmask , sst , flag_sst , &
782	ids , ide , jds , jde , kds , kde , &
783	ims , ime , jms , jme , kms , kme , &
784	its , ite , jts , jte , kts , kte )
785
786	ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
787	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
788	CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
789	st_input , sm_input , sw_input , landmask , sst , &
790	zs , dzs , &
791	st_levels_input , sm_levels_input , sw_levels_input , &
792	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
793	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
794	flag_sst , flag_soil_layers , flag_soil_levels , &
795	ids , ide , jds , jde , kds , kde , &
796	ims , ime , jms , jme , kms , kme , &
797	its , ite , jts , jte , kts , kte )
798	ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
799	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
800	CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
801	st_input , sm_input , landmask , sst , &
802	zs , dzs , &
803	st_levels_input , sm_levels_input , &
804	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
805	num_st_levels_alloc , num_sm_levels_alloc , &
806	flag_sst , flag_soil_layers , flag_soil_levels , &
807	ids , ide , jds , jde , kds , kde , &
808	ims , ime , jms , jme , kms , kme , &
809	its , ite , jts , jte , kts , kte )
810	ELSE IF ( ( sf_surface_physics .EQ. PXLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
811	CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
812	CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
813	st_input , sm_input , sw_input, landmask , sst , &
814	zs , dzs , &
815	st_levels_input , sm_levels_input , sw_levels_input, &
816	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
817	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
818	flag_sst , flag_soil_layers , flag_soil_levels , &
819	ids , ide , jds , jde , kds , kde , &
820	ims , ime , jms , jme , kms , kme , &
821	its , ite , jts , jte , kts , kte )
822	END IF
823
824	END SUBROUTINE process_soil_real
825
826	SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
827	ivgtyp,isltyp,tslb,smois, &
828	tsk,tmn,zs,dzs, &
829	num_soil_layers, &
830	sf_surface_physics , &
831	ids,ide, jds,jde, kds,kde,&
832	ims,ime, jms,jme, kms,kme,&
833	its,ite, jts,jte, kts,kte )
834
835	IMPLICIT NONE
836
837	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
838	ims,ime, jms,jme, kms,kme, &
839	its,ite, jts,jte, kts,kte
840
841	INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
842
843	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
844
845	REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
846
847	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
848	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
849	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
850
851	! Local variables.
852
853	INTEGER :: itf,jtf
854
855	itf=MIN(ite,ide-1)
856	jtf=MIN(jte,jde-1)
857
858	IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
859	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
860	CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
861	ivgtyp,zs,dzs,num_soil_layers, &
862	ids,ide, jds,jde, kds,kde, &
863	ims,ime, jms,jme, kms,kme, &
864	its,ite, jts,jte, kts,kte )
865	ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
866	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
867	CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
868	ivgtyp,isltyp,tslb,smois,tmn, &
869	num_soil_layers, &
870	ids,ide, jds,jde, kds,kde, &
871	ims,ime, jms,jme, kms,kme, &
872	its,ite, jts,jte, kts,kte )
873	ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
874	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
875
876	END IF
877
878	END SUBROUTINE process_soil_ideal
879
880	SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ter , &
881	toposoil , landmask , st_input , st_levels_input, &
882	flag_toposoil , flag_tavgsfc , &
883	flag_soil_layers , flag_soil_levels, &
884	num_st_levels_input, num_st_levels_alloc, &
885	ids , ide , jds , jde , kds , kde , &
886	ims , ime , jms , jme , kms , kme , &
887	its , ite , jts , jte , kts , kte )
888
889	IMPLICIT NONE
890
891	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
892	ims , ime , jms , jme , kms , kme , &
893	its , ite , jts , jte , kts , kte
894	INTEGER , INTENT(IN) :: num_st_levels_input, num_st_levels_alloc
895
896	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
897	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
898	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
899	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(IN) :: st_levels_input
900
901	INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
902	INTEGER , INTENT(IN) :: flag_soil_layers , flag_soil_levels
903
904	INTEGER :: i , j, k , st_near_sfc
905
906	REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
907
908	! Adjust the annual mean temperature as if it is based on from a sea-level elevation
909	! if the value used is from the actual annula mean data set. If the input field to
910	! be used for tmn is one of the first-guess input temp fields, need to do an adjustment
911	! only on the diff in topo from the model terrain and the first-guess terrain.
912
913	SELECT CASE ( sf_surface_physics )
914
915	CASE (LSMSCHEME)
916	DO j = jts , MIN(jde-1,jte)
917	DO i = its , MIN(ide-1,ite)
918	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
919	IF (landmask(i,j) .GT. 0.5 ) THEN
920	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
921	END IF
922	END DO
923	END DO
924
925	CASE (RUCLSMSCHEME)
926	DO j = jts , MIN(jde-1,jte)
927	DO i = its , MIN(ide-1,ite)
928	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
929	IF (landmask(i,j) .GT. 0.5 ) THEN
930	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
931	END IF
932	END DO
933	END DO
934
935	END SELECT
936
937
938	! Do we have a soil field with which to modify soil temperatures?
939
940	IF ( flag_toposoil .EQ. 1 ) THEN
941
942	DO j = jts , MIN(jde-1,jte)
943	DO i = its , MIN(ide-1,ite)
944
945	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
946	! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
947	! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
948	! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
949	! the difference between the soil elevation and the terrain is greater than 3 km means
950	! that the soil data is either all zeros or that the data are inconsistent. Any of these
951	! three conditions is grievous enough to induce a WRF fatality. However, if they are at
952	! a water point, then we can safely ignore them.
953
954	soil_elev_min_val = toposoil(i,j)
955	soil_elev_max_val = toposoil(i,j)
956	soil_elev_min_dif = ter(i,j) - toposoil(i,j)
957	soil_elev_max_dif = ter(i,j) - toposoil(i,j)
958
959	IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
960	CYCLE
961	ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
962	!print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
963	cycle
964	! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
965	ENDIF
966
967	IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
968	CYCLE
969	ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
970	print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
971	cycle
972	CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
973	ENDIF
974
975	IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
976	( landmask(i,j) .LT. 0.5 ) ) THEN
977	CYCLE
978	ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
979	( landmask(i,j) .GT. 0.5 ) ) THEN
980	print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
981	cycle
982	CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
983	ENDIF
984
985	! For each of the fields that we would like to modify, check to see if it came in from the SI.
986	! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
987	! are not at a water point.
988
989	IF (landmask(i,j) .GT. 0.5 ) THEN
990	IF ( sf_surface_physics .EQ. SLABSCHEME ) THEN
991	st_near_sfc = 0 ! Check if there is a soil layer above 40 cm
992	DO k = 1, num_st_levels_input
993	IF ( st_levels_input(k) .LE. 40 ) THEN
994	st_near_sfc = 1
995	END IF
996	END DO
997	IF ( ( flag_tavgsfc == 1 ) .OR. ( st_near_sfc == 1 ) ) THEN
998	tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
999	ELSE
1000	tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
1001	END IF
1002	END IF
1003
1004	tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1005
1006	IF ( flag_soil_layers == 1 ) THEN
1007	DO k = 2, num_st_levels_input+1
1008	st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1009	END DO
1010	ELSE
1011	DO k = 1, num_st_levels_input
1012	st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
1013	END DO
1014	END IF
1015
1016	END IF
1017	END DO
1018	END DO
1019
1020	END IF
1021
1022	END SUBROUTINE adjust_soil_temp_new
1023
1024
1025	SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
1026
1027	IMPLICIT NONE
1028
1029	INTEGER, INTENT(IN) :: num_soil_layers
1030
1031	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1032
1033	INTEGER :: l
1034
1035	! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
1036	! The distance from the ground level to the midpoint of the layer is given by zs.
1037
1038	! ------- Ground Level ---------- \|\| \|\| \|\| \|\|
1039	! \|\| \|\| \|\| \|\| zs(1) = 0.005 m
1040	! -- -- -- -- -- -- -- -- -- \|\| \|\| \|\| \/
1041	! \|\| \|\| \|\|
1042	! ----------------------------------- \|\| \|\| \|\| \/ dzs(1) = 0.01 m
1043	! \|\| \|\| \|\|
1044	! \|\| \|\| \|\| zs(2) = 0.02
1045	! -- -- -- -- -- -- -- -- -- \|\| \|\| \/
1046	! \|\| \|\|
1047	! \|\| \|\|
1048	! ----------------------------------- \|\| \|\| \/ dzs(2) = 0.02 m
1049	! \|\| \|\|
1050	! \|\| \|\|
1051	! \|\| \|\|
1052	! \|\| \|\| zs(3) = 0.05
1053	! -- -- -- -- -- -- -- -- -- \|\| \/
1054	! \|\|
1055	! \|\|
1056	! \|\|
1057	! \|\|
1058	! ----------------------------------- \/ dzs(3) = 0.04 m
1059
1060	IF ( num_soil_layers .NE. 5 ) THEN
1061	PRINT '(A)','Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
1062	CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
1063	END IF
1064
1065	dzs(1)=.01
1066	zs(1)=.5*dzs(1)
1067
1068	DO l=2,num_soil_layers
1069	dzs(l)=2*dzs(l-1)
1070	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
1071	ENDDO
1072
1073	END SUBROUTINE init_soil_depth_1
1074
1075	SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
1076
1077	IMPLICIT NONE
1078
1079	INTEGER, INTENT(IN) :: num_soil_layers
1080
1081	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1082
1083	INTEGER :: l
1084
1085	dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
1086
1087	IF ( num_soil_layers .NE. 4 ) THEN
1088	PRINT '(A)','Usually, the LSM uses 4 layers. Change this in the namelist.'
1089	CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
1090	END IF
1091
1092	zs(1)=.5*dzs(1)
1093
1094	DO l=2,num_soil_layers
1095	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
1096	ENDDO
1097
1098	END SUBROUTINE init_soil_depth_2
1099
1100	SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
1101
1102	IMPLICIT NONE
1103
1104	INTEGER, INTENT(IN) :: num_soil_layers
1105
1106	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1107
1108	INTEGER :: l
1109
1110	CHARACTER (LEN=132) :: message
1111
1112	! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
1113	! ZS is specified in the namelist: num_soil_layers = 6 or 9.
1114	! Other options with number of levels are possible, but
1115	! WRF users should change consistently the namelist entry with the
1116	! ZS array in this subroutine.
1117
1118	IF ( num_soil_layers .EQ. 6) THEN
1119	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
1120	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
1121	ELSEIF ( num_soil_layers .EQ. 9) THEN
1122	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
1123	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
1124	ENDIF
1125
1126	IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
1127	write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
1128	CALL wrf_error_fatal ( message )
1129	END IF
1130
1131	END SUBROUTINE init_soil_depth_3
1132
1133	SUBROUTINE init_soil_depth_7 ( zs , dzs , num_soil_layers )
1134
1135	IMPLICIT NONE
1136
1137	INTEGER, INTENT(IN) :: num_soil_layers
1138
1139	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
1140
1141	INTEGER :: l
1142
1143	dzs = (/ 0.01 , 0.99 /)
1144
1145	IF ( num_soil_layers .NE. 2 ) THEN
1146	PRINT '(A)','Usually, the PX LSM uses 2 layers. Change this in the namelist.'
1147	CALL wrf_error_fatal ( 'PXLSM_uses_2_layers' )
1148	END IF
1149
1150	zs(1) = 0.5 * dzs(1)
1151	zs(2) = dzs(1) + 0.5 * dzs(2)
1152
1153	END SUBROUTINE init_soil_depth_7
1154
1155	SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
1156	num_soil_layers , real_data_init_type , &
1157	landmask , sst , flag_sst , &
1158	ids , ide , jds , jde , kds , kde , &
1159	ims , ime , jms , jme , kms , kme , &
1160	its , ite , jts , jte , kts , kte )
1161
1162	IMPLICIT NONE
1163
1164	INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
1165	ids , ide , jds , jde , kds , kde , &
1166	ims , ime , jms , jme , kms , kme , &
1167	its , ite , jts , jte , kts , kte
1168
1169	INTEGER , INTENT(IN) :: flag_sst
1170
1171	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1172	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
1173
1174	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1175
1176	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
1177
1178	INTEGER :: i , j , l
1179
1180	! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
1181	! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
1182	! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
1183	! annual mean temperature.
1184
1185	DO j = jts , MIN(jde-1,jte)
1186	DO i = its , MIN(ide-1,ite)
1187	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1188	IF ( landmask(i,j) .GT. 0.5 ) THEN
1189	DO l = 1 , num_soil_layers
1190	tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
1191	tmn(i,j) * ( zs( l) - zs(1) ) ) / &
1192	( zs(num_soil_layers) - zs(1) )
1193	END DO
1194	ELSE
1195	IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
1196	DO l = 1 , num_soil_layers
1197	tslb(i,l,j)= sst(i,j)
1198	END DO
1199	ELSE
1200	DO l = 1 , num_soil_layers
1201	tslb(i,l,j)= tsk(i,j)
1202	END DO
1203	END IF
1204	END IF
1205	END DO
1206	END DO
1207
1208	END SUBROUTINE init_soil_1_real
1209
1210	SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
1211	ivgtyp,ZS,DZS,num_soil_layers, &
1212	ids,ide, jds,jde, kds,kde, &
1213	ims,ime, jms,jme, kms,kme, &
1214	its,ite, jts,jte, kts,kte )
1215
1216	IMPLICIT NONE
1217
1218	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1219	ims,ime, jms,jme, kms,kme, &
1220	its,ite, jts,jte, kts,kte
1221
1222	INTEGER, INTENT(IN ) :: num_soil_layers
1223
1224	REAL, DIMENSION( ims:ime , num_soil_layers , jms:jme ), INTENT(OUT) :: tslb
1225	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
1226	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
1227
1228	REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
1229
1230	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
1231
1232	! Lcal variables.
1233
1234	INTEGER :: l,j,i,itf,jtf
1235
1236	itf=MIN(ite,ide-1)
1237	jtf=MIN(jte,jde-1)
1238
1239	IF (num_soil_layers.NE.1)THEN
1240	DO j=jts,jtf
1241	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1242	DO l=1,num_soil_layers
1243	DO i=its,itf
1244	tslb(i,l,j)=( tsk(i,j)(zs(num_soil_layers)-zs(l)) + tmn(i,j)(zs(l)-zs(1)) ) / &
1245	( zs(num_soil_layers)-zs(1) )
1246	ENDDO
1247	ENDDO
1248	ENDDO
1249	ENDIF
1250
1251	! DO j=jts,jtf
1252	! DO i=its,itf
1253	! xland(i,j) = 2
1254	! ivgtyp(i,j) = 7
1255	! ENDDO
1256	! ENDDO
1257
1258	END SUBROUTINE init_soil_1_ideal
1259
1260	SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
1261	st_input , sm_input , sw_input , landmask , sst , &
1262	zs , dzs , &
1263	st_levels_input , sm_levels_input , sw_levels_input , &
1264	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1265	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1266	flag_sst , flag_soil_layers , flag_soil_levels , &
1267	ids , ide , jds , jde , kds , kde , &
1268	ims , ime , jms , jme , kms , kme , &
1269	its , ite , jts , jte , kts , kte )
1270
1271	IMPLICIT NONE
1272
1273	INTEGER , INTENT(IN) :: num_soil_layers , &
1274	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1275	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1276	ids , ide , jds , jde , kds , kde , &
1277	ims , ime , jms , jme , kms , kme , &
1278	its , ite , jts , jte , kts , kte
1279
1280	INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
1281
1282	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1283	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1284	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1285
1286	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1287	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1288	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1289	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1290
1291	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1292	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1293	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1294
1295	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1296
1297	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1298
1299	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1300	REAL :: temp
1301	LOGICAL :: found_levels
1302
1303	CHARACTER (LEN=132) :: message
1304
1305	! Are there any soil temp and moisture levels - ya know, they are mandatory.
1306
1307	num = num_st_levels_input * num_sm_levels_input
1308
1309	IF ( num .GE. 1 ) THEN
1310
1311	! Ordered levels that we have data for.
1312
1313	!tgs add option to initialize from RUCLSM
1314	IF ( flag_soil_levels == 1 ) THEN
1315	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1316	CALL wrf_message ( message )
1317	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
1318	ELSE
1319	write(message, FMT='(A)') ' Assume Noah LSM input'
1320	CALL wrf_message ( message )
1321	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1322	END IF
1323
1324
1325	! Sort the levels for temperature.
1326
1327	outert : DO lout = 1 , num_st_levels_input-1
1328	innert : DO lin = lout+1 , num_st_levels_input
1329	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1330	temp = st_levels_input(lout)
1331	st_levels_input(lout) = st_levels_input(lin)
1332	st_levels_input(lin) = NINT(temp)
1333	DO j = jts , MIN(jde-1,jte)
1334	DO i = its , MIN(ide-1,ite)
1335	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1336	temp = st_input(i,lout+1,j)
1337	st_input(i,lout+1,j) = st_input(i,lin+1,j)
1338	st_input(i,lin+1,j) = temp
1339	END DO
1340	END DO
1341	END IF
1342	END DO innert
1343	END DO outert
1344	!tgs add IF
1345	IF ( flag_soil_layers == 1 ) THEN
1346	DO j = jts , MIN(jde-1,jte)
1347	DO i = its , MIN(ide-1,ite)
1348	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1349	st_input(i,1,j) = tsk(i,j)
1350	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1351	END DO
1352	END DO
1353	ENDIF
1354
1355	! Sort the levels for moisture.
1356
1357	outerm: DO lout = 1 , num_sm_levels_input-1
1358	innerm : DO lin = lout+1 , num_sm_levels_input
1359	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1360	temp = sm_levels_input(lout)
1361	sm_levels_input(lout) = sm_levels_input(lin)
1362	sm_levels_input(lin) = NINT(temp)
1363	DO j = jts , MIN(jde-1,jte)
1364	DO i = its , MIN(ide-1,ite)
1365	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1366	temp = sm_input(i,lout+1,j)
1367	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1368	sm_input(i,lin+1,j) = temp
1369	END DO
1370	END DO
1371	END IF
1372	END DO innerm
1373	END DO outerm
1374	!tgs add IF
1375	IF ( flag_soil_layers == 1 ) THEN
1376	DO j = jts , MIN(jde-1,jte)
1377	DO i = its , MIN(ide-1,ite)
1378	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1379	sm_input(i,1,j) = sm_input(i,2,j)
1380	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1381	END DO
1382	END DO
1383	ENDIF
1384
1385	! Sort the levels for liquid moisture.
1386
1387	outerw: DO lout = 1 , num_sw_levels_input-1
1388	innerw : DO lin = lout+1 , num_sw_levels_input
1389	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1390	temp = sw_levels_input(lout)
1391	sw_levels_input(lout) = sw_levels_input(lin)
1392	sw_levels_input(lin) = NINT(temp)
1393	DO j = jts , MIN(jde-1,jte)
1394	DO i = its , MIN(ide-1,ite)
1395	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1396	temp = sw_input(i,lout+1,j)
1397	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1398	sw_input(i,lin+1,j) = temp
1399	END DO
1400	END DO
1401	END IF
1402	END DO innerw
1403	END DO outerw
1404	IF ( num_sw_levels_input .GT. 1 ) THEN
1405	DO j = jts , MIN(jde-1,jte)
1406	DO i = its , MIN(ide-1,ite)
1407	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1408	sw_input(i,1,j) = sw_input(i,2,j)
1409	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1410	END DO
1411	END DO
1412	END IF
1413
1414	found_levels = .TRUE.
1415
1416	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
1417
1418	found_levels = .FALSE.
1419
1420	ELSE
1421	CALL wrf_error_fatal ( &
1422	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
1423	END IF
1424
1425	! Is it OK to continue?
1426
1427	IF ( found_levels ) THEN
1428
1429	!tgs: Here are the levels that we have from the input for temperature.
1430
1431	IF ( flag_soil_levels == 1 ) THEN
1432	DO l = 1 , num_st_levels_input
1433	zhave(l) = st_levels_input(l) / 100.
1434	END DO
1435
1436	! Interpolate between the layers we have (zhave) and those that we want (zs).
1437
1438	z_wantt : DO lwant = 1 , num_soil_layers
1439	z_havet : DO lhave = 1 , num_st_levels_input -1
1440	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1441	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1442	DO j = jts , MIN(jde-1,jte)
1443	DO i = its , MIN(ide-1,ite)
1444	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1445	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1446	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1447	( zhave(lhave+1) - zhave(lhave) )
1448	END DO
1449	END DO
1450	EXIT z_havet
1451	END IF
1452	END DO z_havet
1453	END DO z_wantt
1454
1455	ELSE
1456
1457	! Here are the levels that we have from the input for temperature. The input levels plus
1458	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1459
1460	zhave(1) = 0.
1461	DO l = 1 , num_st_levels_input
1462	zhave(l+1) = st_levels_input(l) / 100.
1463	END DO
1464	zhave(num_st_levels_input+2) = 300. / 100.
1465
1466	! Interpolate between the layers we have (zhave) and those that we want (zs).
1467
1468	z_wantt_2: DO lwant = 1 , num_soil_layers
1469	z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
1470	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1471	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1472	DO j = jts , MIN(jde-1,jte)
1473	DO i = its , MIN(ide-1,ite)
1474	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1475	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1476	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1477	( zhave(lhave+1) - zhave(lhave) )
1478	END DO
1479	END DO
1480	EXIT z_havet_2
1481	END IF
1482	END DO z_havet_2
1483	END DO z_wantt_2
1484
1485	END IF
1486
1487	!tgs:
1488	IF ( flag_soil_levels == 1 ) THEN
1489	DO l = 1 , num_sm_levels_input
1490	zhave(l) = sm_levels_input(l) / 100.
1491	END DO
1492
1493	! Interpolate between the layers we have (zhave) and those that we want (zs).
1494
1495	z_wantm : DO lwant = 1 , num_soil_layers
1496	z_havem : DO lhave = 1 , num_sm_levels_input -1
1497	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1498	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1499	DO j = jts , MIN(jde-1,jte)
1500	DO i = its , MIN(ide-1,ite)
1501	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1502	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1503	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1504	( zhave(lhave+1) - zhave(lhave) )
1505	END DO
1506	END DO
1507	EXIT z_havem
1508	END IF
1509	END DO z_havem
1510	END DO z_wantm
1511
1512	ELSE
1513	! Here are the levels that we have from the input for moisture. The input levels plus
1514	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
1515	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
1516	! as the most deep layer's value.
1517
1518	zhave(1) = 0.
1519	DO l = 1 , num_sm_levels_input
1520	zhave(l+1) = sm_levels_input(l) / 100.
1521	END DO
1522	zhave(num_sm_levels_input+2) = 300. / 100.
1523
1524	! Interpolate between the layers we have (zhave) and those that we want (zs).
1525
1526	z_wantm_2 : DO lwant = 1 , num_soil_layers
1527	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
1528	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1529	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1530	DO j = jts , MIN(jde-1,jte)
1531	DO i = its , MIN(ide-1,ite)
1532	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1533	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1534	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1535	( zhave(lhave+1) - zhave(lhave) )
1536	END DO
1537	END DO
1538	EXIT z_havem_2
1539	END IF
1540	END DO z_havem_2
1541	END DO z_wantm_2
1542	ENDIF
1543
1544	! Any liquid soil moisture to worry about?
1545
1546	IF ( num_sw_levels_input .GT. 1 ) THEN
1547
1548	zhave(1) = 0.
1549	DO l = 1 , num_sw_levels_input
1550	zhave(l+1) = sw_levels_input(l) / 100.
1551	END DO
1552	zhave(num_sw_levels_input+2) = 300. / 100.
1553
1554	! Interpolate between the layers we have (zhave) and those that we want (zs).
1555
1556	z_wantw : DO lwant = 1 , num_soil_layers
1557	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1558	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1559	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1560	DO j = jts , MIN(jde-1,jte)
1561	DO i = its , MIN(ide-1,ite)
1562	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1563	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
1564	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1565	( zhave(lhave+1) - zhave(lhave) )
1566	END DO
1567	END DO
1568	EXIT z_havew
1569	END IF
1570	END DO z_havew
1571	END DO z_wantw
1572
1573	END IF
1574
1575
1576	! Over water, put in reasonable values for soil temperature and moisture. These won't be
1577	! used, but they will make a more continuous plot.
1578
1579	IF ( flag_sst .EQ. 1 ) THEN
1580	DO j = jts , MIN(jde-1,jte)
1581	DO i = its , MIN(ide-1,ite)
1582	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1583	IF ( landmask(i,j) .LT. 0.5 ) THEN
1584	DO l = 1 , num_soil_layers
1585	tslb(i,l,j)= sst(i,j)
1586	smois(i,l,j)= 1.0
1587	sh2o (i,l,j)= 1.0
1588	END DO
1589	END IF
1590	END DO
1591	END DO
1592	ELSE
1593	DO j = jts , MIN(jde-1,jte)
1594	DO i = its , MIN(ide-1,ite)
1595	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1596	IF ( landmask(i,j) .LT. 0.5 ) THEN
1597	DO l = 1 , num_soil_layers
1598	tslb(i,l,j)= tsk(i,j)
1599	smois(i,l,j)= 1.0
1600	sh2o (i,l,j)= 1.0
1601	END DO
1602	END IF
1603	END DO
1604	END DO
1605	END IF
1606
1607	DEALLOCATE (zhave)
1608
1609	END IF
1610
1611	END SUBROUTINE init_soil_2_real
1612
1613	SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
1614	ivgtyp,isltyp,tslb,smois,tmn, &
1615	num_soil_layers, &
1616	ids,ide, jds,jde, kds,kde, &
1617	ims,ime, jms,jme, kms,kme, &
1618	its,ite, jts,jte, kts,kte )
1619
1620	IMPLICIT NONE
1621
1622	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
1623	ims,ime, jms,jme, kms,kme, &
1624	its,ite, jts,jte, kts,kte
1625
1626	INTEGER, INTENT(IN) ::num_soil_layers
1627
1628	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
1629
1630	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: xland
1631	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: snow, canwat, xice, vegfra, tmn
1632
1633	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1634
1635	INTEGER :: icm,jcm,itf,jtf
1636	INTEGER :: i,j,l
1637
1638	itf=min0(ite,ide-1)
1639	jtf=min0(jte,jde-1)
1640
1641	icm = ide/2
1642	jcm = jde/2
1643
1644	DO j=jts,jtf
1645	DO l=1,num_soil_layers
1646	DO i=its,itf
1647	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1648	if (xland(i,j) .lt. 1.5) then
1649	smois(i,1,j)=0.30
1650	smois(i,2,j)=0.30
1651	smois(i,3,j)=0.30
1652	smois(i,4,j)=0.30
1653
1654	tslb(i,1,j)=290.
1655	tslb(i,2,j)=290.
1656	tslb(i,3,j)=290.
1657	tslb(i,4,j)=290.
1658	endif
1659	ENDDO
1660	ENDDO
1661	ENDDO
1662
1663	END SUBROUTINE init_soil_2_ideal
1664
1665	SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
1666	st_input , sm_input , landmask, sst, &
1667	zs , dzs , &
1668	st_levels_input , sm_levels_input , &
1669	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
1670	num_st_levels_alloc , num_sm_levels_alloc , &
1671	flag_sst , flag_soil_layers , flag_soil_levels , &
1672	ids , ide , jds , jde , kds , kde , &
1673	ims , ime , jms , jme , kms , kme , &
1674	its , ite , jts , jte , kts , kte )
1675
1676	IMPLICIT NONE
1677
1678	INTEGER , INTENT(IN) :: num_soil_layers , &
1679	num_st_levels_input , num_sm_levels_input , &
1680	num_st_levels_alloc , num_sm_levels_alloc , &
1681	ids , ide , jds , jde , kds , kde , &
1682	ims , ime , jms , jme , kms , kme , &
1683	its , ite , jts , jte , kts , kte
1684
1685	INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
1686
1687	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1688	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1689
1690	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1691	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1692	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1693
1694	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1695	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1696	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1697
1698	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
1699
1700	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1701
1702	INTEGER :: i , j , l , lout , lin , lwant , lhave, k
1703	REAL :: temp
1704
1705	CHARACTER (LEN=132) :: message
1706
1707	! Allocate the soil layer array used for interpolating.
1708
1709	IF ( ( num_st_levels_input .LE. 0 ) .OR. &
1710	( num_sm_levels_input .LE. 0 ) ) THEN
1711	write (message, FMT='(A)')&
1712	'No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
1713	CALL wrf_error_fatal ( message )
1714	ELSE
1715	IF ( flag_soil_levels == 1 ) THEN
1716	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1717	CALL wrf_message ( message )
1718	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
1719	ELSE
1720	write(message, FMT='(A)') ' Assume non-RUC LSM input'
1721	CALL wrf_message ( message )
1722	ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
1723	END IF
1724	END IF
1725
1726	! Sort the levels for temperature.
1727
1728	outert : DO lout = 1 , num_st_levels_input-1
1729	innert : DO lin = lout+1 , num_st_levels_input
1730	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1731	temp = st_levels_input(lout)
1732	st_levels_input(lout) = st_levels_input(lin)
1733	st_levels_input(lin) = NINT(temp)
1734	DO j = jts , MIN(jde-1,jte)
1735	DO i = its , MIN(ide-1,ite)
1736	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1737	temp = st_input(i,lout,j)
1738	st_input(i,lout,j) = st_input(i,lin,j)
1739	st_input(i,lin,j) = temp
1740	END DO
1741	END DO
1742	END IF
1743	END DO innert
1744	END DO outert
1745
1746	IF ( flag_soil_layers == 1 ) THEN
1747	DO j = jts , MIN(jde-1,jte)
1748	DO i = its , MIN(ide-1,ite)
1749	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1750	st_input(i,1,j) = tsk(i,j)
1751	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1752	END DO
1753	END DO
1754	END IF
1755
1756	! Sort the levels for moisture.
1757
1758	outerm: DO lout = 1 , num_sm_levels_input-1
1759	innerm : DO lin = lout+1 , num_sm_levels_input
1760	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1761	temp = sm_levels_input(lout)
1762	sm_levels_input(lout) = sm_levels_input(lin)
1763	sm_levels_input(lin) = NINT(temp)
1764	DO j = jts , MIN(jde-1,jte)
1765	DO i = its , MIN(ide-1,ite)
1766	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1767	temp = sm_input(i,lout,j)
1768	sm_input(i,lout,j) = sm_input(i,lin,j)
1769	sm_input(i,lin,j) = temp
1770	END DO
1771	END DO
1772	END IF
1773	END DO innerm
1774	END DO outerm
1775
1776	IF ( flag_soil_layers == 1 ) THEN
1777	DO j = jts , MIN(jde-1,jte)
1778	DO i = its , MIN(ide-1,ite)
1779	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1780	sm_input(i,1,j) = (sm_input(i,2,j)-sm_input(i,3,j))/ &
1781	(st_levels_input(2)-st_levels_input(1))*st_levels_input(1)+ &
1782	sm_input(i,2,j)
1783	! sm_input(i,1,j) = sm_input(i,2,j)
1784	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1785	END DO
1786	END DO
1787	END IF
1788
1789	! Here are the levels that we have from the input for temperature.
1790
1791	IF ( flag_soil_levels == 1 ) THEN
1792	DO l = 1 , num_st_levels_input
1793	zhave(l) = st_levels_input(l) / 100.
1794	END DO
1795
1796	! Interpolate between the layers we have (zhave) and those that we want (zs).
1797
1798
1799	z_wantt : DO lwant = 1 , num_soil_layers
1800	z_havet : DO lhave = 1 , num_st_levels_input -1
1801	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1802	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1803	DO j = jts , MIN(jde-1,jte)
1804	DO i = its , MIN(ide-1,ite)
1805	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1806	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1807	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1808	( zhave(lhave+1) - zhave(lhave) )
1809	END DO
1810	END DO
1811	EXIT z_havet
1812	END IF
1813	END DO z_havet
1814	END DO z_wantt
1815
1816	ELSE
1817
1818	zhave(1) = 0.
1819	DO l = 1 , num_st_levels_input
1820	zhave(l+1) = st_levels_input(l) / 100.
1821	END DO
1822	zhave(num_st_levels_input+2) = 300. / 100.
1823
1824	! Interpolate between the layers we have (zhave) and those that we want (zs).
1825
1826	z_wantt_2 : DO lwant = 1 , num_soil_layers
1827	z_havet_2 : DO lhave = 1 , num_st_levels_input +2
1828	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1829	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1830	DO j = jts , MIN(jde-1,jte)
1831	DO i = its , MIN(ide-1,ite)
1832	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1833	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1834	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1835	( zhave(lhave+1) - zhave(lhave) )
1836	END DO
1837	END DO
1838	EXIT z_havet_2
1839	END IF
1840	END DO z_havet_2
1841	END DO z_wantt_2
1842
1843	END IF
1844
1845	! Here are the levels that we have from the input for moisture.
1846
1847	IF ( flag_soil_levels .EQ. 1 ) THEN
1848	DO l = 1 , num_sm_levels_input
1849	zhave(l) = sm_levels_input(l) / 100.
1850	END DO
1851
1852	! Interpolate between the layers we have (zhave) and those that we want (zs).
1853
1854	z_wantm : DO lwant = 1 , num_soil_layers
1855	z_havem : DO lhave = 1 , num_sm_levels_input -1
1856	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1857	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1858	DO j = jts , MIN(jde-1,jte)
1859	DO i = its , MIN(ide-1,ite)
1860	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1861	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1862	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1863	( zhave(lhave+1) - zhave(lhave) )
1864	END DO
1865	END DO
1866	EXIT z_havem
1867	END IF
1868	END DO z_havem
1869	END DO z_wantm
1870
1871	ELSE
1872
1873	zhave(1) = 0.
1874	DO l = 1 , num_sm_levels_input
1875	zhave(l+1) = sm_levels_input(l) / 100.
1876	END DO
1877	zhave(num_sm_levels_input+2) = 300. / 100.
1878
1879	z_wantm_2 : DO lwant = 1 , num_soil_layers
1880	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
1881	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
1882	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1883	DO j = jts , MIN(jde-1,jte)
1884	DO i = its , MIN(ide-1,ite)
1885	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1886	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
1887	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
1888	( zhave(lhave+1) - zhave(lhave) )
1889	END DO
1890	END DO
1891	EXIT z_havem_2
1892	END IF
1893	END DO z_havem_2
1894	END DO z_wantm_2
1895
1896	END IF
1897	! Over water, put in reasonable values for soil temperature and moisture. These won't be
1898	! used, but they will make a more continuous plot.
1899
1900	IF ( flag_sst .EQ. 1 ) THEN
1901	DO j = jts , MIN(jde-1,jte)
1902	DO i = its , MIN(ide-1,ite)
1903	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1904	IF ( landmask(i,j) .LT. 0.5 ) THEN
1905	DO l = 1 , num_soil_layers
1906	tslb(i,l,j) = sst(i,j)
1907	tsk(i,j) = sst(i,j)
1908	smois(i,l,j)= 1.0
1909	END DO
1910	END IF
1911	END DO
1912	END DO
1913	ELSE
1914	DO j = jts , MIN(jde-1,jte)
1915	DO i = its , MIN(ide-1,ite)
1916	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1917	IF ( landmask(i,j) .LT. 0.5 ) THEN
1918	DO l = 1 , num_soil_layers
1919	tslb(i,l,j)= tsk(i,j)
1920	smois(i,l,j)= 1.0
1921	END DO
1922	END IF
1923	END DO
1924	END DO
1925	END IF
1926
1927	DEALLOCATE (zhave)
1928
1929	END SUBROUTINE init_soil_3_real
1930	SUBROUTINE init_soil_7_real ( tsk , tmn , smois , sh2o , tslb , &
1931	st_input , sm_input , sw_input , landmask , sst , &
1932	zs , dzs , &
1933	st_levels_input , sm_levels_input , sw_levels_input , &
1934	num_soil_layers , num_st_levels_input , &
1935	num_sm_levels_input , num_sw_levels_input , &
1936	num_st_levels_alloc , num_sm_levels_alloc , &
1937	num_sw_levels_alloc , &
1938	flag_sst , flag_soil_layers , flag_soil_levels , &
1939	ids , ide , jds , jde , kds , kde , &
1940	ims , ime , jms , jme , kms , kme , &
1941	its , ite , jts , jte , kts , kte )
1942
1943	! for soil temperature and moisture initialization for the PX LSM
1944
1945	IMPLICIT NONE
1946
1947	INTEGER , INTENT(IN) :: num_soil_layers , &
1948	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1949	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1950	ids , ide , jds , jde , kds , kde , &
1951	ims , ime , jms , jme , kms , kme , &
1952	its , ite , jts , jte , kts , kte
1953
1954	INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
1955
1956	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1957	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1958	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1959
1960	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1961	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1962	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1963	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1964
1965	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1966	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1967	REAL , DIMENSION(num_soil_layers) :: zs , dzs
1968
1969	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1970
1971	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1972
1973	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1974	REAL :: temp
1975	LOGICAL :: found_levels
1976
1977	! Are there any soil temp and moisture levels - ya know, they are mandatory.
1978
1979	num = num_st_levels_input * num_sm_levels_input
1980
1981	IF ( num .GE. 1 ) THEN
1982
1983	! Ordered levels that we have data for.
1984
1985	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1986
1987	! Sort the levels for temperature.
1988	outert : DO lout = 1 , num_st_levels_input-1
1989	innert : DO lin = lout+1 , num_st_levels_input
1990	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1991	temp = st_levels_input(lout)
1992	st_levels_input(lout) = st_levels_input(lin)
1993	st_levels_input(lin) = NINT(temp)
1994	DO j = jts , MIN(jde-1,jte)
1995	DO i = its , MIN(ide-1,ite)
1996	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
1997	temp = st_input(i,lout+1,j)
1998	st_input(i,lout+1,j) = st_input(i,lin+1,j)
1999	st_input(i,lin+1,j) = temp
2000	END DO
2001	END DO
2002	END IF
2003	END DO innert
2004	END DO outert
2005	DO j = jts , MIN(jde-1,jte)
2006	DO i = its , MIN(ide-1,ite)
2007	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2008	st_input(i,1,j) = tsk(i,j)
2009	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2010	END DO
2011	END DO
2012
2013	! Sort the levels for moisture.
2014
2015	outerm: DO lout = 1 , num_sm_levels_input-1
2016	innerm : DO lin = lout+1 , num_sm_levels_input
2017	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2018	temp = sm_levels_input(lout)
2019	sm_levels_input(lout) = sm_levels_input(lin)
2020	sm_levels_input(lin) = NINT(temp)
2021	DO j = jts , MIN(jde-1,jte)
2022	DO i = its , MIN(ide-1,ite)
2023	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2024	temp = sm_input(i,lout+1,j)
2025	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2026	sm_input(i,lin+1,j) = temp
2027	END DO
2028	END DO
2029	END IF
2030	END DO innerm
2031	END DO outerm
2032	DO j = jts , MIN(jde-1,jte)
2033	DO i = its , MIN(ide-1,ite)
2034	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2035	sm_input(i,1,j) = sm_input(i,2,j)
2036	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2037	END DO
2038	END DO
2039
2040	! Sort the levels for liquid moisture.
2041
2042	outerw: DO lout = 1 , num_sw_levels_input-1
2043	innerw : DO lin = lout+1 , num_sw_levels_input
2044	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2045	temp = sw_levels_input(lout)
2046	sw_levels_input(lout) = sw_levels_input(lin)
2047	sw_levels_input(lin) = NINT(temp)
2048	DO j = jts , MIN(jde-1,jte)
2049	DO i = its , MIN(ide-1,ite)
2050	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2051	temp = sw_input(i,lout+1,j)
2052	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2053	sw_input(i,lin+1,j) = temp
2054	END DO
2055	END DO
2056	END IF
2057	END DO innerw
2058	END DO outerw
2059	IF ( num_sw_levels_input .GT. 1 ) THEN
2060	DO j = jts , MIN(jde-1,jte)
2061	DO i = its , MIN(ide-1,ite)
2062	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2063	sw_input(i,1,j) = sw_input(i,2,j)
2064	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2065	END DO
2066	END DO
2067	END IF
2068
2069	found_levels = .TRUE.
2070
2071	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
2072
2073	found_levels = .FALSE.
2074
2075	ELSE
2076	CALL wrf_error_fatal ( &
2077	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for PX LSM.' )
2078	END IF
2079
2080	! Is it OK to continue?
2081
2082	IF ( found_levels ) THEN
2083
2084	! Here are the levels that we have from the input for temperature. The input levels plus
2085	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2086
2087	zhave(1) = 0.
2088	DO l = 1 , num_st_levels_input
2089	zhave(l+1) = st_levels_input(l) / 100.
2090	END DO
2091	zhave(num_st_levels_input+2) = 300. / 100.
2092
2093	! Interpolate between the layers we have (zhave) and those that we want (zs).
2094
2095	z_wantt : DO lwant = 1 , num_soil_layers
2096	z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2097	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2098	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2099	DO j = jts , MIN(jde-1,jte)
2100	DO i = its , MIN(ide-1,ite)
2101	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2102	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2103	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2104	( zhave(lhave+1) - zhave(lhave) )
2105	END DO
2106	END DO
2107	EXIT z_havet
2108	END IF
2109	END DO z_havet
2110	END DO z_wantt
2111
2112	! Here are the levels that we have from the input for moisture. The input levels plus
2113	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
2114	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
2115	! as the most deep layer's value.
2116
2117	zhave(1) = 0.
2118	DO l = 1 , num_sm_levels_input
2119	zhave(l+1) = sm_levels_input(l) / 100.
2120	END DO
2121	zhave(num_sm_levels_input+2) = 300. / 100.
2122
2123	! Interpolate between the layers we have (zhave) and those that we want (zs).
2124
2125	z_wantm : DO lwant = 1 , num_soil_layers
2126	z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2127	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2128	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2129	DO j = jts , MIN(jde-1,jte)
2130	DO i = its , MIN(ide-1,ite)
2131	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2132	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2133	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2134	( zhave(lhave+1) - zhave(lhave) )
2135	END DO
2136	END DO
2137	EXIT z_havem
2138	END IF
2139	END DO z_havem
2140	END DO z_wantm
2141
2142	! Any liquid soil moisture to worry about?
2143
2144	IF ( num_sw_levels_input .GT. 1 ) THEN
2145
2146	zhave(1) = 0.
2147	DO l = 1 , num_sw_levels_input
2148	zhave(l+1) = sw_levels_input(l) / 100.
2149	END DO
2150	zhave(num_sw_levels_input+2) = 300. / 100.
2151
2152	! Interpolate between the layers we have (zhave) and those that we want (zs).
2153
2154	z_wantw : DO lwant = 1 , num_soil_layers
2155	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
2156	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
2157	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2158	DO j = jts , MIN(jde-1,jte)
2159	DO i = its , MIN(ide-1,ite)
2160	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2161	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
2162	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
2163	( zhave(lhave+1) - zhave(lhave) )
2164	END DO
2165	END DO
2166	EXIT z_havew
2167	END IF
2168	END DO z_havew
2169	END DO z_wantw
2170
2171	END IF
2172
2173
2174	! Over water, put in reasonable values for soil temperature and moisture. These won't be
2175	! used, but they will make a more continuous plot.
2176
2177	DO j = jts , MIN(jde-1,jte)
2178	DO i = its , MIN(ide-1,ite)
2179	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2180	tslb(i,1,j)= tsk(i,j)
2181	tslb(i,2,j)= tmn(i,j)
2182	END DO
2183	END DO
2184
2185	IF ( flag_sst .EQ. 1 ) THEN
2186	DO j = jts , MIN(jde-1,jte)
2187	DO i = its , MIN(ide-1,ite)
2188	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2189	IF ( landmask(i,j) .LT. 0.5 ) THEN
2190	DO l = 1 , num_soil_layers
2191	tslb(i,l,j)= sst(i,j)
2192	smois(i,l,j)= 1.0
2193	sh2o (i,l,j)= 1.0
2194	END DO
2195	END IF
2196	END DO
2197	END DO
2198	ELSE
2199	DO j = jts , MIN(jde-1,jte)
2200	DO i = its , MIN(ide-1,ite)
2201	IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
2202	IF ( landmask(i,j) .LT. 0.5 ) THEN
2203	DO l = 1 , num_soil_layers
2204	tslb(i,l,j)= tsk(i,j)
2205	smois(i,l,j)= 1.0
2206	sh2o (i,l,j)= 1.0
2207	END DO
2208	END IF
2209	END DO
2210	END DO
2211	END IF
2212
2213	DEALLOCATE (zhave)
2214
2215	END IF
2216
2217	END SUBROUTINE init_soil_7_real
2218
2219	SUBROUTINE aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu )
2220	IMPLICIT NONE
2221
2222	INTEGER , INTENT(IN) :: iswater
2223	INTEGER , INTENT(IN) :: num_veg_cat
2224	CHARACTER (LEN=4) , INTENT(IN) :: mminlu
2225	REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
2226
2227	! Local variables
2228
2229	INTEGER , PARAMETER :: num_special_bins = 3 ! grass, shrubs, trees; add a new line in the "cib" array if updating this
2230	INTEGER , PARAMETER :: max_cats_per_bin = 8 ! larger than max num of cats inside each of the special bins
2231	! So, no more than 8 total tree categories that we need to consider.
2232	INTEGER , PARAMETER :: fl = -1 ! Flag value so that we know this is not a valid category to consider.
2233	INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_usgs = &
2234	(/ 2 , 3 , 4 , 5 , 6 , 7 , 10 , fl , & ! grass
2235	8 , 9 , fl , fl , fl , fl , fl , fl , & ! shrubs
2236	11 , 12 , 13 , 14 , 15 , fl , fl , fl /) ! trees
2237	INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_modis = &
2238	(/ 9 , 10 , 12 , 14 , fl , fl , fl , fl , & ! grass
2239	6 , 7 , 8 , fl , fl , fl , fl , fl , & ! shrubs
2240	1 , 2 , 3 , 4 , 5 , fl , fl , fl /) ! trees
2241
2242	IF ( mminlu(1:4) .EQ. 'USGS' ) THEN
2243	CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2244	num_special_bins , max_cats_per_bin , fl , cib_usgs )
2245	ELSE IF ( mminlu(1:4) .EQ. 'MODI' ) THEN
2246	CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2247	num_special_bins , max_cats_per_bin , fl , cib_modis )
2248	END IF
2249
2250	END SUBROUTINE aggregate_categories_part1
2251
2252	SUBROUTINE aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
2253	num_special_bins , max_cats_per_bin , fl , cib )
2254	IMPLICIT NONE
2255
2256	INTEGER , INTENT(IN) :: iswater
2257	INTEGER , INTENT(IN) :: num_veg_cat
2258	REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
2259	INTEGER , INTENT(IN) :: num_special_bins , max_cats_per_bin , fl
2260	INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) , INTENT(IN) :: cib
2261
2262	! Local variables
2263
2264	REAL , DIMENSION(1:num_veg_cat) :: landuse_frac_work ! copy of the input array, allows us to be wreckless
2265	INTEGER , DIMENSION ( max_cats_per_bin , num_special_bins ) :: cats_in_bin
2266	INTEGER :: ib , ic ! indexes for the bin and the category loops
2267	REAL , DIMENSION ( num_special_bins ) :: bin_max_val , bin_sum ! max category value in this bin, sum of all cats in this bin
2268	INTEGER , DIMENSION ( num_special_bins ) :: bin_max_idx ! index of this maximum category in this bin
2269	INTEGER :: bin_work , bin_orig ! the bin from whence the maximum hails, respectively for the aggregated and the original data
2270	INTEGER :: max_cat_orig , max_cat_work
2271	REAL :: max_val_orig , max_val_work
2272
2273	! Find the max in the original. If it is >= 50%, no need to even be in here.
2274
2275	DO ic = 1 , num_veg_cat
2276	IF ( landuse_frac(ic) .GE. 0.5 ) THEN
2277	RETURN
2278	END IF
2279	END DO
2280
2281	! Put the categories in the bin into a 2d array.
2282
2283	cats_in_bin = RESHAPE ( cib , (/ max_cats_per_bin , num_special_bins /) )
2284
2285	! Make a copy of the incoming array so that we can eventually diff our working copy with
2286	! the original.
2287
2288	landuse_frac_work = landuse_frac
2289
2290	! Loop over each of the special bins that we know about. Find the max values and the locations of such.
2291
2292	DO ib = 1 , num_special_bins
2293
2294	! For this bin, we know about specific categories. We get the sum of all of those
2295	! categories that we know about for this bin, and we keep track of which category
2296	! is dominant.
2297
2298	bin_sum (ib) = 0
2299	bin_max_val(ib) = fl
2300
2301	cat_loop_accumulate : DO ic = 1 , max_cats_per_bin
2302
2303	! Have we run out of valid categories in this bin? For example, we would not necessarily have
2304	! the same number of tree categories as there are grass categories.
2305
2306	IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
2307	EXIT cat_loop_accumulate
2308	END IF
2309
2310	bin_sum(ib) = bin_sum(ib) + landuse_frac(cats_in_bin(ic,ib))
2311
2312	IF ( landuse_frac(cats_in_bin(ic,ib)) .GT. bin_max_val(ib) ) THEN
2313	bin_max_val(ib) = landuse_frac(cats_in_bin(ic,ib))
2314	bin_max_idx(ib) = cats_in_bin(ic,ib)
2315	END IF
2316
2317	END DO cat_loop_accumulate
2318
2319	cat_loop_assign : DO ic = 1 , max_cats_per_bin
2320
2321	! Plow through each cat in the bin. If we find the dominant one, he gets the total sum. If we land on
2322	! the other guys in the bin, they get set back to zero to maintain the original aggregate influence.
2323
2324	IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
2325	EXIT cat_loop_assign
2326	ELSE IF ( cats_in_bin(ic,ib) .EQ. bin_max_idx(ib) ) THEN
2327	landuse_frac_work(cats_in_bin(ic,ib)) = bin_sum(ib)
2328	ELSE
2329	landuse_frac_work(cats_in_bin(ic,ib)) = 0
2330	END IF
2331
2332	END DO cat_loop_assign
2333
2334	END DO
2335
2336	! Now we loop through the categorical data, and get the max+location for the original input data and the
2337	! modified work data. Water is not allowed to be a "max" category unless it is greater than 50%. We hit
2338	! that test up at the top already, so we can toss out water willy nilly here.
2339
2340	max_cat_orig = fl
2341	max_val_orig = 0
2342	max_cat_work = fl
2343	max_val_work = 0
2344
2345	DO ic = 1 , num_veg_cat
2346	IF ( ic .EQ. iswater ) THEN
2347	CYCLE
2348	END IF
2349	IF ( landuse_frac(ic) .GT. max_val_orig ) THEN
2350	max_val_orig = landuse_frac(ic)
2351	max_cat_orig = ic
2352	END IF
2353	IF ( landuse_frac_work(ic) .GT. max_val_work ) THEN
2354	max_val_work = landuse_frac_work(ic)
2355	max_cat_work = ic
2356	END IF
2357	END DO
2358
2359	! Find the bin for the maximimum value of the original data.
2360
2361	bin_orig = -1
2362	bin_loop_orig : DO ib = 1 , num_special_bins
2363	cat_loop_orig : DO ic = 1 , max_cats_per_bin
2364	IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
2365	EXIT cat_loop_orig
2366	ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_orig ) THEN
2367	bin_orig = ib
2368	EXIT bin_loop_orig
2369	END IF
2370	END DO cat_loop_orig
2371	END DO bin_loop_orig
2372
2373	! Find the bin for the maximimum value of the aggregated data.
2374
2375	bin_work = -1
2376	bin_loop_work : DO ib = 1 , num_special_bins
2377	cat_loop_work : DO ic = 1 , max_cats_per_bin
2378	IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
2379	EXIT cat_loop_work
2380	ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_work ) THEN
2381	bin_work = ib
2382	EXIT bin_loop_work
2383	END IF
2384	END DO cat_loop_work
2385	END DO bin_loop_work
2386
2387	! So the big question is, did the aggregation change the bin? If the aggregation does not change the resulting
2388	! bin, then we leave everything alone. However, if there would be a change in the eventual bin chosen by
2389	! the simple dominant-category method, then we become pro-active interventionists and rectify this heinous
2390	! injustice foisted upon the lowly flora.
2391
2392	IF ( bin_work .EQ. bin_orig ) THEN
2393	! No op, we do nothing.
2394	ELSE
2395	DO ic = 1 , max_cats_per_bin
2396	landuse_frac(cats_in_bin(ic,bin_work)) = landuse_frac_work(cats_in_bin(ic,bin_work))
2397	END DO
2398	END IF
2399
2400	END SUBROUTINE aggregate_categories_part2
2401
2402	END MODULE module_soil_pre
2403
2404	FUNCTION skip_middle_points_t ( ids , ide , jds , jde , &
2405	i_in , j_in , width , &
2406	subtleties_exist ) &
2407	RESULT ( skip_it )
2408
2409	IMPLICIT NONE
2410
2411	INTEGER , INTENT(IN) :: ids , ide , jds , jde
2412	INTEGER , INTENT(IN) :: i_in , j_in , width
2413	LOGICAL , INTENT(IN) :: subtleties_exist
2414
2415	LOGICAL :: skip_it
2416
2417	INTEGER , PARAMETER :: slop = 0
2418
2419	IF ( subtleties_exist ) THEN
2420	skip_it = .FALSE.
2421	ELSE
2422	IF ( ( i_in .GE. ids+width+slop ) .AND. ( i_in .LE. ide-1-width-slop ) .AND. &
2423	( j_in .GE. jds+width+slop ) .AND. ( j_in .LE. jde-1-width-slop ) ) THEN
2424	skip_it = .TRUE.
2425	ELSE
2426	skip_it = .FALSE.
2427	END IF
2428	END IF
2429
2430	END FUNCTION skip_middle_points_t
2431
2432	#else
2433
2434	MODULE module_soil_pre
2435
2436	USE module_date_time
2437	USE module_state_description
2438
2439	CONTAINS
2440
2441	SUBROUTINE process_percent_cat_new ( landmask , &
2442	landuse_frac , soil_top_cat , soil_bot_cat , &
2443	isltyp , ivgtyp , &
2444	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
2445	ids , ide , jds , jde , kds , kde , &
2446	ims , ime , jms , jme , kms , kme , &
2447	its , ite , jts , jte , kts , kte , &
2448	iswater )
2449
2450	IMPLICIT NONE
2451
2452	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2453	ims , ime , jms , jme , kms , kme , &
2454	its , ite , jts , jte , kts , kte , &
2455	iswater
2456	INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
2457	REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
2458	REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
2459	REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
2460	INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
2461	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
2462
2463	INTEGER :: i , j , l , ll, dominant_index
2464	REAL :: dominant_value
2465
2466	#ifdef WRF_CHEM
2467	! REAL :: lwthresh = .99
2468	REAL :: lwthresh = .50
2469	#else
2470	REAL :: lwthresh = .50
2471	#endif
2472
2473	INTEGER , PARAMETER :: iswater_soil = 14
2474	INTEGER :: iforce
2475	CHARACTER (LEN=1024) :: message
2476
2477	! Sanity check on the 50/50 points
2478
2479	DO j = jts , MIN(jde-1,jte)
2480	DO i = its , MIN(ide-1,ite)
2481	dominant_value = landuse_frac(i,iswater,j)
2482	IF ( dominant_value .EQ. lwthresh ) THEN
2483	DO l = 1 , num_veg_cat
2484	IF ( l .EQ. iswater ) CYCLE
2485	IF ( landuse_frac(i,l,j) .EQ. lwthresh ) THEN
2486	write(message, FMT='(I4,I4,A,I4,A,F12.4)') i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
2487	call wrf_message ( message )
2488	landuse_frac(i,l,j) = lwthresh - .01
2489	landuse_frac(i,iswater,j) = 1. - (lwthresh + 0.01)
2490	END IF
2491	END DO
2492	END IF
2493	END DO
2494	END DO
2495
2496	! Compute the dominant VEGETATION INDEX.
2497
2498	DO j = jts , MIN(jde-1,jte)
2499	DO i = its , MIN(ide-1,ite)
2500	dominant_value = landuse_frac(i,1,j)
2501	dominant_index = 1
2502	DO l = 2 , num_veg_cat
2503	IF ( l .EQ. iswater ) THEN
2504	! wait a bit
2505	ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
2506	dominant_value = landuse_frac(i,l,j)
2507	dominant_index = l
2508	END IF
2509	END DO
2510	IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
2511	dominant_value = landuse_frac(i,iswater,j)
2512	dominant_index = iswater
2513	#if 0
2514	si needs to beef up consistency checks before we can use this part
2515	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
2516	( dominant_value .EQ. lwthresh) ) THEN
2517	! no op
2518	#else
2519	ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.LT.iswater))THEN
2520	write(message,*)'temporary SI landmask fix'
2521	call wrf_message(trim(message))
2522	! no op
2523	ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.GT.iswater))THEN
2524	write(message,*)'temporary SI landmask fix'
2525	call wrf_message(trim(message))
2526	dominant_value = landuse_frac(i,iswater,j)
2527	dominant_index = iswater
2528	#endif
2529	ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
2530	( dominant_value .LT. lwthresh ) ) THEN
2531	dominant_value = landuse_frac(i,iswater,j)
2532	dominant_index = iswater
2533	END IF
2534	IF ( dominant_index .EQ. iswater ) THEN
2535	if(landmask(i,j).gt.lwthresh) then
2536	write(message,*) 'changing to water at point ',i,j
2537	call wrf_message(trim(message))
2538	write(message,) nint(landuse_frac(i,:,j)100)
2539	call wrf_message(trim(message))
2540	endif
2541	landmask(i,j) = 0
2542	ELSE IF ( dominant_index .NE. iswater ) THEN
2543	if(landmask(i,j).lt.lwthresh) then
2544	write(message,*) 'changing to land at point ',i,j
2545	call wrf_message(trim(message))
2546	write(message,) nint(landuse_frac(i,:,j)100)
2547	call wrf_message(trim(message))
2548	endif
2549	landmask(i,j) = 1
2550	END IF
2551	ivgtyp(i,j) = dominant_index
2552	END DO
2553	END DO
2554
2555	! Compute the dominant SOIL TEXTURE INDEX, TOP.
2556
2557	iforce = 0
2558	DO i = its , MIN(ide-1,ite)
2559	DO j = jts , MIN(jde-1,jte)
2560	dominant_value = soil_top_cat(i,1,j)
2561	dominant_index = 1
2562	IF ( landmask(i,j) .GT. lwthresh ) THEN
2563	DO l = 2 , num_soil_top_cat
2564	IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
2565	dominant_value = soil_top_cat(i,l,j)
2566	dominant_index = l
2567	END IF
2568	END DO
2569	IF ( dominant_value .LT. 0.01 ) THEN
2570	iforce = iforce + 1
2571	WRITE ( message , FMT = '(A,I4,I4)' ) &
2572	'based on landuse, changing soil to land at point ',i,j
2573	CALL wrf_debug(1,message)
2574	!atec WRITE ( message , FMT = '(16(i3,1x))' ) &
2575	!atec 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
2576	!atec CALL wrf_debug(1,message)
2577	!atec WRITE ( message , FMT = '(16(i3,1x))' ) &
2578	!atec nint(soil_top_cat(i,:,j)*100)
2579	!atec CALL wrf_debug(1,message)
2580	dominant_index = 8
2581	END IF
2582	ELSE
2583	dominant_index = iswater_soil
2584	END IF
2585	isltyp(i,j) = dominant_index
2586	END DO
2587	END DO
2588
2589	if(iforce.ne.0)then
2590	WRITE(message,FMT='(A,I4,A,I6)' ) &
2591	'forcing artificial silty clay loam at ',iforce,' points, out of ',&
2592	(MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
2593	CALL wrf_debug(0,message)
2594	endif
2595
2596	END SUBROUTINE process_percent_cat_new
2597
2598	SUBROUTINE process_soil_real ( tsk , tmn , &
2599	landmask , sst , &
2600	st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
2601	zs , dzs , tslb , smois , sh2o , &
2602	flag_sst , flag_soilt000, flag_soilm000, &
2603	ids , ide , jds , jde , kds , kde , &
2604	ims , ime , jms , jme , kms , kme , &
2605	its , ite , jts , jte , kts , kte , &
2606	sf_surface_physics , num_soil_layers , real_data_init_type , &
2607	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2608	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
2609
2610	IMPLICIT NONE
2611
2612	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2613	ims , ime , jms , jme , kms , kme , &
2614	its , ite , jts , jte , kts , kte , &
2615	sf_surface_physics , num_soil_layers , real_data_init_type , &
2616	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2617	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
2618
2619	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2620
2621	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2622
2623	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2624	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2625	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2626	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2627	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2628	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2629
2630	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2631	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2632	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2633	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2634
2635	INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
2636	REAL :: dominant_value
2637
2638	! Initialize the soil depth, and the soil temperature and moisture.
2639
2640	IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2641	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2642	CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
2643	landmask , sst , flag_sst , &
2644	ids , ide , jds , jde , kds , kde , &
2645	ims , ime , jms , jme , kms , kme , &
2646	its , ite , jts , jte , kts , kte )
2647
2648	ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME .or. sf_surface_physics .eq. 88 ) .AND. &
2649	( num_soil_layers .GT. 1 ) ) THEN
2650	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2651	CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2652	st_input , sm_input , sw_input , landmask , sst , &
2653	zs , dzs , &
2654	st_levels_input , sm_levels_input , sw_levels_input , &
2655	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2656	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2657	flag_sst , flag_soilt000 , flag_soilm000 , &
2658	ids , ide , jds , jde , kds , kde , &
2659	ims , ime , jms , jme , kms , kme , &
2660	its , ite , jts , jte , kts , kte )
2661	! CALL init_soil_old ( tsk , tmn , &
2662	! smois , tslb , zs , dzs , num_soil_layers , &
2663	! st000010_input , st010040_input , st040100_input , st100200_input , &
2664	! st010200_input , &
2665	! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
2666	! sm010200_input , &
2667	! landmask_input , sst_input , &
2668	! ids , ide , jds , jde , kds , kde , &
2669	! ims , ime , jms , jme , kms , kme , &
2670	! its , ite , jts , jte , kts , kte )
2671	ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2672	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2673	CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
2674	st_input , sm_input , landmask , sst , &
2675	zs , dzs , &
2676	st_levels_input , sm_levels_input , &
2677	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
2678	num_st_levels_alloc , num_sm_levels_alloc , &
2679	flag_sst , flag_soilt000 , flag_soilm000 , &
2680	ids , ide , jds , jde , kds , kde , &
2681	ims , ime , jms , jme , kms , kme , &
2682	its , ite , jts , jte , kts , kte )
2683	END IF
2684
2685	END SUBROUTINE process_soil_real
2686
2687	SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
2688	ivgtyp,isltyp,tslb,smois, &
2689	tsk,tmn,zs,dzs, &
2690	num_soil_layers, &
2691	sf_surface_physics , &
2692	ids,ide, jds,jde, kds,kde,&
2693	ims,ime, jms,jme, kms,kme,&
2694	its,ite, jts,jte, kts,kte )
2695
2696	IMPLICIT NONE
2697
2698	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
2699	ims,ime, jms,jme, kms,kme, &
2700	its,ite, jts,jte, kts,kte
2701
2702	INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
2703
2704	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
2705
2706	REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
2707
2708	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
2709	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
2710	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
2711
2712	! Local variables.
2713
2714	INTEGER :: itf,jtf
2715
2716	itf=MIN(ite,ide-1)
2717	jtf=MIN(jte,jde-1)
2718
2719	IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2720	CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2721	CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
2722	ivgtyp,zs,dzs,num_soil_layers, &
2723	ids,ide, jds,jde, kds,kde, &
2724	ims,ime, jms,jme, kms,kme, &
2725	its,ite, jts,jte, kts,kte )
2726	ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2727	CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2728	CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
2729	ivgtyp,isltyp,tslb,smois,tmn, &
2730	num_soil_layers, &
2731	ids,ide, jds,jde, kds,kde, &
2732	ims,ime, jms,jme, kms,kme, &
2733	its,ite, jts,jte, kts,kte )
2734	ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2735	CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2736
2737	END IF
2738
2739	END SUBROUTINE process_soil_ideal
2740
2741	SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
2742	tsk , ter , toposoil , landmask , flag_toposoil , &
2743	st000010 , st010040 , st040100 , st100200 , st010200 , &
2744	flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
2745	soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
2746	flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
2747	ids , ide , jds , jde , kds , kde , &
2748	ims , ime , jms , jme , kms , kme , &
2749	its , ite , jts , jte , kts , kte )
2750
2751	IMPLICIT NONE
2752
2753	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2754	ims , ime , jms , jme , kms , kme , &
2755	its , ite , jts , jte , kts , kte
2756
2757	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
2758	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
2759	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
2760	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
2761
2762	INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
2763	INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
2764	INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
2765
2766	INTEGER :: i , j
2767
2768	REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
2769
2770	! Do we have a soil field with which to modify soil temperatures?
2771
2772	IF ( flag_toposoil .EQ. 1 ) THEN
2773
2774	DO j = jts , MIN(jde-1,jte)
2775	DO i = its , MIN(ide-1,ite)
2776
2777	! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
2778	! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
2779	! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
2780	! the difference between the soil elevation and the terrain is greater than 3 km means
2781	! that the soil data is either all zeros or that the data are inconsistent. Any of these
2782	! three conditions is grievous enough to induce a WRF fatality. However, if they are at
2783	! a water point, then we can safely ignore them.
2784
2785	soil_elev_min_val = toposoil(i,j)
2786	soil_elev_max_val = toposoil(i,j)
2787	soil_elev_min_dif = ter(i,j) - toposoil(i,j)
2788	soil_elev_max_dif = ter(i,j) - toposoil(i,j)
2789
2790	IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2791	CYCLE
2792	ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2793	!print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
2794	cycle
2795	! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
2796	ENDIF
2797
2798	IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2799	CYCLE
2800	ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2801	print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
2802	cycle
2803	CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
2804	ENDIF
2805
2806	IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2807	( landmask(i,j) .LT. 0.5 ) ) THEN
2808	CYCLE
2809	ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2810	( landmask(i,j) .GT. 0.5 ) ) THEN
2811	print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
2812	cycle
2813	CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
2814	ENDIF
2815
2816	! For each of the fields that we would like to modify, check to see if it came in from the SI.
2817	! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
2818	! are not at a water point.
2819
2820	IF (landmask(i,j) .GT. 0.5 ) THEN
2821
2822	IF ( sf_surface_physics .EQ. SLABSCHEME .OR. sf_surface_physics .EQ. PXLSMSCHEME ) THEN
2823	tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2824	END IF
2825
2826	tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2827
2828	IF ( flag_st000010 .EQ. 1 ) THEN
2829	st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2830	END IF
2831	IF ( flag_st010040 .EQ. 1 ) THEN
2832	st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2833	END IF
2834	IF ( flag_st040100 .EQ. 1 ) THEN
2835	st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2836	END IF
2837	IF ( flag_st100200 .EQ. 1 ) THEN
2838	st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2839	END IF
2840	IF ( flag_st010200 .EQ. 1 ) THEN
2841	st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2842	END IF
2843
2844	IF ( flag_soilt000 .EQ. 1 ) THEN
2845	soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2846	END IF
2847	IF ( flag_soilt005 .EQ. 1 ) THEN
2848	soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2849	END IF
2850	IF ( flag_soilt020 .EQ. 1 ) THEN
2851	soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2852	END IF
2853	IF ( flag_soilt040 .EQ. 1 ) THEN
2854	soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2855	END IF
2856	IF ( flag_soilt160 .EQ. 1 ) THEN
2857	soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2858	END IF
2859	IF ( flag_soilt300 .EQ. 1 ) THEN
2860	soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2861	END IF
2862
2863	END IF
2864	END DO
2865	END DO
2866
2867	END IF
2868
2869	END SUBROUTINE adjust_soil_temp_new
2870
2871
2872	SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
2873
2874	IMPLICIT NONE
2875
2876	INTEGER, INTENT(IN) :: num_soil_layers
2877
2878	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2879
2880	INTEGER :: l
2881
2882	CHARACTER (LEN=132) :: message
2883
2884	! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
2885	! The distance from the ground level to the midpoint of the layer is given by zs.
2886
2887	! ------- Ground Level ---------- \|\| \|\| \|\| \|\|
2888	! \|\| \|\| \|\| \|\| zs(1) = 0.005 m
2889	! -- -- -- -- -- -- -- -- -- \|\| \|\| \|\| \/
2890	! \|\| \|\| \|\|
2891	! ----------------------------------- \|\| \|\| \|\| \/ dzs(1) = 0.01 m
2892	! \|\| \|\| \|\|
2893	! \|\| \|\| \|\| zs(2) = 0.02
2894	! -- -- -- -- -- -- -- -- -- \|\| \|\| \/
2895	! \|\| \|\|
2896	! \|\| \|\|
2897	! ----------------------------------- \|\| \|\| \/ dzs(2) = 0.02 m
2898	! \|\| \|\|
2899	! \|\| \|\|
2900	! \|\| \|\|
2901	! \|\| \|\| zs(3) = 0.05
2902	! -- -- -- -- -- -- -- -- -- \|\| \/
2903	! \|\|
2904	! \|\|
2905	! \|\|
2906	! \|\|
2907	! ----------------------------------- \/ dzs(3) = 0.04 m
2908
2909	IF ( num_soil_layers .NE. 5 ) THEN
2910	write(message,FMT= '(A)') 'Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
2911	CALL wrf_error_fatal ( message )
2912	END IF
2913
2914	dzs(1)=.01
2915	zs(1)=.5*dzs(1)
2916
2917	DO l=2,num_soil_layers
2918	dzs(l)=2*dzs(l-1)
2919	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
2920	ENDDO
2921
2922	END SUBROUTINE init_soil_depth_1
2923
2924	SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
2925
2926	IMPLICIT NONE
2927
2928	INTEGER, INTENT(IN) :: num_soil_layers
2929
2930	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2931
2932	INTEGER :: l
2933
2934	CHARACTER (LEN=132) :: message
2935
2936	dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
2937
2938	IF ( num_soil_layers .NE. 4 ) THEN
2939	write(message,FMT='(A)') 'Usually, the LSM uses 4 layers. Change this in the namelist.'
2940	CALL wrf_error_fatal ( message )
2941	END IF
2942
2943	zs(1)=.5*dzs(1)
2944
2945	DO l=2,num_soil_layers
2946	zs(l)=zs(l-1)+.5dzs(l-1)+.5dzs(l)
2947	ENDDO
2948
2949	END SUBROUTINE init_soil_depth_2
2950
2951	SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
2952
2953	IMPLICIT NONE
2954
2955	INTEGER, INTENT(IN) :: num_soil_layers
2956
2957	REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
2958
2959	INTEGER :: l
2960
2961	CHARACTER (LEN=132) :: message
2962
2963	! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
2964	! ZS is specified in the namelist: num_soil_layers = 6 or 9.
2965	! Other options with number of levels are possible, but
2966	! WRF users should change consistently the namelist entry with the
2967	! ZS array in this subroutine.
2968
2969	IF ( num_soil_layers .EQ. 6) THEN
2970	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
2971	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2972	ELSEIF ( num_soil_layers .EQ. 9) THEN
2973	zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
2974	! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2975	ENDIF
2976
2977	IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
2978	WRITE(message,FMT= '(A)')'Usually, the RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
2979	CALL wrf_error_fatal ( message )
2980	END IF
2981
2982	END SUBROUTINE init_soil_depth_3
2983
2984	SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
2985	num_soil_layers , real_data_init_type , &
2986	landmask , sst , flag_sst , &
2987	ids , ide , jds , jde , kds , kde , &
2988	ims , ime , jms , jme , kms , kme , &
2989	its , ite , jts , jte , kts , kte )
2990
2991	IMPLICIT NONE
2992
2993	INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
2994	ids , ide , jds , jde , kds , kde , &
2995	ims , ime , jms , jme , kms , kme , &
2996	its , ite , jts , jte , kts , kte
2997
2998	INTEGER , INTENT(IN) :: flag_sst
2999
3000	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
3001	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
3002
3003	REAL , DIMENSION(num_soil_layers) :: zs , dzs
3004
3005	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
3006
3007	INTEGER :: i , j , l
3008
3009	! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
3010	! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
3011	! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
3012	! annual mean temperature.
3013
3014	DO j = jts , MIN(jde-1,jte)
3015	DO i = its , MIN(ide-1,ite)
3016	IF ( landmask(i,j) .GT. 0.5 ) THEN
3017	DO l = 1 , num_soil_layers
3018	tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
3019	tmn(i,j) * ( zs( l) - zs(1) ) ) / &
3020	( zs(num_soil_layers) - zs(1) )
3021	END DO
3022	ELSE
3023	IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
3024	DO l = 1 , num_soil_layers
3025	tslb(i,l,j)= sst(i,j)
3026	END DO
3027	ELSE
3028	DO l = 1 , num_soil_layers
3029	tslb(i,l,j)= tsk(i,j)
3030	END DO
3031	END IF
3032	END IF
3033	END DO
3034	END DO
3035
3036	END SUBROUTINE init_soil_1_real
3037
3038	SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
3039	ivgtyp,ZS,DZS,num_soil_layers, &
3040	ids,ide, jds,jde, kds,kde, &
3041	ims,ime, jms,jme, kms,kme, &
3042	its,ite, jts,jte, kts,kte )
3043
3044	IMPLICIT NONE
3045
3046	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
3047	ims,ime, jms,jme, kms,kme, &
3048	its,ite, jts,jte, kts,kte
3049
3050	INTEGER, INTENT(IN ) :: num_soil_layers
3051
3052	REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
3053	REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
3054	INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
3055
3056	REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
3057
3058	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
3059
3060	! Lcal variables.
3061
3062	INTEGER :: l,j,i,itf,jtf
3063
3064	itf=MIN(ite,ide-1)
3065	jtf=MIN(jte,jde-1)
3066
3067	IF (num_soil_layers.NE.1)THEN
3068	DO j=jts,jtf
3069	DO l=1,num_soil_layers
3070	DO i=its,itf
3071	tslb(i,l,j)=( tsk(i,j)(zs(num_soil_layers)-zs(l)) + tmn(i,j)(zs(l)-zs(1)) ) / &
3072	( zs(num_soil_layers)-zs(1) )
3073	ENDDO
3074	ENDDO
3075	ENDDO
3076	ENDIF
3077	DO j=jts,jtf
3078	DO i=its,itf
3079	xland(i,j) = 2
3080	ivgtyp(i,j) = 7
3081	ENDDO
3082	ENDDO
3083
3084	END SUBROUTINE init_soil_1_ideal
3085
3086	SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
3087	st_input , sm_input , sw_input , landmask , sst , &
3088	zs , dzs , &
3089	st_levels_input , sm_levels_input , sw_levels_input , &
3090	num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
3091	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
3092	flag_sst , flag_soilt000 , flag_soilm000 , &
3093	ids , ide , jds , jde , kds , kde , &
3094	ims , ime , jms , jme , kms , kme , &
3095	its , ite , jts , jte , kts , kte )
3096
3097	IMPLICIT NONE
3098
3099	INTEGER , INTENT(IN) :: num_soil_layers , &
3100	num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
3101	num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
3102	ids , ide , jds , jde , kds , kde , &
3103	ims , ime , jms , jme , kms , kme , &
3104	its , ite , jts , jte , kts , kte
3105
3106	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
3107
3108	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
3109	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
3110	INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
3111
3112	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
3113	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
3114	REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
3115	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
3116
3117	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
3118	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
3119	REAL , DIMENSION(num_soil_layers) :: zs , dzs
3120
3121	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
3122
3123	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
3124
3125	INTEGER :: i , j , l , lout , lin , lwant , lhave , num
3126	REAL :: temp
3127	LOGICAL :: found_levels
3128
3129	CHARACTER (LEN=132) :: message
3130
3131	! Are there any soil temp and moisture levels - ya know, they are mandatory.
3132
3133	num = num_st_levels_input * num_sm_levels_input
3134
3135	IF ( num .GE. 1 ) THEN
3136
3137	! Ordered levels that we have data for.
3138	IF ( flag_soilt000 .eq. 1 ) THEN
3139	write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
3140	CALL wrf_message ( message )
3141	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
3142	ELSE
3143	write(message, FMT='(A)') ' Assume Noah LSM input'
3144	CALL wrf_message ( message )
3145	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
3146	END IF
3147
3148
3149	! Sort the levels for temperature.
3150	!print*,'num_st_levels_input, st_levels_input', num_st_levels_input, st_levels_input
3151	!print*,'num_sm_levels_input,num_sw_levels_input',num_sm_levels_input,num_sw_levels_input
3152
3153	outert : DO lout = 1 , num_st_levels_input-1
3154	innert : DO lin = lout+1 , num_st_levels_input
3155	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
3156	temp = st_levels_input(lout)
3157	st_levels_input(lout) = st_levels_input(lin)
3158	st_levels_input(lin) = NINT(temp)
3159	DO j = jts , MIN(jde-1,jte)
3160	DO i = its , MIN(ide-1,ite)
3161	temp = st_input(i,lout+1,j)
3162	st_input(i,lout+1,j) = st_input(i,lin+1,j)
3163	st_input(i,lin+1,j) = temp
3164	END DO
3165	END DO
3166	END IF
3167	END DO innert
3168	END DO outert
3169	!tgs add IF
3170	IF ( flag_soilt000 .NE. 1 ) THEN
3171	DO j = jts , MIN(jde-1,jte)
3172	DO i = its , MIN(ide-1,ite)
3173	st_input(i,1,j) = tsk(i,j)
3174	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
3175	END DO
3176	END DO
3177	ENDIF
3178
3179
3180	#if ( NMM_CORE == 1 )
3181	!new
3182	! write(0,*) 'st_input(1) in init_soil_2_real'
3183	! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
3184	! write(0,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
3185	! ENDDO
3186
3187	! write(0,*) 'st_input(2) in init_soil_2_real'
3188	! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
3189	! write(0,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
3190	! ENDDO
3191
3192	616 format(20(f4.0,1x))
3193	#endif
3194
3195	! Sort the levels for moisture.
3196
3197	outerm: DO lout = 1 , num_sm_levels_input-1
3198	innerm : DO lin = lout+1 , num_sm_levels_input
3199	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
3200	temp = sm_levels_input(lout)
3201	sm_levels_input(lout) = sm_levels_input(lin)
3202	sm_levels_input(lin) = NINT(temp)
3203	DO j = jts , MIN(jde-1,jte)
3204	DO i = its , MIN(ide-1,ite)
3205	temp = sm_input(i,lout+1,j)
3206	sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
3207	sm_input(i,lin+1,j) = temp
3208	END DO
3209	END DO
3210	END IF
3211	END DO innerm
3212	END DO outerm
3213	!tgs add IF
3214	IF ( flag_soilm000 .NE. 1 ) THEN
3215	DO j = jts , MIN(jde-1,jte)
3216	DO i = its , MIN(ide-1,ite)
3217	sm_input(i,1,j) = sm_input(i,2,j)
3218	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
3219	END DO
3220	END DO
3221	ENDIF
3222
3223	! Sort the levels for liquid moisture.
3224
3225	outerw: DO lout = 1 , num_sw_levels_input-1
3226	innerw : DO lin = lout+1 , num_sw_levels_input
3227	IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
3228	temp = sw_levels_input(lout)
3229	sw_levels_input(lout) = sw_levels_input(lin)
3230	sw_levels_input(lin) = NINT(temp)
3231	DO j = jts , MIN(jde-1,jte)
3232	DO i = its , MIN(ide-1,ite)
3233	temp = sw_input(i,lout+1,j)
3234	sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
3235	sw_input(i,lin+1,j) = temp
3236	END DO
3237	END DO
3238	END IF
3239	END DO innerw
3240	END DO outerw
3241	IF ( num_sw_levels_input .GT. 1 ) THEN
3242	DO j = jts , MIN(jde-1,jte)
3243	DO i = its , MIN(ide-1,ite)
3244	sw_input(i,1,j) = sw_input(i,2,j)
3245	sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
3246	END DO
3247	END DO
3248	END IF
3249
3250	found_levels = .TRUE.
3251
3252	ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
3253
3254	found_levels = .FALSE.
3255
3256	ELSE
3257	CALL wrf_error_fatal ( &
3258	'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
3259	END IF
3260
3261	! Is it OK to continue?
3262
3263	IF ( found_levels ) THEN
3264
3265	!tgs add IF
3266	IF ( flag_soilt000 .NE.1) THEN
3267	!tgs initialize from Noah data
3268	! Here are the levels that we have from the input for temperature. The input levels plus
3269	! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
3270
3271	zhave(1) = 0.
3272	DO l = 1 , num_st_levels_input
3273	zhave(l+1) = st_levels_input(l) / 100.
3274	END DO
3275	zhave(num_st_levels_input+2) = 300. / 100.
3276
3277	! Interpolate between the layers we have (zhave) and those that we want (zs).
3278
3279	z_wantt : DO lwant = 1 , num_soil_layers
3280	z_havet : DO lhave = 1 , num_st_levels_input +2 -1
3281	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3282	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3283	DO j = jts , MIN(jde-1,jte)
3284	DO i = its , MIN(ide-1,ite)
3285	tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
3286	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3287	( zhave(lhave+1) - zhave(lhave) )
3288	END DO
3289	END DO
3290	EXIT z_havet
3291	END IF
3292	END DO z_havet
3293	END DO z_wantt
3294
3295	ELSE
3296
3297	!tgs initialize from RUCLSM data
3298	DO l = 1 , num_st_levels_input
3299	zhave(l) = st_levels_input(l) / 100.
3300	END DO
3301
3302	! Interpolate between the layers we have (zhave) and those that we want (zs).
3303
3304
3305	z_wantt_2 : DO lwant = 1 , num_soil_layers
3306	z_havet_2 : DO lhave = 1 , num_st_levels_input -1
3307	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3308	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3309	DO j = jts , MIN(jde-1,jte)
3310	DO i = its , MIN(ide-1,ite)
3311	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3312	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3313	( zhave(lhave+1) - zhave(lhave) )
3314	END DO
3315	END DO
3316	EXIT z_havet_2
3317	END IF
3318	END DO z_havet_2
3319	END DO z_wantt_2
3320
3321	ENDIF
3322
3323
3324	IF ( flag_soilm000 .NE. 1 ) THEN
3325	!tgs initialize from Noah
3326	! Here are the levels that we have from the input for moisture. The input levels plus
3327	! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
3328	! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
3329	! as the most deep layer's value.
3330
3331	zhave(1) = 0.
3332	DO l = 1 , num_sm_levels_input
3333	zhave(l+1) = sm_levels_input(l) / 100.
3334	END DO
3335	zhave(num_sm_levels_input+2) = 300. / 100.
3336
3337	! Interpolate between the layers we have (zhave) and those that we want (zs).
3338
3339	z_wantm : DO lwant = 1 , num_soil_layers
3340	z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
3341	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3342	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3343	DO j = jts , MIN(jde-1,jte)
3344	DO i = its , MIN(ide-1,ite)
3345	smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
3346	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3347	( zhave(lhave+1) - zhave(lhave) )
3348	END DO
3349	END DO
3350	EXIT z_havem
3351	END IF
3352	END DO z_havem
3353	END DO z_wantm
3354
3355	ELSE
3356
3357	!tgs initialize from RUCLSM data
3358	DO l = 1 , num_sm_levels_input
3359	zhave(l) = sm_levels_input(l) / 100.
3360	END DO
3361
3362	! Interpolate between the layers we have (zhave) and those that we want (zs).
3363
3364	z_wantm_2 : DO lwant = 1 , num_soil_layers
3365	z_havem_2 : DO lhave = 1 , num_sm_levels_input -1
3366	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3367	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3368	DO j = jts , MIN(jde-1,jte)
3369	DO i = its , MIN(ide-1,ite)
3370	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3371	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3372	( zhave(lhave+1) - zhave(lhave) )
3373	END DO
3374	END DO
3375	EXIT z_havem_2
3376	END IF
3377	END DO z_havem_2
3378	END DO z_wantm_2
3379
3380	ENDIF
3381	! Any liquid soil moisture to worry about?
3382
3383	IF ( num_sw_levels_input .GT. 1 ) THEN
3384
3385	zhave(1) = 0.
3386	DO l = 1 , num_sw_levels_input
3387	zhave(l+1) = sw_levels_input(l) / 100.
3388	END DO
3389	zhave(num_sw_levels_input+2) = 300. / 100.
3390
3391	! Interpolate between the layers we have (zhave) and those that we want (zs).
3392
3393	z_wantw : DO lwant = 1 , num_soil_layers
3394	z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
3395	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3396	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3397	DO j = jts , MIN(jde-1,jte)
3398	DO i = its , MIN(ide-1,ite)
3399	sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
3400	sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3401	( zhave(lhave+1) - zhave(lhave) )
3402	END DO
3403	END DO
3404	EXIT z_havew
3405	END IF
3406	END DO z_havew
3407	END DO z_wantw
3408
3409	END IF
3410
3411
3412	! Over water, put in reasonable values for soil temperature and moisture. These won't be
3413	! used, but they will make a more continuous plot.
3414
3415	IF ( flag_sst .EQ. 1 ) THEN
3416	DO j = jts , MIN(jde-1,jte)
3417	DO i = its , MIN(ide-1,ite)
3418	IF ( landmask(i,j) .LT. 0.5 ) THEN
3419	DO l = 1 , num_soil_layers
3420	tslb(i,l,j)= sst(i,j)
3421	!tgs add line for tsk
3422	tsk(i,j) = sst(i,j)
3423	smois(i,l,j)= 1.0
3424	sh2o (i,l,j)= 1.0
3425	END DO
3426	END IF
3427	END DO
3428	END DO
3429	ELSE
3430	DO j = jts , MIN(jde-1,jte)
3431	DO i = its , MIN(ide-1,ite)
3432	IF ( landmask(i,j) .LT. 0.5 ) THEN
3433	DO l = 1 , num_soil_layers
3434	tslb(i,l,j)= tsk(i,j)
3435	smois(i,l,j)= 1.0
3436	sh2o (i,l,j)= 1.0
3437	END DO
3438	END IF
3439	END DO
3440	END DO
3441	END IF
3442
3443	DEALLOCATE (zhave)
3444
3445	END IF
3446
3447	END SUBROUTINE init_soil_2_real
3448
3449	SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
3450	ivgtyp,isltyp,tslb,smois,tmn, &
3451	num_soil_layers, &
3452	ids,ide, jds,jde, kds,kde, &
3453	ims,ime, jms,jme, kms,kme, &
3454	its,ite, jts,jte, kts,kte )
3455
3456	IMPLICIT NONE
3457
3458	INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
3459	ims,ime, jms,jme, kms,kme, &
3460	its,ite, jts,jte, kts,kte
3461
3462	INTEGER, INTENT(IN) ::num_soil_layers
3463
3464	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
3465
3466	REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
3467
3468	INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
3469
3470	INTEGER :: icm,jcm,itf,jtf
3471	INTEGER :: i,j,l
3472
3473	itf=min0(ite,ide-1)
3474	jtf=min0(jte,jde-1)
3475
3476	icm = ide/2
3477	jcm = jde/2
3478
3479	DO j=jts,jtf
3480	DO l=1,num_soil_layers
3481	DO i=its,itf
3482
3483	smois(i,1,j)=0.10
3484	smois(i,2,j)=0.10
3485	smois(i,3,j)=0.10
3486	smois(i,4,j)=0.10
3487
3488	tslb(i,1,j)=295.
3489	tslb(i,2,j)=297.
3490	tslb(i,3,j)=293.
3491	tslb(i,4,j)=293.
3492
3493	ENDDO
3494	ENDDO
3495	ENDDO
3496
3497	DO j=jts,jtf
3498	DO i=its,itf
3499	xland(i,j) = 2
3500	tmn(i,j) = 294.
3501	xice(i,j) = 0.
3502	vegfra(i,j) = 0.
3503	snow(i,j) = 0.
3504	canwat(i,j) = 0.
3505	ivgtyp(i,j) = 7
3506	isltyp(i,j) = 8
3507	ENDDO
3508	ENDDO
3509
3510	END SUBROUTINE init_soil_2_ideal
3511
3512	SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
3513	st_input , sm_input , landmask, sst, &
3514	zs , dzs , &
3515	st_levels_input , sm_levels_input , &
3516	num_soil_layers , num_st_levels_input , num_sm_levels_input , &
3517	num_st_levels_alloc , num_sm_levels_alloc , &
3518	flag_sst , flag_soilt000 , flag_soilm000 , &
3519	ids , ide , jds , jde , kds , kde , &
3520	ims , ime , jms , jme , kms , kme , &
3521	its , ite , jts , jte , kts , kte )
3522
3523	IMPLICIT NONE
3524
3525	INTEGER , INTENT(IN) :: num_soil_layers , &
3526	num_st_levels_input , num_sm_levels_input , &
3527	num_st_levels_alloc , num_sm_levels_alloc , &
3528	ids , ide , jds , jde , kds , kde , &
3529	ims , ime , jms , jme , kms , kme , &
3530	its , ite , jts , jte , kts , kte
3531
3532	INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
3533
3534	INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
3535	INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
3536
3537	REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
3538	REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
3539	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
3540
3541	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
3542	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
3543	REAL , DIMENSION(num_soil_layers) :: zs , dzs
3544
3545	REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
3546
3547	REAL , ALLOCATABLE , DIMENSION(:) :: zhave
3548
3549	INTEGER :: i , j , l , lout , lin , lwant , lhave
3550	REAL :: temp
3551
3552	! Allocate the soil layer array used for interpolating.
3553
3554	IF ( ( num_st_levels_input .LE. 0 ) .OR. &
3555	( num_sm_levels_input .LE. 0 ) ) THEN
3556	PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
3557	CALL wrf_error_fatal ( 'no soil data' )
3558	ELSE
3559	IF ( flag_soilt000 .eq. 1 ) THEN
3560	PRINT '(A)',' Assume RUC LSM 6-level input'
3561	ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
3562	ELSE
3563	PRINT '(A)',' Assume non-RUC LSM input'
3564	ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
3565	END IF
3566	END IF
3567
3568	! Sort the levels for temperature.
3569
3570	outert : DO lout = 1 , num_st_levels_input-1
3571	innert : DO lin = lout+1 , num_st_levels_input
3572	IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
3573	temp = st_levels_input(lout)
3574	st_levels_input(lout) = st_levels_input(lin)
3575	st_levels_input(lin) = NINT(temp)
3576	DO j = jts , MIN(jde-1,jte)
3577	DO i = its , MIN(ide-1,ite)
3578	temp = st_input(i,lout,j)
3579	st_input(i,lout,j) = st_input(i,lin,j)
3580	st_input(i,lin,j) = temp
3581	END DO
3582	END DO
3583	END IF
3584	END DO innert
3585	END DO outert
3586
3587	IF ( flag_soilt000 .NE. 1 ) THEN
3588	DO j = jts , MIN(jde-1,jte)
3589	DO i = its , MIN(ide-1,ite)
3590	st_input(i,1,j) = tsk(i,j)
3591	st_input(i,num_st_levels_input+2,j) = tmn(i,j)
3592	END DO
3593	END DO
3594	END IF
3595
3596	! Sort the levels for moisture.
3597
3598	outerm: DO lout = 1 , num_sm_levels_input-1
3599	innerm : DO lin = lout+1 , num_sm_levels_input
3600	IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
3601	temp = sm_levels_input(lout)
3602	sm_levels_input(lout) = sm_levels_input(lin)
3603	sm_levels_input(lin) = NINT(temp)
3604	DO j = jts , MIN(jde-1,jte)
3605	DO i = its , MIN(ide-1,ite)
3606	temp = sm_input(i,lout,j)
3607	sm_input(i,lout,j) = sm_input(i,lin,j)
3608	sm_input(i,lin,j) = temp
3609	END DO
3610	END DO
3611	END IF
3612	END DO innerm
3613	END DO outerm
3614
3615	IF ( flag_soilm000 .NE. 1 ) THEN
3616	DO j = jts , MIN(jde-1,jte)
3617	DO i = its , MIN(ide-1,ite)
3618	sm_input(i,1,j) = sm_input(i,2,j)
3619	sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
3620	END DO
3621	END DO
3622	END IF
3623
3624	! Here are the levels that we have from the input for temperature.
3625
3626	IF ( flag_soilt000 .EQ. 1 ) THEN
3627	DO l = 1 , num_st_levels_input
3628	zhave(l) = st_levels_input(l) / 100.
3629	END DO
3630
3631	! Interpolate between the layers we have (zhave) and those that we want (zs).
3632
3633	z_wantt : DO lwant = 1 , num_soil_layers
3634	z_havet : DO lhave = 1 , num_st_levels_input -1
3635	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3636	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3637	DO j = jts , MIN(jde-1,jte)
3638	DO i = its , MIN(ide-1,ite)
3639	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3640	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3641	( zhave(lhave+1) - zhave(lhave) )
3642	END DO
3643	END DO
3644	EXIT z_havet
3645	END IF
3646	END DO z_havet
3647	END DO z_wantt
3648
3649	ELSE
3650
3651	zhave(1) = 0.
3652	DO l = 1 , num_st_levels_input
3653	zhave(l+1) = st_levels_input(l) / 100.
3654	END DO
3655	zhave(num_st_levels_input+2) = 300. / 100.
3656
3657	! Interpolate between the layers we have (zhave) and those that we want (zs).
3658
3659	z_wantt_2 : DO lwant = 1 , num_soil_layers
3660	z_havet_2 : DO lhave = 1 , num_st_levels_input +2
3661	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3662	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3663	DO j = jts , MIN(jde-1,jte)
3664	DO i = its , MIN(ide-1,ite)
3665	tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3666	st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3667	( zhave(lhave+1) - zhave(lhave) )
3668	END DO
3669	END DO
3670	EXIT z_havet_2
3671	END IF
3672	END DO z_havet_2
3673	END DO z_wantt_2
3674
3675	END IF
3676
3677	! Here are the levels that we have from the input for moisture.
3678
3679	IF ( flag_soilm000 .EQ. 1 ) THEN
3680	DO l = 1 , num_sm_levels_input
3681	zhave(l) = sm_levels_input(l) / 100.
3682	END DO
3683
3684	! Interpolate between the layers we have (zhave) and those that we want (zs).
3685
3686	z_wantm : DO lwant = 1 , num_soil_layers
3687	z_havem : DO lhave = 1 , num_sm_levels_input -1
3688	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3689	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3690	DO j = jts , MIN(jde-1,jte)
3691	DO i = its , MIN(ide-1,ite)
3692	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3693	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3694	( zhave(lhave+1) - zhave(lhave) )
3695	END DO
3696	END DO
3697	EXIT z_havem
3698	END IF
3699	END DO z_havem
3700	END DO z_wantm
3701
3702	ELSE
3703
3704	zhave(1) = 0.
3705	DO l = 1 , num_sm_levels_input
3706	zhave(l+1) = sm_levels_input(l) / 100.
3707	END DO
3708	zhave(num_sm_levels_input+2) = 300. / 100.
3709
3710	z_wantm_2 : DO lwant = 1 , num_soil_layers
3711	z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
3712	IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
3713	( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3714	DO j = jts , MIN(jde-1,jte)
3715	DO i = its , MIN(ide-1,ite)
3716	smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
3717	sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
3718	( zhave(lhave+1) - zhave(lhave) )
3719	END DO
3720	END DO
3721	EXIT z_havem_2
3722	END IF
3723	END DO z_havem_2
3724	END DO z_wantm_2
3725
3726	END IF
3727	! Over water, put in reasonable values for soil temperature and moisture. These won't be
3728	! used, but they will make a more continuous plot.
3729
3730	IF ( flag_sst .EQ. 1 ) THEN
3731	DO j = jts , MIN(jde-1,jte)
3732	DO i = its , MIN(ide-1,ite)
3733	IF ( landmask(i,j) .LT. 0.5 ) THEN
3734	DO l = 1 , num_soil_layers
3735	tslb(i,l,j) = sst(i,j)
3736	tsk(i,j) = sst(i,j)
3737	smois(i,l,j)= 1.0
3738	END DO
3739	END IF
3740	END DO
3741	END DO
3742	ELSE
3743	DO j = jts , MIN(jde-1,jte)
3744	DO i = its , MIN(ide-1,ite)
3745	IF ( landmask(i,j) .LT. 0.5 ) THEN
3746	DO l = 1 , num_soil_layers
3747	tslb(i,l,j)= tsk(i,j)
3748	smois(i,l,j)= 1.0
3749	END DO
3750	END IF
3751	END DO
3752	END DO
3753	END IF
3754
3755	DEALLOCATE (zhave)
3756
3757	END SUBROUTINE init_soil_3_real
3758
3759	END MODULE module_soil_pre
3760
3761	#endif

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

Download in other formats:

Original Format