Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diff_module_initialize_real @ 2884

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

Line
1	9,11d8
2	< !****MARS: modified May 2007
3	<
4	<
5	224,225d220
6	< read (20+grid%id) grid%em_albedo_gcm
7	< read (20+grid%id) grid%em_therm_inert
8	243,244d237
9	< write (20+grid%id) grid%em_albedo_gcm
10	< write (20+grid%id) grid%em_therm_inert
11	258d250
12	< !!****MARS: tsk is surface temperature
13	264,268d255
14	< !!****MARS
15	< !!un peu artificiel, mais u10 et v10 sont des bons intermediaires (facultatifs de plus)
16	< grid%u10(i,j) = grid%em_albedo_gcm(i,j)
17	< grid%v10(i,j) = grid%em_therm_inert(i,j)
18	< !!****MARS
19	272d258
20	<
21	275,286d260
22	< !!****MARS
23	< !!fix pour être certain d'être avec les bons flag
24	< print *,flag_psfc
25	< flag_psfc=1
26	< print *,flag_soilhgt
27	< flag_soilhgt=1
28	< print *,flag_metgrid
29	< flag_metgrid=1
30	< !!**** TODO: trouver quand même pourquoi ça donne 0 :)
31	< !!****MARS
32	<
33	<
34	317,330c291,302
35	< !****MARS
36	< ! DO j = jts, min(jde-1,jte)
37	< ! DO i = its, min(ide,ite)
38	< ! grid%u10(i,j)=grid%em_u_gc(i,1,j)
39	< ! END DO
40	< ! END DO
41	< !
42	< ! DO j = jts, min(jde,jte)
43	< ! DO i = its, min(ide-1,ite)
44	< ! grid%v10(i,j)=grid%em_v_gc(i,1,j)
45	< ! END DO
46	< ! END DO
47	< !****MARS
48	<
49	---
50	> DO j = jts, min(jde-1,jte)
51	> DO i = its, min(ide,ite)
52	> grid%u10(i,j)=grid%em_u_gc(i,1,j)
53	> END DO
54	> END DO
55	>
56	> DO j = jts, min(jde,jte)
57	> DO i = its, min(ide-1,ite)
58	> grid%v10(i,j)=grid%em_v_gc(i,1,j)
59	> END DO
60	> END DO
61	>
62	467,477c439,443
63	< !!****MARS: decide to switch off this option
64	< !!****MARS: --> cf sfcprs2 and geopotential function at 500mb
65	< ! IF ( config_flags%adjust_heights ) THEN
66	< ! we_have_tavgsfc = ( flag_tavgsfc == 1 )
67	< ! ELSE
68	< ! we_have_tavgsfc = .FALSE.
69	< ! END IF
70	< !****MARS:
71	< we_have_tavgsfc = .FALSE.
72	<
73	<
74	---
75	> IF ( config_flags%adjust_heights ) THEN
76	> we_have_tavgsfc = ( flag_tavgsfc == 1 )
77	> ELSE
78	> we_have_tavgsfc = .FALSE.
79	> END IF
80	479d444
81	< !****MARS: hi-res psfc is done if the flag 'sfcp_to_sfcp' is active
82	482d446
83	< print *,'compute psfc from hi-res topography'
84	488,497c452,457
85	<
86	< !****MARS: no sea-level pressure inputs possible
87	< ! ELSE
88	< ! CALL sfcprs (grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_pslv_gc, grid%ht, &
89	< ! grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, &
90	< ! ids , ide , jds , jde , 1 , num_metgrid_levels , &
91	< ! ims , ime , jms , jme , 1 , num_metgrid_levels , &
92	< ! its , ite , jts , jte , 1 , num_metgrid_levels )
93	< !****MARS: no sea-level pressure inputs possible
94	<
95	---
96	> ELSE
97	> CALL sfcprs (grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_pslv_gc, grid%ht, &
98	> grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, &
99	> ids , ide , jds , jde , 1 , num_metgrid_levels , &
100	> ims , ime , jms , jme , 1 , num_metgrid_levels , &
101	> its , ite , jts , jte , 1 , num_metgrid_levels )
102	509d468
103	<
104	533,534d491
105	< !****MARS: em_dhs seems OK
106	<
107	567d523
108	<
109	580,583c536
110	< !****MARS: normalement c'est vert_interp
111	< !****MARS: mais les résultats sont trop discontinus > retour à une
112	< !****MARS: interpolation plus classique
113	< CALL vert_interp_old ( grid%em_qv_gc , grid%em_pd_gc , moist(:,:,:,P_QV) , grid%em_pb , &
114	---
115	> CALL vert_interp ( grid%em_qv_gc , grid%em_pd_gc , moist(:,:,:,P_QV) , grid%em_pb , &
116	590,592c543,544
117	<
118	< !****MARS: normalement c'est vert_interp
119	< CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2 , grid%em_pb , &
120	---
121	>
122	> CALL vert_interp ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2 , grid%em_pb , &
123	599d550
124	<
125	686,688c637,638
126	<
127	< !****MARS: normalement c'est vert_interp
128	< CALL vert_interp_old ( grid%em_u_gc , grid%em_pd_gc , grid%em_u_2 , grid%em_pb , &
129	---
130	>
131	> CALL vert_interp ( grid%em_u_gc , grid%em_pd_gc , grid%em_u_2 , grid%em_pb , &
132	695,696c645,646
133	< !****MARS: normalement c'est vert_interp
134	< CALL vert_interp_old ( grid%em_v_gc , grid%em_pd_gc , grid%em_v_2 , grid%em_pb , &
135	---
136	>
137	> CALL vert_interp ( grid%em_v_gc , grid%em_pd_gc , grid%em_v_2 , grid%em_pb , &
138	705a656,657
139	> ! Protect against bad grid%em_tsk values over water by supplying grid%sst (if it is
140	> ! available, and if the grid%sst is reasonable).
141	707,951c659,666
142	< !****MARS: no need
143	< ! ! Protect against bad grid%em_tsk values over water by supplying grid%sst (if it is
144	< ! ! available, and if the grid%sst is reasonable).
145	< !
146	< ! DO j = jts, MIN(jde-1,jte)
147	< ! DO i = its, MIN(ide-1,ite)
148	< ! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
149	< ! ( grid%sst(i,j) .GT. 200. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
150	< ! grid%tsk(i,j) = grid%sst(i,j)
151	< ! ENDIF
152	< ! END DO
153	< ! END DO
154	< !
155	< ! ! Save the grid%em_tsk field for later use in the sea ice surface temperature
156	< ! ! for the Noah LSM scheme.
157	< !
158	< ! DO j = jts, MIN(jte,jde-1)
159	< ! DO i = its, MIN(ite,ide-1)
160	< ! grid%tsk_save(i,j) = grid%tsk(i,j)
161	< ! END DO
162	< ! END DO
163	< !
164	< !!****MARS: no need
165	< ! ! Take the data from the input file and store it in the variables that
166	< ! ! use the WRF naming and ordering conventions.
167	< !
168	< ! DO j = jts, MIN(jte,jde-1)
169	< ! DO i = its, MIN(ite,ide-1)
170	< ! IF ( grid%snow(i,j) .GE. 10. ) then
171	< ! grid%snowc(i,j) = 1.
172	< ! ELSE
173	< ! grid%snowc(i,j) = 0.0
174	< ! END IF
175	< ! END DO
176	< ! END DO
177	< !
178	< ! ! Set flag integers for presence of snowh and soilw fields
179	< !
180	< ! grid%ifndsnowh = flag_snowh
181	< ! IF (num_sw_levels_input .GE. 1) THEN
182	< ! grid%ifndsoilw = 1
183	< ! ELSE
184	< ! grid%ifndsoilw = 0
185	< ! END IF
186	< !
187	< !****MARS: no need
188	< ! ! We require input data for the various LSM schemes.
189	< !
190	< ! enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
191	< !
192	< ! CASE (LSMSCHEME)
193	< ! IF ( num_st_levels_input .LT. 2 ) THEN
194	< ! CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
195	< ! END IF
196	< !
197	< ! CASE (RUCLSMSCHEME)
198	< ! IF ( num_st_levels_input .LT. 2 ) THEN
199	< ! CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
200	< ! END IF
201	< !
202	< ! END SELECT enough_data
203	< !
204	< ! ! For sf_surface_physics = 1, we want to use close to a 30 cm value
205	< ! ! for the bottom level of the soil temps.
206	< !
207	< ! fix_bottom_level_for_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
208	< !
209	< ! CASE (SLABSCHEME)
210	< ! IF ( flag_tavgsfc .EQ. 1 ) THEN
211	< ! DO j = jts , MIN(jde-1,jte)
212	< ! DO i = its , MIN(ide-1,ite)
213	< ! grid%tmn(i,j) = grid%em_tavgsfc(i,j)
214	< ! END DO
215	< ! END DO
216	< ! ELSE IF ( flag_st010040 .EQ. 1 ) THEN
217	< ! DO j = jts , MIN(jde-1,jte)
218	< ! DO i = its , MIN(ide-1,ite)
219	< ! grid%tmn(i,j) = grid%st010040(i,j)
220	< ! END DO
221	< ! END DO
222	< ! ELSE IF ( flag_st000010 .EQ. 1 ) THEN
223	< ! DO j = jts , MIN(jde-1,jte)
224	< ! DO i = its , MIN(ide-1,ite)
225	< ! grid%tmn(i,j) = grid%st000010(i,j)
226	< ! END DO
227	< ! END DO
228	< ! ELSE IF ( flag_soilt020 .EQ. 1 ) THEN
229	< ! DO j = jts , MIN(jde-1,jte)
230	< ! DO i = its , MIN(ide-1,ite)
231	< ! grid%tmn(i,j) = grid%soilt020(i,j)
232	< ! END DO
233	< ! END DO
234	< ! ELSE IF ( flag_st007028 .EQ. 1 ) THEN
235	< ! DO j = jts , MIN(jde-1,jte)
236	< ! DO i = its , MIN(ide-1,ite)
237	< ! grid%tmn(i,j) = grid%st007028(i,j)
238	< ! END DO
239	< ! END DO
240	< ! ELSE
241	< ! CALL wrf_debug ( 0 , 'No 10-40 cm, 0-10 cm, 7-28, or 20 cm soil temperature data for grid%em_tmn')
242	< ! CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
243	< ! END IF
244	< !
245	< ! CASE (LSMSCHEME)
246	< !
247	< ! CASE (RUCLSMSCHEME)
248	< !
249	< ! END SELECT fix_bottom_level_for_temp
250	< !
251	< ! ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
252	< ! ! is for the 5-layer scheme.
253	< !
254	< ! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
255	< ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
256	< ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
257	< ! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
258	< ! CALL nl_get_isice ( grid%id , grid%isice )
259	< ! CALL nl_get_iswater ( grid%id , grid%iswater )
260	< ! CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
261	< ! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
262	< ! grid%soilcbot , grid%tmn , &
263	< ! grid%seaice_threshold , &
264	< ! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
265	< ! grid%iswater , grid%isice , &
266	< ! model_config_rec%sf_surface_physics(grid%id) , &
267	< ! ids , ide , jds , jde , kds , kde , &
268	< ! ims , ime , jms , jme , kms , kme , &
269	< ! its , ite , jts , jte , kts , kte )
270	< !
271	< ! ! surface_input_source=1 => use data from static file (fractional category as input)
272	< ! ! surface_input_source=2 => use data from grib file (dominant category as input)
273	< !
274	< ! IF ( config_flags%surface_input_source .EQ. 1 ) THEN
275	< ! grid%vegcat (its,jts) = 0
276	< ! grid%soilcat(its,jts) = 0
277	< ! END IF
278	< !
279	< ! ! Generate the vegetation and soil category information from the fractional input
280	< ! ! data, or use the existing dominant category fields if they exist.
281	< !
282	< ! IF ( ( grid%soilcat(its,jts) .LT. 0.5 ) .AND. ( grid%vegcat(its,jts) .LT. 0.5 ) ) THEN
283	< !
284	< ! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
285	< ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
286	< ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
287	< !
288	< ! CALL process_percent_cat_new ( grid%landmask , &
289	< ! grid%landusef , grid%soilctop , grid%soilcbot , &
290	< ! grid%isltyp , grid%ivgtyp , &
291	< ! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
292	< ! ids , ide , jds , jde , kds , kde , &
293	< ! ims , ime , jms , jme , kms , kme , &
294	< ! its , ite , jts , jte , kts , kte , &
295	< ! model_config_rec%iswater(grid%id) )
296	< !
297	< ! ! Make all the veg/soil parms the same so as not to confuse the developer.
298	< !
299	< ! DO j = jts , MIN(jde-1,jte)
300	< ! DO i = its , MIN(ide-1,ite)
301	< ! grid%vegcat(i,j) = grid%ivgtyp(i,j)
302	< ! grid%soilcat(i,j) = grid%isltyp(i,j)
303	< ! END DO
304	< ! END DO
305	< !
306	< ! ELSE
307	< !
308	< ! ! Do we have dominant soil and veg data from the input already?
309	< !
310	< ! IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
311	< ! DO j = jts, MIN(jde-1,jte)
312	< ! DO i = its, MIN(ide-1,ite)
313	< ! grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
314	< ! END DO
315	< ! END DO
316	< ! END IF
317	< ! IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
318	< ! DO j = jts, MIN(jde-1,jte)
319	< ! DO i = its, MIN(ide-1,ite)
320	< ! grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
321	< ! END DO
322	< ! END DO
323	< ! END IF
324	< !
325	< ! END IF
326	< !
327	< ! ! Land use assignment.
328	< !
329	< ! DO j = jts, MIN(jde-1,jte)
330	< ! DO i = its, MIN(ide-1,ite)
331	< ! grid%lu_index(i,j) = grid%ivgtyp(i,j)
332	< ! IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
333	< ! grid%landmask(i,j) = 1
334	< ! grid%xland(i,j) = 1
335	< ! ELSE
336	< ! grid%landmask(i,j) = 0
337	< ! grid%xland(i,j) = 2
338	< ! END IF
339	< ! END DO
340	< ! END DO
341	< !
342	< ! ! Adjust the various soil temperature values depending on the difference in
343	< ! ! in elevation between the current model's elevation and the incoming data's
344	< ! ! orography.
345	< !
346	< ! IF ( flag_soilhgt .EQ. 1 ) THEN
347	< ! adjust_soil : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
348	< !
349	< ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
350	< ! CALL adjust_soil_temp_new ( grid%tmn , model_config_rec%sf_surface_physics(grid%id) , &
351	< ! grid%tsk , grid%ht , grid%toposoil , grid%landmask , flag_soilhgt , &
352	< ! grid%st000010 , grid%st010040 , grid%st040100 , grid%st100200 , grid%st010200 , &
353	< ! flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
354	< ! grid%st000007 , grid%st007028 , grid%st028100 , grid%st100255 , &
355	< ! flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
356	< ! grid%soilt000 , grid%soilt005 , grid%soilt020 , grid%soilt040 , grid%soilt160 , &
357	< ! grid%soilt300 , &
358	< ! flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , &
359	< ! flag_soilt160 , flag_soilt300 , &
360	< ! ids , ide , jds , jde , kds , kde , &
361	< ! ims , ime , jms , jme , kms , kme , &
362	< ! its , ite , jts , jte , kts , kte )
363	< !
364	< ! END SELECT adjust_soil
365	< ! END IF
366	< !
367	< ! ! Fix grid%em_tmn and grid%em_tsk.
368	< !
369	< ! fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
370	< !
371	< ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
372	< ! DO j = jts, MIN(jde-1,jte)
373	< ! DO i = its, MIN(ide-1,ite)
374	< ! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
375	< ! ( grid%sst(i,j) .GT. 240. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
376	< ! grid%tmn(i,j) = grid%sst(i,j)
377	< ! grid%tsk(i,j) = grid%sst(i,j)
378	< ! ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
379	< ! grid%tmn(i,j) = grid%tsk(i,j)
380	< ! END IF
381	< ! END DO
382	< ! END DO
383	< ! END SELECT fix_tsk_tmn
384	< !
385	< ! ! Is the grid%em_tsk reasonable?
386	< !
387	---
388	> DO j = jts, MIN(jde-1,jte)
389	> DO i = its, MIN(ide-1,ite)
390	> IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
391	> ( grid%sst(i,j) .GT. 200. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
392	> grid%tsk(i,j) = grid%sst(i,j)
393	> ENDIF
394	> END DO
395	> END DO
396	952a668,669
397	> ! Save the grid%em_tsk field for later use in the sea ice surface temperature
398	> ! for the Noah LSM scheme.
399	954c671,898
400	< !!**** MARS
401	---
402	> DO j = jts, MIN(jte,jde-1)
403	> DO i = its, MIN(ite,ide-1)
404	> grid%tsk_save(i,j) = grid%tsk(i,j)
405	> END DO
406	> END DO
407	>
408	> ! Take the data from the input file and store it in the variables that
409	> ! use the WRF naming and ordering conventions.
410	>
411	> DO j = jts, MIN(jte,jde-1)
412	> DO i = its, MIN(ite,ide-1)
413	> IF ( grid%snow(i,j) .GE. 10. ) then
414	> grid%snowc(i,j) = 1.
415	> ELSE
416	> grid%snowc(i,j) = 0.0
417	> END IF
418	> END DO
419	> END DO
420	>
421	> ! Set flag integers for presence of snowh and soilw fields
422	>
423	> grid%ifndsnowh = flag_snowh
424	> IF (num_sw_levels_input .GE. 1) THEN
425	> grid%ifndsoilw = 1
426	> ELSE
427	> grid%ifndsoilw = 0
428	> END IF
429	>
430	> ! We require input data for the various LSM schemes.
431	>
432	> enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
433	>
434	> CASE (LSMSCHEME)
435	> IF ( num_st_levels_input .LT. 2 ) THEN
436	> CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
437	> END IF
438	>
439	> CASE (RUCLSMSCHEME)
440	> IF ( num_st_levels_input .LT. 2 ) THEN
441	> CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
442	> END IF
443	>
444	> END SELECT enough_data
445	>
446	> ! For sf_surface_physics = 1, we want to use close to a 30 cm value
447	> ! for the bottom level of the soil temps.
448	>
449	> fix_bottom_level_for_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
450	>
451	> CASE (SLABSCHEME)
452	> IF ( flag_tavgsfc .EQ. 1 ) THEN
453	> DO j = jts , MIN(jde-1,jte)
454	> DO i = its , MIN(ide-1,ite)
455	> grid%tmn(i,j) = grid%em_tavgsfc(i,j)
456	> END DO
457	> END DO
458	> ELSE IF ( flag_st010040 .EQ. 1 ) THEN
459	> DO j = jts , MIN(jde-1,jte)
460	> DO i = its , MIN(ide-1,ite)
461	> grid%tmn(i,j) = grid%st010040(i,j)
462	> END DO
463	> END DO
464	> ELSE IF ( flag_st000010 .EQ. 1 ) THEN
465	> DO j = jts , MIN(jde-1,jte)
466	> DO i = its , MIN(ide-1,ite)
467	> grid%tmn(i,j) = grid%st000010(i,j)
468	> END DO
469	> END DO
470	> ELSE IF ( flag_soilt020 .EQ. 1 ) THEN
471	> DO j = jts , MIN(jde-1,jte)
472	> DO i = its , MIN(ide-1,ite)
473	> grid%tmn(i,j) = grid%soilt020(i,j)
474	> END DO
475	> END DO
476	> ELSE IF ( flag_st007028 .EQ. 1 ) THEN
477	> DO j = jts , MIN(jde-1,jte)
478	> DO i = its , MIN(ide-1,ite)
479	> grid%tmn(i,j) = grid%st007028(i,j)
480	> END DO
481	> END DO
482	> ELSE
483	> CALL wrf_debug ( 0 , 'No 10-40 cm, 0-10 cm, 7-28, or 20 cm soil temperature data for grid%em_tmn')
484	> CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
485	> END IF
486	>
487	> CASE (LSMSCHEME)
488	>
489	> CASE (RUCLSMSCHEME)
490	>
491	> END SELECT fix_bottom_level_for_temp
492	>
493	> ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
494	> ! is for the 5-layer scheme.
495	>
496	> num_veg_cat = SIZE ( grid%landusef , DIM=2 )
497	> num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
498	> num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
499	> CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
500	> CALL nl_get_isice ( grid%id , grid%isice )
501	> CALL nl_get_iswater ( grid%id , grid%iswater )
502	> CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
503	> grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
504	> grid%soilcbot , grid%tmn , &
505	> grid%seaice_threshold , &
506	> num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
507	> grid%iswater , grid%isice , &
508	> model_config_rec%sf_surface_physics(grid%id) , &
509	> ids , ide , jds , jde , kds , kde , &
510	> ims , ime , jms , jme , kms , kme , &
511	> its , ite , jts , jte , kts , kte )
512	>
513	> ! surface_input_source=1 => use data from static file (fractional category as input)
514	> ! surface_input_source=2 => use data from grib file (dominant category as input)
515	>
516	> IF ( config_flags%surface_input_source .EQ. 1 ) THEN
517	> grid%vegcat (its,jts) = 0
518	> grid%soilcat(its,jts) = 0
519	> END IF
520	>
521	> ! Generate the vegetation and soil category information from the fractional input
522	> ! data, or use the existing dominant category fields if they exist.
523	>
524	> IF ( ( grid%soilcat(its,jts) .LT. 0.5 ) .AND. ( grid%vegcat(its,jts) .LT. 0.5 ) ) THEN
525	>
526	> num_veg_cat = SIZE ( grid%landusef , DIM=2 )
527	> num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
528	> num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
529	>
530	> CALL process_percent_cat_new ( grid%landmask , &
531	> grid%landusef , grid%soilctop , grid%soilcbot , &
532	> grid%isltyp , grid%ivgtyp , &
533	> num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
534	> ids , ide , jds , jde , kds , kde , &
535	> ims , ime , jms , jme , kms , kme , &
536	> its , ite , jts , jte , kts , kte , &
537	> model_config_rec%iswater(grid%id) )
538	>
539	> ! Make all the veg/soil parms the same so as not to confuse the developer.
540	>
541	> DO j = jts , MIN(jde-1,jte)
542	> DO i = its , MIN(ide-1,ite)
543	> grid%vegcat(i,j) = grid%ivgtyp(i,j)
544	> grid%soilcat(i,j) = grid%isltyp(i,j)
545	> END DO
546	> END DO
547	>
548	> ELSE
549	>
550	> ! Do we have dominant soil and veg data from the input already?
551	>
552	> IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
553	> DO j = jts, MIN(jde-1,jte)
554	> DO i = its, MIN(ide-1,ite)
555	> grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
556	> END DO
557	> END DO
558	> END IF
559	> IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
560	> DO j = jts, MIN(jde-1,jte)
561	> DO i = its, MIN(ide-1,ite)
562	> grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
563	> END DO
564	> END DO
565	> END IF
566	>
567	> END IF
568	>
569	> ! Land use assignment.
570	>
571	> DO j = jts, MIN(jde-1,jte)
572	> DO i = its, MIN(ide-1,ite)
573	> grid%lu_index(i,j) = grid%ivgtyp(i,j)
574	> IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
575	> grid%landmask(i,j) = 1
576	> grid%xland(i,j) = 1
577	> ELSE
578	> grid%landmask(i,j) = 0
579	> grid%xland(i,j) = 2
580	> END IF
581	> END DO
582	> END DO
583	>
584	> ! Adjust the various soil temperature values depending on the difference in
585	> ! in elevation between the current model's elevation and the incoming data's
586	> ! orography.
587	>
588	> IF ( flag_soilhgt .EQ. 1 ) THEN
589	> adjust_soil : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
590	>
591	> CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
592	> CALL adjust_soil_temp_new ( grid%tmn , model_config_rec%sf_surface_physics(grid%id) , &
593	> grid%tsk , grid%ht , grid%toposoil , grid%landmask , flag_soilhgt , &
594	> grid%st000010 , grid%st010040 , grid%st040100 , grid%st100200 , grid%st010200 , &
595	> flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
596	> grid%st000007 , grid%st007028 , grid%st028100 , grid%st100255 , &
597	> flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
598	> grid%soilt000 , grid%soilt005 , grid%soilt020 , grid%soilt040 , grid%soilt160 , &
599	> grid%soilt300 , &
600	> flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , &
601	> flag_soilt160 , flag_soilt300 , &
602	> ids , ide , jds , jde , kds , kde , &
603	> ims , ime , jms , jme , kms , kme , &
604	> its , ite , jts , jte , kts , kte )
605	>
606	> END SELECT adjust_soil
607	> END IF
608	>
609	> ! Fix grid%em_tmn and grid%em_tsk.
610	>
611	> fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
612	>
613	> CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
614	> DO j = jts, MIN(jde-1,jte)
615	> DO i = its, MIN(ide-1,ite)
616	> IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
617	> ( grid%sst(i,j) .GT. 240. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
618	> grid%tmn(i,j) = grid%sst(i,j)
619	> grid%tsk(i,j) = grid%sst(i,j)
620	> ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
621	> grid%tmn(i,j) = grid%tsk(i,j)
622	> END IF
623	> END DO
624	> END DO
625	> END SELECT fix_tsk_tmn
626	>
627	> ! Is the grid%em_tsk reasonable?
628	>
629	> IF ( internal_time_loop .NE. 1 ) THEN
630	957,1248c901,1036
631	< !!grid%tsk(i,j)=200
632	< grid%tmn(i,j)=0
633	< grid%sst(i,j)=0 !!no use on Mars!!
634	< grid%tslb(i,j)=0
635	< END DO
636	< END DO
637	< !!**** MARS
638	<
639	< ! IF ( internal_time_loop .NE. 1 ) THEN
640	< ! DO j = jts, MIN(jde-1,jte)
641	< ! DO i = its, MIN(ide-1,ite)
642	< ! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
643	< ! grid%tsk(i,j) = grid%em_t_2(i,1,j)
644	< ! END IF
645	< ! END DO
646	< ! END DO
647	< ! ELSE
648	< ! DO j = jts, MIN(jde-1,jte)
649	< ! DO i = its, MIN(ide-1,ite)
650	< ! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
651	< ! print *,'error in the grid%em_tsk'
652	< ! print *,'i,j=',i,j
653	< ! print *,'grid%landmask=',grid%landmask(i,j)
654	< ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
655	< ! if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
656	< ! grid%tsk(i,j)=grid%tmn(i,j)
657	< ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
658	< ! grid%tsk(i,j)=grid%sst(i,j)
659	< ! else
660	< ! CALL wrf_error_fatal ( 'grid%em_tsk unreasonable' )
661	< ! end if
662	< ! END IF
663	< ! END DO
664	< ! END DO
665	< ! END IF
666	< !
667	< ! ! Is the grid%em_tmn reasonable?
668	< !
669	< ! DO j = jts, MIN(jde-1,jte)
670	< ! DO i = its, MIN(ide-1,ite)
671	< ! IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
672	< ! .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
673	< ! IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN
674	< ! print *,'error in the grid%em_tmn'
675	< ! print *,'i,j=',i,j
676	< ! print *,'grid%landmask=',grid%landmask(i,j)
677	< ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
678	< ! END IF
679	< !
680	< ! if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
681	< ! grid%tmn(i,j)=grid%tsk(i,j)
682	< ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
683	< ! grid%tmn(i,j)=grid%sst(i,j)
684	< ! else
685	< ! CALL wrf_error_fatal ( 'grid%em_tmn unreasonable' )
686	< ! endif
687	< ! END IF
688	< ! END DO
689	< ! END DO
690	< !
691	< ! interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
692	< !
693	< ! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
694	< ! CALL process_soil_real ( grid%tsk , grid%tmn , &
695	< ! grid%landmask , grid%sst , &
696	< ! st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
697	< ! grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
698	< ! flag_sst , flag_soilt000, flag_soilm000, &
699	< ! ids , ide , jds , jde , kds , kde , &
700	< ! ims , ime , jms , jme , kms , kme , &
701	< ! its , ite , jts , jte , kts , kte , &
702	< ! model_config_rec%sf_surface_physics(grid%id) , &
703	< ! model_config_rec%num_soil_layers , &
704	< ! model_config_rec%real_data_init_type , &
705	< ! num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
706	< ! num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
707	< !
708	< ! END SELECT interpolate_soil_tmw
709	< !
710	< ! ! Minimum soil values, residual, from RUC LSM scheme. For input from Noah and using
711	< ! ! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
712	< ! ! input RUC data and using the Noah LSM scheme, this value must be added to the soil
713	< ! ! moisture input.
714	< !
715	< ! lqmi(1:num_soil_top_cat) = &
716	< ! (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
717	< ! 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
718	< ! 0.004, 0.065 /)
719	< !! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
720	< !
721	< ! ! At the initial time we care about values of soil moisture and temperature, other times are
722	< ! ! ignored by the model, so we ignore them, too.
723	< !
724	< ! IF ( domain_ClockIsStartTime(grid) ) THEN
725	< ! account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
726	< !
727	< ! CASE ( LSMSCHEME )
728	< ! iicount = 0
729	< ! IF ( FLAG_SM000010 .EQ. 1 ) THEN
730	< ! DO j = jts, MIN(jde-1,jte)
731	< ! DO i = its, MIN(ide-1,ite)
732	< ! IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
733	< ! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
734	< ! print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
735	< ! iicount = iicount + 1
736	< ! grid%smois(i,:,j) = 0.005
737	< ! END IF
738	< ! END DO
739	< ! END DO
740	< ! IF ( iicount .GT. 0 ) THEN
741	< ! print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
742	< ! END IF
743	< ! ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
744	< ! DO j = jts, MIN(jde-1,jte)
745	< ! DO i = its, MIN(ide-1,ite)
746	< ! grid%smois(i,:,j) = grid%smois(i,:,j) + lqmi(grid%isltyp(i,j))
747	< ! END DO
748	< ! END DO
749	< ! DO j = jts, MIN(jde-1,jte)
750	< ! DO i = its, MIN(ide-1,ite)
751	< ! IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
752	< ! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
753	< ! print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
754	< ! iicount = iicount + 1
755	< ! grid%smois(i,:,j) = 0.005
756	< ! END IF
757	< ! END DO
758	< ! END DO
759	< ! IF ( iicount .GT. 0 ) THEN
760	< ! print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
761	< ! END IF
762	< ! END IF
763	< !
764	< ! CASE ( RUCLSMSCHEME )
765	< ! iicount = 0
766	< ! IF ( FLAG_SM000010 .EQ. 1 ) THEN
767	< ! DO j = jts, MIN(jde-1,jte)
768	< ! DO i = its, MIN(ide-1,ite)
769	< ! grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
770	< ! END DO
771	< ! END DO
772	< ! ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
773	< ! ! no op
774	< ! END IF
775	< !
776	< ! END SELECT account_for_zero_soil_moisture
777	< ! END IF
778	< !
779	< ! ! Is the grid%tslb reasonable?
780	< !
781	< ! IF ( internal_time_loop .NE. 1 ) THEN
782	< ! DO j = jts, MIN(jde-1,jte)
783	< ! DO ns = 1 , model_config_rec%num_soil_layers
784	< ! DO i = its, MIN(ide-1,ite)
785	< ! IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
786	< ! grid%tslb(i,ns,j) = grid%em_t_2(i,1,j)
787	< ! grid%smois(i,ns,j) = 0.3
788	< ! END IF
789	< ! END DO
790	< ! END DO
791	< ! END DO
792	< ! ELSE
793	< ! DO j = jts, MIN(jde-1,jte)
794	< ! DO i = its, MIN(ide-1,ite)
795	< ! IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
796	< ! ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
797	< ! IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
798	< ! ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ) ) THEN
799	< ! print *,'error in the grid%tslb'
800	< ! print *,'i,j=',i,j
801	< ! print *,'grid%landmask=',grid%landmask(i,j)
802	< ! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
803	< ! print *,'grid%tslb = ',grid%tslb(i,:,j)
804	< ! print *,'old grid%smois = ',grid%smois(i,:,j)
805	< ! grid%smois(i,1,j) = 0.3
806	< ! grid%smois(i,2,j) = 0.3
807	< ! grid%smois(i,3,j) = 0.3
808	< ! grid%smois(i,4,j) = 0.3
809	< ! END IF
810	< !
811	< ! IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
812	< ! (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
813	< ! fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
814	< ! CASE ( SLABSCHEME )
815	< ! DO ns = 1 , model_config_rec%num_soil_layers
816	< ! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
817	< ! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
818	< ! END DO
819	< ! CASE ( LSMSCHEME , RUCLSMSCHEME )
820	< ! CALL wrf_error_fatal ( 'Assigning constant soil moisture, bad idea')
821	< ! DO ns = 1 , model_config_rec%num_soil_layers
822	< ! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
823	< ! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
824	< ! END DO
825	< ! END SELECT fake_soil_temp
826	< ! else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
827	< ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
828	< ! DO ns = 1 , model_config_rec%num_soil_layers
829	< ! grid%tslb(i,ns,j)=grid%tsk(i,j)
830	< ! END DO
831	< ! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
832	< ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
833	< ! DO ns = 1 , model_config_rec%num_soil_layers
834	< ! grid%tslb(i,ns,j)=grid%sst(i,j)
835	< ! END DO
836	< ! else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
837	< ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
838	< ! DO ns = 1 , model_config_rec%num_soil_layers
839	< ! grid%tslb(i,ns,j)=grid%tmn(i,j)
840	< ! END DO
841	< ! else
842	< ! CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
843	< ! endif
844	< ! END IF
845	< ! END DO
846	< ! END DO
847	< ! END IF
848	< !
849	< ! ! Adjustments for the seaice field AFTER the grid%tslb computations. This is
850	< ! ! is for the Noah LSM scheme.
851	< !
852	< ! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
853	< ! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
854	< ! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
855	< ! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
856	< ! CALL nl_get_isice ( grid%id , grid%isice )
857	< ! CALL nl_get_iswater ( grid%id , grid%iswater )
858	< ! CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
859	< ! grid%ivgtyp , grid%vegcat , grid%lu_index , &
860	< ! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
861	< ! grid%soilctop , &
862	< ! grid%soilcbot , grid%tmn , grid%vegfra , &
863	< ! grid%tslb , grid%smois , grid%sh2o , &
864	< ! grid%seaice_threshold , &
865	< ! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
866	< ! model_config_rec%num_soil_layers , &
867	< ! grid%iswater , grid%isice , &
868	< ! model_config_rec%sf_surface_physics(grid%id) , &
869	< ! ids , ide , jds , jde , kds , kde , &
870	< ! ims , ime , jms , jme , kms , kme , &
871	< ! its , ite , jts , jte , kts , kte )
872	< !
873	< ! ! Let us make sure (again) that the grid%landmask and the veg/soil categories match.
874	< !
875	< !oops1=0
876	< !oops2=0
877	< ! DO j = jts, MIN(jde-1,jte)
878	< ! DO i = its, MIN(ide-1,ite)
879	< ! IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
880	< ! ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
881	< ! ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
882	< ! ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
883	< ! IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
884	< !oops1=oops1+1
885	< ! grid%ivgtyp(i,j) = 5
886	< ! grid%isltyp(i,j) = 8
887	< ! grid%landmask(i,j) = 1
888	< ! grid%xland(i,j) = 1
889	< ! ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
890	< !oops2=oops2+1
891	< ! grid%ivgtyp(i,j) = config_flags%iswater
892	< ! grid%isltyp(i,j) = 14
893	< ! grid%landmask(i,j) = 0
894	< ! grid%xland(i,j) = 2
895	< ! ELSE
896	< ! print *,'the grid%landmask and soil/veg cats do not match'
897	< ! print *,'i,j=',i,j
898	< ! print *,'grid%landmask=',grid%landmask(i,j)
899	< ! print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
900	< ! print *,'grid%isltyp=',grid%isltyp(i,j)
901	< ! print *,'iswater=', config_flags%iswater
902	< ! print *,'grid%tslb=',grid%tslb(i,:,j)
903	< ! print *,'grid%sst=',grid%sst(i,j)
904	< ! CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
905	< ! END IF
906	< ! END IF
907	< ! END DO
908	< ! END DO
909	< !if (oops1.gt.0) then
910	< !print *,'points artificially set to land : ',oops1
911	< !endif
912	< !if(oops2.gt.0) then
913	< !print *,'points artificially set to water: ',oops2
914	< !endif
915	< !! fill grid%sst array with grid%em_tsk if missing in real input (needed for time-varying grid%sst in wrf)
916	< ! DO j = jts, MIN(jde-1,jte)
917	< ! DO i = its, MIN(ide-1,ite)
918	< ! IF ( flag_sst .NE. 1 ) THEN
919	< ! grid%sst(i,j) = grid%tsk(i,j)
920	< ! ENDIF
921	< ! END DO
922	< ! END DO
923	---
924	> IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
925	> grid%tsk(i,j) = grid%em_t_2(i,1,j)
926	> END IF
927	> END DO
928	> END DO
929	> ELSE
930	> DO j = jts, MIN(jde-1,jte)
931	> DO i = its, MIN(ide-1,ite)
932	> IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
933	> print *,'error in the grid%em_tsk'
934	> print *,'i,j=',i,j
935	> print *,'grid%landmask=',grid%landmask(i,j)
936	> print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
937	> if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
938	> grid%tsk(i,j)=grid%tmn(i,j)
939	> else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
940	> grid%tsk(i,j)=grid%sst(i,j)
941	> else
942	> CALL wrf_error_fatal ( 'grid%em_tsk unreasonable' )
943	> end if
944	> END IF
945	> END DO
946	> END DO
947	> END IF
948	>
949	> ! Is the grid%em_tmn reasonable?
950	>
951	> DO j = jts, MIN(jde-1,jte)
952	> DO i = its, MIN(ide-1,ite)
953	> IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
954	> .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
955	> IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN
956	> print *,'error in the grid%em_tmn'
957	> print *,'i,j=',i,j
958	> print *,'grid%landmask=',grid%landmask(i,j)
959	> print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
960	> END IF
961	>
962	> if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
963	> grid%tmn(i,j)=grid%tsk(i,j)
964	> else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
965	> grid%tmn(i,j)=grid%sst(i,j)
966	> else
967	> CALL wrf_error_fatal ( 'grid%em_tmn unreasonable' )
968	> endif
969	> END IF
970	> END DO
971	> END DO
972	>
973	> interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
974	>
975	> CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
976	> CALL process_soil_real ( grid%tsk , grid%tmn , &
977	> grid%landmask , grid%sst , &
978	> st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
979	> grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
980	> flag_sst , flag_soilt000, flag_soilm000, &
981	> ids , ide , jds , jde , kds , kde , &
982	> ims , ime , jms , jme , kms , kme , &
983	> its , ite , jts , jte , kts , kte , &
984	> model_config_rec%sf_surface_physics(grid%id) , &
985	> model_config_rec%num_soil_layers , &
986	> model_config_rec%real_data_init_type , &
987	> num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
988	> num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
989	>
990	> END SELECT interpolate_soil_tmw
991	>
992	> ! Minimum soil values, residual, from RUC LSM scheme. For input from Noah and using
993	> ! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
994	> ! input RUC data and using the Noah LSM scheme, this value must be added to the soil
995	> ! moisture input.
996	>
997	> lqmi(1:num_soil_top_cat) = &
998	> (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
999	> 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1000	> 0.004, 0.065 /)
1001	> ! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
1002	>
1003	> ! At the initial time we care about values of soil moisture and temperature, other times are
1004	> ! ignored by the model, so we ignore them, too.
1005	>
1006	> IF ( domain_ClockIsStartTime(grid) ) THEN
1007	> account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1008	>
1009	> CASE ( LSMSCHEME )
1010	> iicount = 0
1011	> IF ( FLAG_SM000010 .EQ. 1 ) THEN
1012	> DO j = jts, MIN(jde-1,jte)
1013	> DO i = its, MIN(ide-1,ite)
1014	> IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
1015	> ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1016	> print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1017	> iicount = iicount + 1
1018	> grid%smois(i,:,j) = 0.005
1019	> END IF
1020	> END DO
1021	> END DO
1022	> IF ( iicount .GT. 0 ) THEN
1023	> print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
1024	> END IF
1025	> ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1026	> DO j = jts, MIN(jde-1,jte)
1027	> DO i = its, MIN(ide-1,ite)
1028	> grid%smois(i,:,j) = grid%smois(i,:,j) + lqmi(grid%isltyp(i,j))
1029	> END DO
1030	> END DO
1031	> DO j = jts, MIN(jde-1,jte)
1032	> DO i = its, MIN(ide-1,ite)
1033	> IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
1034	> ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1035	> print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1036	> iicount = iicount + 1
1037	> grid%smois(i,:,j) = 0.005
1038	> END IF
1039	> END DO
1040	> END DO
1041	> IF ( iicount .GT. 0 ) THEN
1042	> print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
1043	> END IF
1044	> END IF
1045	>
1046	> CASE ( RUCLSMSCHEME )
1047	> iicount = 0
1048	> IF ( FLAG_SM000010 .EQ. 1 ) THEN
1049	> DO j = jts, MIN(jde-1,jte)
1050	> DO i = its, MIN(ide-1,ite)
1051	> grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
1052	> END DO
1053	> END DO
1054	> ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1055	> ! no op
1056	> END IF
1057	>
1058	> END SELECT account_for_zero_soil_moisture
1059	> END IF
1060	1249a1038,1181
1061	> ! Is the grid%tslb reasonable?
1062	>
1063	> IF ( internal_time_loop .NE. 1 ) THEN
1064	> DO j = jts, MIN(jde-1,jte)
1065	> DO ns = 1 , model_config_rec%num_soil_layers
1066	> DO i = its, MIN(ide-1,ite)
1067	> IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
1068	> grid%tslb(i,ns,j) = grid%em_t_2(i,1,j)
1069	> grid%smois(i,ns,j) = 0.3
1070	> END IF
1071	> END DO
1072	> END DO
1073	> END DO
1074	> ELSE
1075	> DO j = jts, MIN(jde-1,jte)
1076	> DO i = its, MIN(ide-1,ite)
1077	> IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
1078	> ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1079	> IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
1080	> ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ) ) THEN
1081	> print *,'error in the grid%tslb'
1082	> print *,'i,j=',i,j
1083	> print *,'grid%landmask=',grid%landmask(i,j)
1084	> print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1085	> print *,'grid%tslb = ',grid%tslb(i,:,j)
1086	> print *,'old grid%smois = ',grid%smois(i,:,j)
1087	> grid%smois(i,1,j) = 0.3
1088	> grid%smois(i,2,j) = 0.3
1089	> grid%smois(i,3,j) = 0.3
1090	> grid%smois(i,4,j) = 0.3
1091	> END IF
1092	>
1093	> IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
1094	> (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
1095	> fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1096	> CASE ( SLABSCHEME )
1097	> DO ns = 1 , model_config_rec%num_soil_layers
1098	> grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1099	> grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1100	> END DO
1101	> CASE ( LSMSCHEME , RUCLSMSCHEME )
1102	> CALL wrf_error_fatal ( 'Assigning constant soil moisture, bad idea')
1103	> DO ns = 1 , model_config_rec%num_soil_layers
1104	> grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1105	> grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1106	> END DO
1107	> END SELECT fake_soil_temp
1108	> else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1109	> CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
1110	> DO ns = 1 , model_config_rec%num_soil_layers
1111	> grid%tslb(i,ns,j)=grid%tsk(i,j)
1112	> END DO
1113	> else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1114	> CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
1115	> DO ns = 1 , model_config_rec%num_soil_layers
1116	> grid%tslb(i,ns,j)=grid%sst(i,j)
1117	> END DO
1118	> else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1119	> CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
1120	> DO ns = 1 , model_config_rec%num_soil_layers
1121	> grid%tslb(i,ns,j)=grid%tmn(i,j)
1122	> END DO
1123	> else
1124	> CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
1125	> endif
1126	> END IF
1127	> END DO
1128	> END DO
1129	> END IF
1130	>
1131	> ! Adjustments for the seaice field AFTER the grid%tslb computations. This is
1132	> ! is for the Noah LSM scheme.
1133	>
1134	> num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1135	> num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1136	> num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1137	> CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1138	> CALL nl_get_isice ( grid%id , grid%isice )
1139	> CALL nl_get_iswater ( grid%id , grid%iswater )
1140	> CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
1141	> grid%ivgtyp , grid%vegcat , grid%lu_index , &
1142	> grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
1143	> grid%soilctop , &
1144	> grid%soilcbot , grid%tmn , grid%vegfra , &
1145	> grid%tslb , grid%smois , grid%sh2o , &
1146	> grid%seaice_threshold , &
1147	> num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1148	> model_config_rec%num_soil_layers , &
1149	> grid%iswater , grid%isice , &
1150	> model_config_rec%sf_surface_physics(grid%id) , &
1151	> ids , ide , jds , jde , kds , kde , &
1152	> ims , ime , jms , jme , kms , kme , &
1153	> its , ite , jts , jte , kts , kte )
1154	>
1155	> ! Let us make sure (again) that the grid%landmask and the veg/soil categories match.
1156	>
1157	> oops1=0
1158	> oops2=0
1159	> DO j = jts, MIN(jde-1,jte)
1160	> DO i = its, MIN(ide-1,ite)
1161	> IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
1162	> ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
1163	> ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
1164	> ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
1165	> IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
1166	> oops1=oops1+1
1167	> grid%ivgtyp(i,j) = 5
1168	> grid%isltyp(i,j) = 8
1169	> grid%landmask(i,j) = 1
1170	> grid%xland(i,j) = 1
1171	> ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
1172	> oops2=oops2+1
1173	> grid%ivgtyp(i,j) = config_flags%iswater
1174	> grid%isltyp(i,j) = 14
1175	> grid%landmask(i,j) = 0
1176	> grid%xland(i,j) = 2
1177	> ELSE
1178	> print *,'the grid%landmask and soil/veg cats do not match'
1179	> print *,'i,j=',i,j
1180	> print *,'grid%landmask=',grid%landmask(i,j)
1181	> print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
1182	> print *,'grid%isltyp=',grid%isltyp(i,j)
1183	> print *,'iswater=', config_flags%iswater
1184	> print *,'grid%tslb=',grid%tslb(i,:,j)
1185	> print *,'grid%sst=',grid%sst(i,j)
1186	> CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
1187	> END IF
1188	> END IF
1189	> END DO
1190	> END DO
1191	> if (oops1.gt.0) then
1192	> print *,'points artificially set to land : ',oops1
1193	> endif
1194	> if(oops2.gt.0) then
1195	> print *,'points artificially set to water: ',oops2
1196	> endif
1197	> ! fill grid%sst array with grid%em_tsk if missing in real input (needed for time-varying grid%sst in wrf)
1198	> DO j = jts, MIN(jde-1,jte)
1199	> DO i = its, MIN(ide-1,ite)
1200	> IF ( flag_sst .NE. 1 ) THEN
1201	> grid%sst(i,j) = grid%tsk(i,j)
1202	> ENDIF
1203	> END DO
1204	> END DO
1205	1348,1351d1279
1206	<
1207	< !****MARS
1208	< !TODO: étudier si une meilleure formule n'existe pas pour Mars
1209	< !****MARS
1210	1357c1285
1211	<
1212	---
1213	>
1214	1457,1469c1385,1391
1215	< !!---------------------------------------------------------------
1216	< !!****MARS: no 500mb adjustment needed
1217	< !!****MARS: must keep however the hydrostatic equation integration performed in this loop !
1218	< !!****MARS: the DO WHILE loop is deactivated, since we will always be in the case
1219	< !!****MARS: ... of "ELSE dpmu = 0."
1220	< !!---------------------------------------------------------------
1221	< ! dpmu = 10001.
1222	< ! loop_count = 0
1223	< !
1224	< ! DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
1225	< ! ( loop_count .LT. 5 ) )
1226	< !
1227	< ! loop_count = loop_count + 1
1228	---
1229	> dpmu = 10001.
1230	> loop_count = 0
1231	>
1232	> DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
1233	> ( loop_count .LT. 5 ) )
1234	>
1235	> loop_count = loop_count + 1
1236	1490c1412
1237	< DO k=kte-2,1,-1
1238	---
1239	> DO k=kte-2,1,-1
1240	1509a1432,1495
1241	>
1242	> ! Adjust the column pressure so that the computed 500 mb height is close to the
1243	> ! input value (of course, not when we are doing hybrid input).
1244	>
1245	> IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. its ) .AND. ( j .EQ. jts ) ) THEN
1246	> DO k = 1 , num_metgrid_levels
1247	> IF ( ABS ( grid%em_p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
1248	> lev500 = k
1249	> EXIT
1250	> END IF
1251	> END DO
1252	> END IF
1253	>
1254	> ! We only do the adjustment of height if we have the input data on pressure
1255	> ! surfaces, and folks have asked to do this option.
1256	>
1257	> IF ( ( flag_metgrid .EQ. 1 ) .AND. &
1258	> ( config_flags%adjust_heights ) .AND. &
1259	> ( lev500 .NE. 0 ) ) THEN
1260	>
1261	> DO k = 2 , kte-1
1262	>
1263	> ! Get the pressures on the full eta levels (grid%em_php is defined above as
1264	> ! the full-lev base pressure, an easy array to use for 3d space).
1265	>
1266	> pl = grid%em_php(i,k ,j) + &
1267	> ( grid%em_p(i,k-1 ,j) * ( grid%em_znw(k ) - grid%em_znu(k ) ) + &
1268	> grid%em_p(i,k ,j) * ( grid%em_znu(k-1 ) - grid%em_znw(k ) ) ) / &
1269	> ( grid%em_znu(k-1 ) - grid%em_znu(k ) )
1270	> pu = grid%em_php(i,k+1,j) + &
1271	> ( grid%em_p(i,k-1+1,j) * ( grid%em_znw(k +1) - grid%em_znu(k+1) ) + &
1272	> grid%em_p(i,k +1,j) * ( grid%em_znu(k-1+1) - grid%em_znw(k+1) ) ) / &
1273	> ( grid%em_znu(k-1+1) - grid%em_znu(k+1) )
1274	>
1275	> ! If these pressure levels trap 500 mb, use them to interpolate
1276	> ! to the 500 mb level of the computed height.
1277	>
1278	> IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
1279	> zl = ( grid%em_ph_2(i,k ,j) + grid%em_phb(i,k ,j) ) / g
1280	> zu = ( grid%em_ph_2(i,k+1,j) + grid%em_phb(i,k+1,j) ) / g
1281	>
1282	> z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
1283	> zu * ( LOG(pl ) - LOG(50000.) ) ) / &
1284	> ( LOG(pl) - LOG(pu) )
1285	> ! z500 = ( zl * ( (50000.) - (pu ) ) + &
1286	> ! zu * ( (pl ) - (50000.) ) ) / &
1287	> ! ( (pl) - (pu) )
1288	>
1289	> ! Compute the difference of the 500 mb heights (computed minus input), and
1290	> ! then the change in grid%em_mu_2. The grid%em_php is still full-levels, base pressure.
1291	>
1292	> dz500 = z500 - grid%em_ght_gc(i,lev500,j)
1293	> tvsfc = ((grid%em_t_2(i,1,j)+t0)((grid%em_p(i,1,j)+grid%em_php(i,1,j))/p1000mb)(r_d/cp)) &
1294	> (1.+0.6*moist(i,1,j,P_QV))
1295	> dpmu = ( grid%em_php(i,1,j) + grid%em_p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
1296	> dpmu = dpmu - ( grid%em_php(i,1,j) + grid%em_p(i,1,j) )
1297	> grid%em_mu_2(i,j) = grid%em_mu_2(i,j) - dpmu
1298	> EXIT
1299	> END IF
1300	>
1301	> END DO
1302	> ELSE
1303	> dpmu = 0.
1304	> END IF
1305	1511,1575c1497
1306	< ! ! Adjust the column pressure so that the computed 500 mb height is close to the
1307	< ! ! input value (of course, not when we are doing hybrid input).
1308	< !
1309	< ! IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. its ) .AND. ( j .EQ. jts ) ) THEN
1310	< ! DO k = 1 , num_metgrid_levels
1311	< ! IF ( ABS ( grid%em_p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
1312	< ! lev500 = k
1313	< ! EXIT
1314	< ! END IF
1315	< ! END DO
1316	< ! END IF
1317	< !
1318	< ! ! We only do the adjustment of height if we have the input data on pressure
1319	< ! ! surfaces, and folks have asked to do this option.
1320	< !
1321	< ! IF ( ( flag_metgrid .EQ. 1 ) .AND. &
1322	< ! ( config_flags%adjust_heights ) .AND. &
1323	< ! ( lev500 .NE. 0 ) ) THEN
1324	< !
1325	< ! DO k = 2 , kte-1
1326	< !
1327	< ! ! Get the pressures on the full eta levels (grid%em_php is defined above as
1328	< ! ! the full-lev base pressure, an easy array to use for 3d space).
1329	< !
1330	< ! pl = grid%em_php(i,k ,j) + &
1331	< ! ( grid%em_p(i,k-1 ,j) * ( grid%em_znw(k ) - grid%em_znu(k ) ) + &
1332	< ! grid%em_p(i,k ,j) * ( grid%em_znu(k-1 ) - grid%em_znw(k ) ) ) / &
1333	< ! ( grid%em_znu(k-1 ) - grid%em_znu(k ) )
1334	< ! pu = grid%em_php(i,k+1,j) + &
1335	< ! ( grid%em_p(i,k-1+1,j) * ( grid%em_znw(k +1) - grid%em_znu(k+1) ) + &
1336	< ! grid%em_p(i,k +1,j) * ( grid%em_znu(k-1+1) - grid%em_znw(k+1) ) ) / &
1337	< ! ( grid%em_znu(k-1+1) - grid%em_znu(k+1) )
1338	< !
1339	< ! ! If these pressure levels trap 500 mb, use them to interpolate
1340	< ! ! to the 500 mb level of the computed height.
1341	< !!**** PB on MARS .... ?
1342	< ! IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
1343	< ! zl = ( grid%em_ph_2(i,k ,j) + grid%em_phb(i,k ,j) ) / g
1344	< ! zu = ( grid%em_ph_2(i,k+1,j) + grid%em_phb(i,k+1,j) ) / g
1345	< !
1346	< ! z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
1347	< ! zu * ( LOG(pl ) - LOG(50000.) ) ) / &
1348	< ! ( LOG(pl) - LOG(pu) )
1349	< !! z500 = ( zl * ( (50000.) - (pu ) ) + &
1350	< !! zu * ( (pl ) - (50000.) ) ) / &
1351	< !! ( (pl) - (pu) )
1352	< !
1353	< ! ! Compute the difference of the 500 mb heights (computed minus input), and
1354	< ! ! then the change in grid%em_mu_2. The grid%em_php is still full-levels, base pressure.
1355	< !
1356	< ! dz500 = z500 - grid%em_ght_gc(i,lev500,j)
1357	< ! tvsfc = ((grid%em_t_2(i,1,j)+t0)((grid%em_p(i,1,j)+grid%em_php(i,1,j))/p1000mb)(r_d/cp)) &
1358	< ! (1.+0.6*moist(i,1,j,P_QV))
1359	< ! dpmu = ( grid%em_php(i,1,j) + grid%em_p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
1360	< ! dpmu = dpmu - ( grid%em_php(i,1,j) + grid%em_p(i,1,j) )
1361	< ! grid%em_mu_2(i,j) = grid%em_mu_2(i,j) - dpmu
1362	< ! EXIT
1363	< ! END IF
1364	< !
1365	< ! END DO
1366	< ! ELSE
1367	< ! dpmu = 0.
1368	< ! END IF
1369	< !
1370	< ! END DO
1371	---
1372	> END DO
1373	1580,1619c1502,1537
1374	< !!****MARS: we use WPS
1375	< !
1376	< ! ! If this is data from the SI, then we probably do not have the original
1377	< ! ! surface data laying around. Note that these are all the lowest levels
1378	< ! ! of the respective 3d arrays. For surface pressure, we assume that the
1379	< ! ! vertical gradient of grid%em_p prime is zilch. This is not all that important.
1380	< ! ! These are filled in so that the various plotting routines have something
1381	< ! ! to play with at the initial time for the model.
1382	< !
1383	< ! IF ( flag_metgrid .NE. 1 ) THEN
1384	< ! DO j = jts, min(jde-1,jte)
1385	< ! DO i = its, min(ide,ite)
1386	< ! grid%u10(i,j)=grid%em_u_2(i,1,j)
1387	< ! END DO
1388	< ! END DO
1389	< !
1390	< ! DO j = jts, min(jde,jte)
1391	< ! DO i = its, min(ide-1,ite)
1392	< ! grid%v10(i,j)=grid%em_v_2(i,1,j)
1393	< ! END DO
1394	< ! END DO
1395	< !
1396	< ! DO j = jts, min(jde-1,jte)
1397	< ! DO i = its, min(ide-1,ite)
1398	< ! p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
1399	< ! grid%psfc(i,j)=p_surf + grid%em_p(i,1,j)
1400	< ! grid%q2(i,j)=moist(i,1,j,P_QV)
1401	< ! grid%th2(i,j)=grid%em_t_2(i,1,j)+300.
1402	< ! grid%t2(i,j)=grid%th2(i,j)(((grid%em_p(i,1,j)+grid%em_pb(i,1,j))/p00)*(r_d/cp))
1403	< ! END DO
1404	< ! END DO
1405	< !
1406	< ! ! If this data is from WPS, then we have previously assigned the surface
1407	< ! ! data for u, v, and t. If we have an input qv, welp, we assigned that one,
1408	< ! ! too. Now we pick up the left overs, and if RH came in - we assign the
1409	< ! ! mixing ratio.
1410	< !
1411	< ! ELSE IF ( flag_metgrid .EQ. 1 ) THEN
1412	< !
1413	< !!****MARS: we use WPS
1414	---
1415	> ! If this is data from the SI, then we probably do not have the original
1416	> ! surface data laying around. Note that these are all the lowest levels
1417	> ! of the respective 3d arrays. For surface pressure, we assume that the
1418	> ! vertical gradient of grid%em_p prime is zilch. This is not all that important.
1419	> ! These are filled in so that the various plotting routines have something
1420	> ! to play with at the initial time for the model.
1421	>
1422	> IF ( flag_metgrid .NE. 1 ) THEN
1423	> DO j = jts, min(jde-1,jte)
1424	> DO i = its, min(ide,ite)
1425	> grid%u10(i,j)=grid%em_u_2(i,1,j)
1426	> END DO
1427	> END DO
1428	>
1429	> DO j = jts, min(jde,jte)
1430	> DO i = its, min(ide-1,ite)
1431	> grid%v10(i,j)=grid%em_v_2(i,1,j)
1432	> END DO
1433	> END DO
1434	>
1435	> DO j = jts, min(jde-1,jte)
1436	> DO i = its, min(ide-1,ite)
1437	> p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
1438	> grid%psfc(i,j)=p_surf + grid%em_p(i,1,j)
1439	> grid%q2(i,j)=moist(i,1,j,P_QV)
1440	> grid%th2(i,j)=grid%em_t_2(i,1,j)+300.
1441	> grid%t2(i,j)=grid%th2(i,j)(((grid%em_p(i,1,j)+grid%em_pb(i,1,j))/p00)*(r_d/cp))
1442	> END DO
1443	> END DO
1444	>
1445	> ! If this data is from WPS, then we have previously assigned the surface
1446	> ! data for u, v, and t. If we have an input qv, welp, we assigned that one,
1447	> ! too. Now we pick up the left overs, and if RH came in - we assign the
1448	> ! mixing ratio.
1449	>
1450	> ELSE IF ( flag_metgrid .EQ. 1 ) THEN
1451	1636c1554
1452	< ! END IF
1453	---
1454	> END IF
1455	2186,2192d2103
1456	< !****MARS
1457	< !big problems ... discontinuity in the interpolated fields ...
1458	< print *, '25/05/2007: decided to use simple linear interpolations'
1459	< stop
1460	< !****MARS
1461	<
1462	<
1463	2616d2526
1464	< !****MARS
1465	2619d2528
1466	< !****MARS
1467	2621c2530
1468	< ! Horizontal loop bounds for different variable types.
1469	---
1470	> ! Horiontal loop bounds for different variable types.
1471	2765d2673
1472	<
1473	2778d2685
1474	<
1475	2782,2794d2688
1476	< !!****MARS
1477	< !!
1478	< !! Pressure level may be OK, however data from the diagfi is possibly missing
1479	< IF (forig(i,ko,j) .EQ. -1.0e+30) THEN
1480	< ko_above_sfc(i) = -1
1481	< END IF
1482	< !! Once the right start level is found, check that it is OK
1483	< !! >> first column should be 1e30 or so, second column should be a realistic value
1484	< !IF ( ko_above_sfc(i) .NE. -1 ) THEN
1485	< ! print *, 'verif', forig(i,ko-1,j), forig(i,ko,j), forig(i,ko+1,j), ko
1486	< !END IF
1487	< !!
1488	< !!****MARS
1489	2797d2690
1490	<
1491	2843,2844d2735
1492	< !****MARS
1493	< !the possible issue is fixed later in the code ...
1494	2847d2737
1495	< !****MARS
1496	2882,2885d2771
1497	< !!****MARS
1498	< !!check: values are usually quite close
1499	< !print *,porig(i,1,j),pnew(i,kn,j)
1500	< !!****MARS
1501	2905,2910d2790
1502	< !!****MARS
1503	< k_above(i,kn) = 1
1504	< ks(i) = 1
1505	< !!"Hopefully, we are not extrapolating too far"
1506	< !!>> true on Mars ??
1507	< !!****MARS
1508	2941,2958d2820
1509	< !!****MARS
1510	< !!ne doit pas arriver avec la temperature si l'on definit bien le champ au sol
1511	< IF (forig(i,1,j) .EQ. -1.0e+30) THEN
1512	< print *,'no data here - surface - var is ...',var_type,i,j,1
1513	< print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j)
1514	< forig(i,1,j) = forig(i,ko_above_sfc(i),j)
1515	< !IF ( ( var_type .EQ. 'U' ) .OR. &
1516	< ! ( var_type .EQ. 'V' ) .OR. &
1517	< ! ( var_type .EQ. 'Q' ) ) THEN
1518	< ! print *,'zero wind at the ground'
1519	< ! forig(i,1,j) = 0
1520	< !ENDIF
1521	< IF (forig(i,1,j) .EQ. -1.0e+30) THEN
1522	< print *,'well ... are you sure ?'
1523	< stop
1524	< ENDIF
1525	< END IF
1526	< !!****MARS
1527	2966,2979d2827
1528	< !!****MARS
1529	< IF (forig(i,k2,j) .EQ. -1.0e+30) THEN
1530	< print *,'no data here - level above - you_d better stop',i,j,k2
1531	< stop
1532	< END IF
1533	< IF (forig(i,k1,j) .EQ. -1.0e+30) THEN
1534	< print *,'no data here - level below - var is ...',var_type,i,j,k1
1535	< print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j)
1536	< forig(i,k1,j) = forig(i,ko_above_sfc(i),j)
1537	< !!!VERIFIER QUE LA TEMPERATURE AU SOL N'EST PAS CONCERNEE
1538	< !!!(montagnes=sources locales de chaleur)
1539	< !!!normalement, pas de souci, et lors de l'exécution rien ne s'affiche
1540	< END IF
1541	< !!****MARS
1542	3026d2873
1543	< print *,'finished with ... ', var_type
1544	3062d2908
1545	<
1546	3089,3097d2934
1547	< !****MARS: check if no errors here
1548	< !print *,'interpolating ... ',var_type
1549	< ! print *,'i,j = ',i,j
1550	< ! print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
1551	< ! DO loop = 1 , all_dim
1552	< ! print *,'column of pressure and value = ',all_x(loop),all_y(loop)
1553	< ! END DO
1554	< !END IF
1555	< !****MARS
1556	3118,3119d2954
1557	< !****MARS: normally, no errors here (otherwise, keep this part commented ?)
1558	< print *, var_type
1559	3125,3126c2960
1560	< CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
1561	< !****MARS: end of 'keep this part commented'
1562	---
1563	> CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
1564	3313,3317c3147,3149
1565	< !****MARS ....
1566	< REAL , PARAMETER :: Rd = 192.
1567	< REAL , PARAMETER :: g = 3.72
1568	< print *,'compute dry, hydrostatic surface pressure'
1569	< !****MARS ....
1570	---
1571	>
1572	> REAL , PARAMETER :: Rd = 287.
1573	> REAL , PARAMETER :: g = 9.8
1574	3325,3334d3156
1575	< !****MARS
1576	< !****MARS cette formule est-elle juste sur Mars ?
1577	< !****MARS >> a première vue, ne donne pas de résultats absurdes
1578	< !****TODO: il y a peut être meilleur !
1579	< !****MARS
1580	<
1581	< !print *,pdhs
1582	< !stop
1583	<
1584	<
1585	3408,3412c3230,3233
1586	< !****MARS
1587	< REAL , PARAMETER :: Rd = 192.
1588	< REAL , PARAMETER :: Cp = 844.6
1589	< !****MARS
1590	<
1591	---
1592	>
1593	> REAL , PARAMETER :: Rd = 287.
1594	> REAL , PARAMETER :: Cp = 1004.
1595	>
1596	3456,3460c3277,3278
1597	< !****MARS
1598	< REAL , PARAMETER :: Rd = 192.
1599	< REAL , PARAMETER :: g = 3.72
1600	< !****MARS
1601	<
1602	---
1603	> REAL , PARAMETER :: Rd = 287.
1604	> REAL , PARAMETER :: g = 9.8
1605	3597,3605d3414
1606	<
1607	< !!****MARS: no water vapor pressure
1608	< ! DO k = level_above_sfc(i)-1,kts+1,-1
1609	< ! pd(i,k) = p(i,k)
1610	< ! END DO
1611	< ! pd(i,kts) = psfc(i,j)
1612	< !!****MARS
1613	<
1614	<
1615	3632,3633c3441
1616	< !****MARS .... à régler si besoin ....
1617	< !****MARS
1618	---
1619	>
1620	3663,3665d3470
1621	< !****MARS
1622	< !****MARS
1623	<
1624	3743,3749d3547
1625	< !****MARS
1626	< !****TEMPORARY
1627	< !****TEMPORARY
1628	< !TODO: change once tracers are activated ?
1629	< q=0.
1630	< !****MARS
1631	<
1632	3788,3796d3585
1633	< !****MARS
1634	< !****MARS
1635	< print *, 'check Mars: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0'
1636	< print *, p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0
1637	< !-----solution alternative: définir dans la namelist les niveaux verticaux
1638	< !****MARS
1639	< !****MARS
1640	<
1641	<
1642	3823,3843c3612,3613
1643	< ! znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
1644	< ! 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
1645	<
1646	< !****MARS
1647	< !****MARS
1648	< ! on Mars, this is important to correctly resolve the surface
1649	< ! -- levels were changed to get closer to the surface
1650	< ! -- values were chosen as done typically in LMD GCM simulations
1651	< !TODO: better repartition ?
1652	< znw_prac = (/ 1.000 , &
1653	< 0.9995 , &
1654	< 0.9980 , &
1655	< 0.9950 , &
1656	< 0.9850 , &
1657	< 0.9700 , &
1658	< 0.9400 , &
1659	< 0.9000 , &
1660	< 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
1661	< !****MARS
1662	< !****MARS
1663	<
1664	---
1665	> znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
1666	> 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
1667	3876a3647
1668	>
1669	3901d3671
1670	<
1671	3911,3920d3680
1672	< !!****MARS
1673	< !!attention 'base_lapse' ne doit pas être trop grand
1674	< !!sinon ... des NaN car températures négatives en haut
1675	< !IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
1676	< ! IF (k .EQ. 8) THEN
1677	< ! print *, 'p,t,z,k'
1678	< ! END IF
1679	< ! print *, pb,temp,znw(k+1),k
1680	< !END IF
1681	< !****MARS
1682	3950,3960d3709
1683	<
1684	< ! ****MARS
1685	< ! Display the computed levels
1686	< print *,'WRF levels are:'
1687	< print *,'z (m) = ',phb(1)/g
1688	< do k = 2 ,kte
1689	< print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
1690	< end do
1691	< ! ****MARS
1692	<
1693	<
1694	4108,4115c3857
1695	< !****MARS
1696	< REAL , PARAMETER :: Rd = 192.
1697	< REAL , PARAMETER :: Cp = 844.6
1698	< REAL, PARAMETER :: g = 3.72
1699	< REAL, PARAMETER :: pconst = 610.
1700	< !****MARS
1701	<
1702	< !****MARS .... to be changed if used
1703	---
1704	> REAL, PARAMETER :: g = 9.8
1705	4116a3859,3860
1706	> REAL, PARAMETER :: pconst = 10000.0
1707	> REAL, PARAMETER :: Rd = 287.
1708	4117a3862
1709	>
1710	4120d3864
1711	< !****MARS .... to be changed if used
1712	4158,4163d3901
1713	< !****MARS ....
1714	< !****MARS .... the mean sea level method is abandoned
1715	< print *, 'no sea level pressure on Mars, please'
1716	< stop
1717	< !****MARS ....
1718	<
1719	4412,4415c4150,4151
1720	< !****MARS
1721	< REAL , PARAMETER :: Rd = 192.
1722	< REAL, PARAMETER :: g = 3.72
1723	< !****MARS
1724	---
1725	> REAL, PARAMETER :: g = 9.8
1726	> REAL, PARAMETER :: Rd = 287.
1727	4435,4444c4171,4172
1728	<
1729	< !****MARS: as is done in MCD/pres0 with the MOLA topography :)
1730	<
1731	<
1732	< !! del_z = diff in surface topo, lo-res vs hi-res
1733	< !grid%em_ght_gc - grid%ht
1734	< !!* em_ght_gc: surface geopotential height from the GCM
1735	< !!* ht: hi-res altimetry
1736	< ! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
1737	<
1738	---
1739	> ! del_z = diff in surface topo, lo-res vs hi-res
1740	> ! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
1741	4448,4450d4175
1742	< !!
1743	< !!****MARS: 'ez_method' is 'we_have_tavgsfc', hard-coded as false
1744	< !!
1745	4459,4462d4183
1746	< !!****MARS .... here is what is done
1747	< !!****TODO:
1748	< !!****MARS .... toujours 0.608 ???
1749	< !!****MARS .... changer pour la température à 1 km ???
1750	4468,4478d4188
1751	< ! !****MARS .... check of the altimetry differences
1752	< ! print *,del_z, tv_sfc
1753	<
1754	<
1755	< !****MARS
1756	< !****MARS .... which temperature is used in the Laplace formula ?
1757	< !!****TODO: change the 220K value ??
1758	< !!****NB: pas d'influence énorme cependant de la valeur de T
1759	< psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * 220 ) )
1760	<
1761	<

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

Download in other formats:

Original Format