Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real_method2_1param1interp.F @ 2756

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

Line
1	!REAL:MODEL_LAYER:INITIALIZATION
2
3	#ifndef VERT_UNIT
4	! This MODULE holds the routines which are used to perform various initializations
5	! for the individual domains, specifically for the Eulerian, mass-based coordinate.
6
7	!-----------------------------------------------------------------------
8
9	!****MARS: modified May 2007
10
11
12	MODULE module_initialize
13
14	USE module_bc
15	USE module_configure
16	USE module_domain
17	USE module_io_domain
18	USE module_model_constants
19	USE module_state_description
20	USE module_timing
21	USE module_soil_pre
22	USE module_date_time
23	#ifdef DM_PARALLEL
24	USE module_dm
25	#endif
26
27	REAL , SAVE :: p_top_save
28	INTEGER :: internal_time_loop
29
30	CONTAINS
31
32	!-------------------------------------------------------------------
33
34	SUBROUTINE init_domain ( grid )
35
36	IMPLICIT NONE
37
38	! Input space and data. No gridded meteorological data has been stored, though.
39
40	! TYPE (domain), POINTER :: grid
41	TYPE (domain) :: grid
42
43	! Local data.
44
45	INTEGER :: dyn_opt
46	INTEGER :: idum1, idum2
47
48	CALL nl_get_dyn_opt ( 1, dyn_opt )
49
50	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
51
52	IF ( dyn_opt .eq. 1 &
53	.or. dyn_opt .eq. 2 &
54	.or. dyn_opt .eq. 3 &
55	) THEN
56	CALL init_domain_rk( grid &
57	!
58	#include "em_actual_new_args.inc"
59	!
60	)
61
62	ELSE
63	WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
64	CALL wrf_error_fatal ( 'ERROR-dyn_opt-wrong-in-namelist' )
65	ENDIF
66
67	END SUBROUTINE init_domain
68
69	!-------------------------------------------------------------------
70
71	SUBROUTINE init_domain_rk ( grid &
72	!
73	#include "em_dummy_new_args.inc"
74	!
75	)
76
77	USE module_optional_si_input
78	IMPLICIT NONE
79
80	! Input space and data. No gridded meteorological data has been stored, though.
81
82	! TYPE (domain), POINTER :: grid
83	TYPE (domain) :: grid
84
85	#include "em_dummy_new_decl.inc"
86
87	TYPE (grid_config_rec_type) :: config_flags
88
89	! Local domain indices and counters.
90
91	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
92	INTEGER :: loop , num_seaice_changes
93
94	INTEGER :: &
95	ids, ide, jds, jde, kds, kde, &
96	ims, ime, jms, jme, kms, kme, &
97	its, ite, jts, jte, kts, kte, &
98	ips, ipe, jps, jpe, kps, kpe, &
99	i, j, k
100	INTEGER :: ns
101
102	! Local data
103
104	INTEGER :: error
105	REAL :: p_surf, p_level
106	REAL :: cof1, cof2
107	REAL :: qvf , qvf1 , qvf2 , pd_surf
108	REAL :: p00 , t00 , a
109	REAL :: hold_znw
110	LOGICAL :: were_bad
111
112	LOGICAL :: stretch_grid, dry_sounding, debug
113	INTEGER IICOUNT
114
115	REAL :: p_top_requested , temp
116	INTEGER :: num_metgrid_levels
117	REAL , DIMENSION(max_eta) :: eta_levels
118	REAL :: max_dz
119
120	! INTEGER , PARAMETER :: nl_max = 1000
121	! REAL , DIMENSION(nl_max) :: grid%em_dn
122
123	integer::oops1,oops2
124
125	REAL :: zap_close_levels
126	INTEGER :: force_sfc_in_vinterp
127	INTEGER :: interp_type , lagrange_order
128	LOGICAL :: lowest_lev_from_sfc
129	LOGICAL :: we_have_tavgsfc
130
131	INTEGER :: lev500 , loop_count
132	REAL :: zl , zu , pl , pu , z500 , dz500 , tvsfc , dpmu
133
134	!-- Carsel and Parrish [1988]
135	REAL , DIMENSION(100) :: lqmi
136
137
138	!****MARS
139	INTEGER :: sizegcm, kold, knew,inew,jnew
140	REAL :: pa, indic, p1, p2, pn
141	REAL, ALLOCATABLE, DIMENSION (:,:,:) :: sig, ap, bp, box
142	!****MARS
143
144
145	#ifdef DM_PARALLEL
146	# include "em_data_calls.inc"
147	#endif
148
149	SELECT CASE ( model_data_order )
150	CASE ( DATA_ORDER_ZXY )
151	kds = grid%sd31 ; kde = grid%ed31 ;
152	ids = grid%sd32 ; ide = grid%ed32 ;
153	jds = grid%sd33 ; jde = grid%ed33 ;
154
155	kms = grid%sm31 ; kme = grid%em31 ;
156	ims = grid%sm32 ; ime = grid%em32 ;
157	jms = grid%sm33 ; jme = grid%em33 ;
158
159	kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
160	its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
161	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
162
163	CASE ( DATA_ORDER_XYZ )
164	ids = grid%sd31 ; ide = grid%ed31 ;
165	jds = grid%sd32 ; jde = grid%ed32 ;
166	kds = grid%sd33 ; kde = grid%ed33 ;
167
168	ims = grid%sm31 ; ime = grid%em31 ;
169	jms = grid%sm32 ; jme = grid%em32 ;
170	kms = grid%sm33 ; kme = grid%em33 ;
171
172	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
173	jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
174	kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
175
176	CASE ( DATA_ORDER_XZY )
177	ids = grid%sd31 ; ide = grid%ed31 ;
178	kds = grid%sd32 ; kde = grid%ed32 ;
179	jds = grid%sd33 ; jde = grid%ed33 ;
180
181	ims = grid%sm31 ; ime = grid%em31 ;
182	kms = grid%sm32 ; kme = grid%em32 ;
183	jms = grid%sm33 ; jme = grid%em33 ;
184
185	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
186	kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
187	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
188
189	END SELECT
190
191	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
192
193	! Check to see if the boundary conditions are set properly in the namelist file.
194	! This checks for sufficiency and redundancy.
195
196	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
197
198	! Some sort of "this is the first time" initialization. Who knows.
199
200	grid%step_number = 0
201	grid%itimestep=0
202
203	! Pull in the info in the namelist to compare it to the input data.
204
205	grid%real_data_init_type = model_config_rec%real_data_init_type
206
207	! To define the base state, we call a USER MODIFIED routine to set the three
208	! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
209	! and A (temperature difference, from 1000 mb to 300 mb, K).
210
211	CALL const_module_initialize ( p00 , t00 , a )
212
213	#if 0
214	!KLUDGE, this is for testing only
215	if ( flag_metgrid .eq. 1 ) then
216	read (20+grid%id) grid%em_ht_gc
217	read (20+grid%id) grid%em_xlat_gc
218	read (20+grid%id) grid%em_xlong_gc
219	read (20+grid%id) msft
220	read (20+grid%id) msfu
221	read (20+grid%id) msfv
222	read (20+grid%id) f
223	read (20+grid%id) e
224	read (20+grid%id) sina
225	read (20+grid%id) cosa
226	read (20+grid%id) grid%landmask
227	read (20+grid%id) grid%landusef
228	read (20+grid%id) grid%soilctop
229	read (20+grid%id) grid%soilcbot
230	read (20+grid%id) grid%vegcat
231	read (20+grid%id) grid%soilcat
232	read (20+grid%id) grid%em_albedo_gcm
233	read (20+grid%id) grid%em_therm_inert
234	else
235	write (20+grid%id) grid%em_ht
236	write (20+grid%id) grid%em_xlat
237	write (20+grid%id) grid%em_xlong
238	write (20+grid%id) msft
239	write (20+grid%id) msfu
240	write (20+grid%id) msfv
241	write (20+grid%id) f
242	write (20+grid%id) e
243	write (20+grid%id) sina
244	write (20+grid%id) cosa
245	write (20+grid%id) grid%landmask
246	write (20+grid%id) grid%landusef
247	write (20+grid%id) grid%soilctop
248	write (20+grid%id) grid%soilcbot
249	write (20+grid%id) grid%vegcat
250	write (20+grid%id) grid%soilcat
251	write (20+grid%id) grid%em_albedo_gcm
252	write (20+grid%id) grid%em_therm_inert
253	endif
254	#endif
255
256
257	! Is there any vertical interpolation to do? The "old" data comes in on the correct
258	! vertical locations already.
259
260	IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
261
262	! Variables that are named differently between SI and WPS.
263
264	DO j = jts, MIN(jte,jde-1)
265	DO i = its, MIN(ite,ide-1)
266	!!****MARS: tsk is surface temperature
267	grid%tsk(i,j) = grid%em_tsk_gc(i,j)
268	grid%tmn(i,j) = grid%em_tmn_gc(i,j)
269	grid%xlat(i,j) = grid%em_xlat_gc(i,j)
270	grid%xlong(i,j) = grid%em_xlong_gc(i,j)
271	grid%ht(i,j) = grid%em_ht_gc(i,j)
272	!!****MARS
273	!!un peu artificiel, mais u10 et v10 sont des bons intermediaires (facultatifs de plus)
274	grid%u10(i,j) = grid%em_albedo_gcm(i,j)
275	grid%v10(i,j) = grid%em_therm_inert(i,j)
276	include "../custom_static.inc"
277	!!****MARS
278	END DO
279	END DO
280
281
282	!!****MARS
283	!!
284	!! case with idealized topography
285	!!
286	!!CALL ideal_topo ( grid%ht , 2000., 6., &
287	!CALL ideal_topo ( grid%ht , 2000., 3., &
288	! ids , ide , jds , jde , kds , kde , &
289	! ims , ime , jms , jme , kms , kme , &
290	! its , ite , jts , jte , kts , kte )
291	!!****MARS
292	!****MARS
293	!!
294	!! idealized atmospheric state
295	!!
296	!!
297	!! adjust_heights ne sert pas :)
298	print *,config_flags%adjust_heights
299	IF (config_flags%adjust_heights) THEN
300	print *, 'setting uniform values and profiles'
301	print *, 'u'
302	print *, grid%em_u_gc(its+1,:,jts+1)
303	print *, 'v'
304	print *, grid%em_v_gc(its+1,:,jts+1)
305	print *, 't'
306	print *, grid%em_t_gc(its+1,:,jts+1)
307	print *, 'p'
308	print *, grid%em_rh_gc(its+1,:,jts+1)
309	print *, 'geop'
310	print *, grid%em_ght_gc(its+1,:,jts+1)
311	print *, 'albedo'
312	print *, grid%u10(its+1,jts+1)
313	print *, 'thermal inertia'
314	print *, grid%v10(its+1,jts+1)
315	print *, 'topography'
316	print *, grid%ht(its+1,jts+1)
317	print *, 'toposoil'
318	print *, grid%toposoil(its+1,jts+1)
319	print *, 'surface temperature'
320	print *, grid%tsk(its+1,jts+1)
321	print *, 'surface pressure'
322	print *, grid%psfc(its+1,jts+1)
323	print *, grid%em_psfc_gc(its+1,jts+1)
324
325	DO j = jts, MIN(jte,jde-1)
326	DO i = its, MIN(ite,ide-1)
327	grid%em_u_gc(i,:,j)=grid%em_u_gc(its+1,:,jts+1)
328	grid%em_v_gc(i,:,j)=grid%em_v_gc(its+1,:,jts+1)
329	grid%em_t_gc(i,:,j)=grid%em_t_gc(its+1,:,jts+1)
330	grid%em_rh_gc(i,:,j)=grid%em_rh_gc(its+1,:,jts+1)
331	grid%em_ght_gc(i,:,j) = grid%em_ght_gc(its+1,:,jts+1)
332	grid%u10(i,j) = grid%u10(its+1,jts+1)
333	grid%v10(i,j) = grid%v10(its+1,jts+1)
334	grid%ht(i,j) = grid%ht(its+1,jts+1)
335	grid%toposoil(i,j) = grid%toposoil(its+1,jts+1)
336	grid%tsk(i,j) = grid%tsk(its+1,jts+1)
337	grid%psfc(i,j) = grid%psfc(its+1,jts+1)
338	grid%em_psfc_gc(i,j) = grid%em_psfc_gc(its+1,jts+1)
339	ENDDO
340	ENDDO
341
342	!!
343	!!
344	print *, 'u divided by 10'
345	grid%em_u_gc=grid%em_u_gc/10.
346	grid%em_v_gc=grid%em_v_gc/10.
347	!!
348	!!
349
350
351	ENDIF
352	!!
353	!!
354	!! idealized atmospheric state
355	!!
356	!****MARS
357
358
359	!!!!caca caca
360	!print *, 'u multiplied by 3'
361	!grid%em_u_gc=grid%em_u_gc*3.
362	!grid%em_v_gc=grid%em_v_gc*3.
363	!!!!caca caca
364
365
366
367
368	! If we have any input low-res surface pressure, we store it.
369
370	!!****MARS
371	!!fix pour être certain d'être avec les bons flag
372	print *,flag_psfc
373	print *,flag_soilhgt
374	print *,flag_metgrid
375
376	flag_psfc=1
377	flag_soilhgt=1
378	flag_metgrid=1
379	!!**** TODO: trouver quand même pourquoi ça donne 0 :)
380	pa=999999.
381	!!****MARS
382
383	IF ( flag_psfc .EQ. 1 ) THEN
384	DO j = jts, MIN(jte,jde-1)
385	DO i = its, MIN(ite,ide-1)
386	grid%em_psfc_gc(i,j) = grid%psfc(i,j)
387	!!!****MARS: em_p_gc is only a way to count vertical levels in WPS :)
388	!!!****MARS: is filled here with real pressure levels
389	grid%em_p_gc(i,:,j) = grid%em_rh_gc(i,:,j)
390	!!!****MARS
391	grid%em_p_gc(i,1,j) = grid%psfc(i,j)
392	!!!****MARS
393	IF (pa .gt. grid%em_p_gc(i,1,j)) pa=grid%em_p_gc(i,1,j)
394	!!!****MARS
395	END DO
396	END DO
397	END IF
398	print *, 'found minimum pressure (Pa) :',pa
399
400
401	!!!!****MARS
402	!!!!****MARS
403	!!!! define new hybrid coordinate levels
404	!!!! with transition level between sigma and pressure
405	!!!! lower than input data
406	!
407	!
408	! !--get vertical size of the GCM input array
409	! sizegcm=SIZE(grid%em_rh_gc(1,:,1))
410	! ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
411	! ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
412	! ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
413	! !ALLOCATE(box(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
414	!
415	!
416	!
417	! !--define sigma levels,
418	! !--then derive new sigma levels, and new pressure levels
419	! DO j = jts, MIN(jte,jde-1)
420	! DO i = its, MIN(ite,ide-1)
421	!
422	! ! old sigma levels
423	! sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j)
424	!! sig(i,:,j)=grid%em_p_gc(20,:,20)/grid%em_psfc_gc(20,20)
425	! ! new pressure levels
426	! ! - pressure_new = ap_new + bp_new * ps_gcm
427	! ! - bp_new is converging more rapidly than bp
428	! ! ... while conserving the same structure near the surface
429	! !
430	! ! NB: grid%zap_close_levels ne sert pas dans vert_interp_old :)
431	! ! NB: peut donc servir pour préciser une constante reelle
432	! ! NB: qui permet de rehausser la zone de transition
433	! !
434	! bp(i,:,j)=sqrt(sqrt(exp(1.-1./(sig(i,:,j)**4))))
435	! ap(i,:,j)=paexp(-grid%zap_close_levels/10.)(sig(i,:,j)-bp(i,:,j))
436	! grid%em_rh_gc(i,:,j)=ap(i,:,j)+bp(i,:,j)*grid%em_psfc_gc(i,j)
437	!
438	! ! avoid extrapolation at the top
439	! ! -- the last level is thus unsignificant
440	! grid%em_p_gc(i,sizegcm,j)=grid%em_p_gc(i,sizegcm,j)/100.
441	!! grid%em_p_gc(i,sizegcm,j)=grid%em_p_gc(i,sizegcm,j)/10000.
442	!
443	! ENDDO
444	! ENDDO
445	!
446	!
447	!
448	!
449	! !-- check that the biggest differences are higher
450	! print *, 'sigma levels'
451	! print *, sig(its+1,:,jts+1)
452	! print *, 'old pressure levels'
453	! print *, grid%em_p_gc(its+1,:,jts+1)
454	! print *, 'new pressure levels'
455	! print *, grid%em_rh_gc(its+1,:,jts+1)
456	!
457	!!print *, 't_gc', SIZE(grid%em_t_gc(1,:,1))
458	!!print *, 'p_gc', SIZE(grid%em_p_gc(1,:,1))
459	!!print *, 't_2', SIZE(grid%em_t_2(1,:,1))
460	!!print *, 'rh_gc', SIZE(grid%em_rh_gc(1,:,1))
461	!
462	!
463	!!--------
464	!!-- interpolate on the new pressure levels
465	!!--------
466	!
467	! DO j = jts, MIN(jte,jde-1)
468	! DO i = its, MIN(ite,ide-1)
469	!
470	!DO knew = 1,sizegcm ! loop on each level of the new grid
471	!
472	! DO kold = 1,sizegcm-1 ! find the two enclosing levels
473	!
474	! ! indic becomes negative when the two levels are found
475	! indic=(grid%em_p_gc(i,kold,j)-grid%em_rh_gc(i,knew,j))&
476	! *(grid%em_p_gc(i,kold+1,j)-grid%em_rh_gc(i,knew,j))
477	!
478	! ! 1. the two levels are found - define p1,p2,pn and exit the loop
479	! IF (indic < 0.) THEN
480	! !IF ((i == its) .AND. (j == jts)) THEN !just a check
481	! ! print *, 'new level', grid%em_rh_gc(i,knew,j)
482	! ! print *, 'interp levels', grid%em_p_gc(i,kold,j), &
483	! ! grid%em_p_gc(i,kold+1,j)
484	! !ENDIF
485	! p1 = ALOG(grid%em_p_gc(i,kold,j))
486	! p2 = ALOG(grid%em_p_gc(i,kold+1,j))
487	! pn = ALOG(grid%em_rh_gc(i,knew,j))
488	! EXIT
489	!
490	! ! 2. must handle the case (usually close to the surface)
491	! ! of similar new/old levels - then exit with the right kold value
492	! ELSE IF (1-abs(grid%em_rh_gc(i,knew,j)/grid%em_p_gc(i,kold,j)) .lt. 1e-8) THEN
493	! !print *,grid%em_p_gc(i,kold,j),grid%em_rh_gc(i,knew,j)
494	! EXIT
495	! ELSE IF (1-abs(grid%em_rh_gc(i,knew,j)/grid%em_p_gc(i,kold+1,j)) .lt. 1e-8) THEN
496	! !print *,grid%em_p_gc(i,kold+1,j),grid%em_rh_gc(i,knew,j)
497	! kold=kold+1
498	! EXIT
499	!
500	! ! 3. continue looping if the two levels are not found ....
501	! ENDIF
502	! ENDDO
503	!
504	! ! this is an initialization, useful for case 2 (and erased just below if case 1)
505	! grid%em_t_2(i,knew,j)= grid%em_t_gc(i,kold,j)
506	! grid%em_u_2(i,knew,j)= grid%em_u_gc(i,kold,j)
507	! grid%em_v_2(i,knew,j)= grid%em_v_gc(i,kold,j)
508	!
509	! ! case 1: OK, in the previous loop, the two levels were found, and stored in p1 and p2
510	! ! ... thus interpolation can be performed
511	! IF (indic < 0.) THEN
512	! grid%em_t_2(i,knew,j)= ( grid%em_t_gc(i,kold,j) * ( p2 - pn ) + &
513	! grid%em_t_gc(i,kold+1,j) * ( pn - p1 ) ) / &
514	! ( p2 - p1 )
515	! grid%em_u_2(i,knew,j)= ( grid%em_u_gc(i,kold,j) * ( p2 - pn ) + &
516	! grid%em_u_gc(i,kold+1,j) * ( pn - p1 ) ) / &
517	! ( p2 - p1 )
518	! grid%em_v_2(i,knew,j)= ( grid%em_v_gc(i,kold,j) * ( p2 - pn ) + &
519	! grid%em_v_gc(i,kold+1,j) * ( pn - p1 ) ) / &
520	! ( p2 - p1 )
521	! ENDIF
522	!
523	!
524	!ENDDO
525	!
526	! ENDDO
527	! ENDDO
528	!grid%em_u_2(MIN(ite,ide-1)+1,:,:)=grid%em_u_2(MIN(ite,ide-1),:,:)
529	!grid%em_v_2(:,:,MIN(jte,jde-1)+1)=grid%em_v_2(:,:,MIN(jte,jde-1))
530	!!--------
531	!!-- end - interpolate on the new pressure levels
532	!!--------
533
534
535	! !-- interpolate on the new pressure levels
536	! CALL vert_interp_old ( grid%em_t_gc , & ! --- interpolate this field
537	! grid%em_p_gc, & ! --- with coordinates
538	! grid%em_t_2, & ! --- to obtain the new field
539	! grid%em_rh_gc, & ! --- on coordinates
540	! sizegcm, &
541	! 'T', & ! --- no staggering (will be done later)
542	! 2, & ! --- log p interpolation
543	! 1, & ! --- (0) lagrange_order
544	! .false., & ! --- (0) lowest_lev_from_sfc
545	! 0., & ! --- (0) zap_close_levels
546	! 0, & ! --- (0) force_sfc_in_vinterp
547	! ids , ide , jds , jde , kds , kde , &
548	! ims , ime , jms , jme , kms , kme , &
549	! its , ite , jts , jte , kts , kte )
550	! CALL vert_interp_old ( grid%em_u_gc , & ! --- interpolate this field
551	! grid%em_p_gc, & ! --- with coordinates
552	! grid%em_u_2 , & ! --- to obtain the new field
553	! grid%em_rh_gc, & ! --- on coordinates
554	! sizegcm, &
555	! 'U', & ! --- no staggering (will be done later)
556	! 2, & ! --- log p interpolation
557	! 1, & ! --- (0) lagrange_order
558	! .false., & ! --- (0) lowest_lev_from_sfc
559	! 0., & ! --- (0) zap_close_levels
560	! 0, & ! --- (0) force_sfc_in_vinterp
561	! ids , ide , jds , jde , kds , kde , &
562	! ims , ime , jms , jme , kms , kme , &
563	! its , ite , jts , jte , kts , kte )
564	! CALL vert_interp_old ( grid%em_v_gc , & ! --- interpolate this field
565	! grid%em_p_gc, & ! --- with coordinates
566	! grid%em_v_2 , & ! --- to obtain the new field
567	! grid%em_rh_gc, & ! --- on coordinates
568	! sizegcm, &
569	! 'V', & ! --- no staggering (will be done later)
570	! 2, & ! --- log p interpolation
571	! 1, & ! --- (0) lagrange_order
572	! .false., & ! --- (0) lowest_lev_from_sfc
573	! 0., & ! --- (0) zap_close_levels
574	! 0, & ! --- (0) force_sfc_in_vinterp
575	! !ids , ide , jds , jde , kds , sizegcm , &
576	! !ims , ime , jms , jme , kms , sizegcm , &
577	! !its , ite , jts , jte , kts , sizegcm )
578	! ids , ide , jds , jde , kds , kde , &
579	! ims , ime , jms , jme , kms , kme , &
580	! its , ite , jts , jte , kts , kte )
581	!
582	!
583	! !-- save the new field and the new pressure coordinates
584	! !-- these will be regarded now as the inputs from the GCM
585	! grid%em_t_gc=grid%em_t_2
586	! grid%em_t_2(:,:,:)=0.
587	! grid%em_u_gc=grid%em_u_2
588	! grid%em_u_2(:,:,:)=0.
589	! grid%em_v_gc=grid%em_v_2
590	! grid%em_v_2(:,:,:)=0.
591	! grid%em_p_gc=grid%em_rh_gc
592	! grid%em_rh_gc(:,:,:)=0.
593	!!!!****MARS
594	!!!****MARS
595
596
597
598
599	! If we have the low-resolution surface elevation, stick that in the
600	! "input" locations of the 3d height. We still have the "hi-res" topo
601	! stuck in the grid%em_ht array. The grid%landmask if test is required as some sources
602	! have ZERO elevation over water (thank you very much).
603
604	IF ( flag_soilhgt .EQ. 1) THEN
605	DO j = jts, MIN(jte,jde-1)
606	DO i = its, MIN(ite,ide-1)
607	IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
608	grid%em_ght_gc(i,1,j) = grid%toposoil(i,j)
609	grid%em_ht_gc(i,j)= grid%toposoil(i,j)
610	END IF
611	END DO
612	END DO
613	END IF
614
615	! Assign surface fields with original input values. If this is hybrid data,
616	! the values are not exactly representative. However - this is only for
617	! plotting purposes and such at the 0h of the forecast, so we are not all that
618	! worried.
619
620	!****MARS
621	! DO j = jts, min(jde-1,jte)
622	! DO i = its, min(ide,ite)
623	! grid%u10(i,j)=grid%em_u_gc(i,1,j)
624	! END DO
625	! END DO
626	!
627	! DO j = jts, min(jde,jte)
628	! DO i = its, min(ide-1,ite)
629	! grid%v10(i,j)=grid%em_v_gc(i,1,j)
630	! END DO
631	! END DO
632	!****MARS
633
634	DO j = jts, min(jde-1,jte)
635	DO i = its, min(ide-1,ite)
636	grid%t2(i,j)=grid%em_t_gc(i,1,j)
637	END DO
638	END DO
639
640	IF ( flag_qv .EQ. 1 ) THEN
641	DO j = jts, min(jde-1,jte)
642	DO i = its, min(ide-1,ite)
643	grid%q2(i,j)=grid%em_qv_gc(i,1,j)
644	END DO
645	END DO
646	END IF
647
648	! The number of vertical levels in the input data. There is no staggering for
649	! different variables.
650
651	num_metgrid_levels = grid%num_metgrid_levels
652
653	! The requested ptop for real data cases.
654
655	p_top_requested = grid%p_top_requested
656
657	! Compute the top pressure, grid%p_top. For isobaric data, this is just the
658	! top level. For the generalized vertical coordinate data, we find the
659	! max pressure on the top level. We have to be careful of two things:
660	! 1) the value has to be communicated, 2) the value can not increase
661	! at subsequent times from the initial value.
662
663	IF ( internal_time_loop .EQ. 1 ) THEN
664
665	CALL find_p_top ( grid%em_p_gc , grid%p_top , &
666	ids , ide , jds , jde , 1 , num_metgrid_levels , &
667	ims , ime , jms , jme , 1 , num_metgrid_levels , &
668	its , ite , jts , jte , 1 , num_metgrid_levels )
669
670	!! ^---- equivalent to:
671	!!grid%ptop=MINVAL(grid%em_p_gc(:,:,:))
672
673
674
675	!!!!obsolete
676	!print *,'ptop GCM',grid%em_rh_gc(2,1,2)
677	!IF (grid%em_rh_gc(2,1,2) == 0) THEN
678	! print *,'ptop cannot be 0'
679	! stop
680	!ENDIF
681	!grid%p_top=grid%em_rh_gc(2,1,2)
682	!!!!obsolete
683
684
685	#ifdef DM_PARALLEL
686	grid%p_top = wrf_dm_max_real ( grid%p_top )
687	#endif
688
689	! Compare the requested grid%p_top with the value available from the input data.
690
691	print *,'p_top_requested = ',p_top_requested
692	print *,'allowable grid%p_top in data = ',grid%p_top
693	IF ( p_top_requested .LT. grid%p_top ) THEN
694	CALL wrf_error_fatal ( 'p_top_requested < grid%p_top possible from data' )
695	END IF
696
697	! The grid%p_top valus is the max of what is available from the data and the
698	! requested value. We have already compared <, so grid%p_top is directly set to
699	! the value in the namelist.
700
701	grid%p_top = p_top_requested
702
703	! For subsequent times, we have to remember what the grid%p_top for the first
704	! time was. Why? If we have a generalized vert coordinate, the grid%p_top value
705	! could fluctuate.
706
707	p_top_save = grid%p_top
708
709	ELSE
710	CALL find_p_top ( grid%em_p_gc , grid%p_top , &
711	ids , ide , jds , jde , 1 , num_metgrid_levels , &
712	ims , ime , jms , jme , 1 , num_metgrid_levels , &
713	its , ite , jts , jte , 1 , num_metgrid_levels )
714
715	#ifdef DM_PARALLEL
716	grid%p_top = wrf_dm_max_real ( grid%p_top )
717	#endif
718	IF ( grid%p_top .GT. p_top_save ) THEN
719	print *,'grid%p_top from last time period = ',p_top_save
720	print *,'grid%p_top from this time period = ',grid%p_top
721	CALL wrf_error_fatal ( 'grid%p_top > previous value' )
722	END IF
723	grid%p_top = p_top_save
724	ENDIF
725
726
727	!****MARS
728	!****MARS
729	print *,'ptop GCM', grid%p_top
730	print *,'sample: pressure at its jts'
731	print *,grid%em_p_gc(its,:,jts)
732	!****MARS
733	!****MARS
734
735
736
737	!****MARS: useless
738	!****MARS:
739	! ! Get the monthly values interpolated to the current date for the traditional monthly
740	! ! fields of green-ness fraction and background albedo.
741	!
742	! CALL monthly_interp_to_date ( grid%em_greenfrac , current_date , grid%vegfra , &
743	! ids , ide , jds , jde , kds , kde , &
744	! ims , ime , jms , jme , kms , kme , &
745	! its , ite , jts , jte , kts , kte )
746	!
747	! CALL monthly_interp_to_date ( grid%em_albedo12m , current_date , grid%albbck , &
748	! ids , ide , jds , jde , kds , kde , &
749	! ims , ime , jms , jme , kms , kme , &
750	! its , ite , jts , jte , kts , kte )
751	!
752	! ! Get the min/max of each i,j for the monthly green-ness fraction.
753	!
754	! CALL monthly_min_max ( grid%em_greenfrac , grid%shdmin , grid%shdmax , &
755	! ids , ide , jds , jde , kds , kde , &
756	! ims , ime , jms , jme , kms , kme , &
757	! its , ite , jts , jte , kts , kte )
758	!
759	! ! The model expects the green-ness values in percent, not fraction.
760	!
761	! DO j = jts, MIN(jte,jde-1)
762	! DO i = its, MIN(ite,ide-1)
763	! grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
764	! grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
765	! grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
766	! END DO
767	! END DO
768	!
769	! ! The model expects the albedo fields as a fraction, not a percent. Set the
770	! ! water values to 8%.
771	!
772	! DO j = jts, MIN(jte,jde-1)
773	! DO i = its, MIN(ite,ide-1)
774	! grid%albbck(i,j) = grid%albbck(i,j) / 100.
775	! grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
776	! IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
777	! grid%albbck(i,j) = 0.08
778	! grid%snoalb(i,j) = 0.08
779	! END IF
780	! END DO
781	! END DO
782	!!****MARS:
783	!!****MARS: useless
784
785
786
787
788	!!****MARS:
789	! ! Compute the mixing ratio from the input relative humidity.
790	!
791	! IF ( flag_qv .NE. 1 ) THEN
792	! CALL rh_to_mxrat (grid%em_rh_gc, grid%em_t_gc, grid%em_p_gc, grid%em_qv_gc , .TRUE. , &
793	! ids , ide , jds , jde , 1 , num_metgrid_levels , &
794	! ims , ime , jms , jme , 1 , num_metgrid_levels , &
795	! its , ite , jts , jte , 1 , num_metgrid_levels )
796	! END IF
797	!!****MARS:
798	!!grid%em_rh_gc are GCM equivalent eta_levels
799	!!****MARS
800
801
802
803
804	! Two ways to get the surface pressure. 1) If we have the low-res input surface
805	! pressure and the low-res topography, then we can do a simple hydrostatic
806	! relation. 2) Otherwise we compute the surface pressure from the sea-level
807	! pressure.
808	! Note that on output, grid%em_psfc is now hi-res. The low-res surface pressure and
809	! elevation are grid%em_psfc_gc and grid%em_ht_gc (same as grid%em_ght_gc(k=1)).
810
811	!!****MARS: switch off this option
812	!!****MARS: --> cf sfcprs2 and geopotential function at 500mb
813	! IF ( config_flags%adjust_heights ) THEN
814	! we_have_tavgsfc = ( flag_tavgsfc == 1 )
815	! ELSE
816	! we_have_tavgsfc = .FALSE.
817	! END IF
818	!****MARS:
819	we_have_tavgsfc = .FALSE.
820
821
822
823	!****MARS: hi-res psfc is done if the flag 'sfcp_to_sfcp' is active (recommended)
824	IF ( ( flag_psfc .EQ. 1 ) .AND. ( flag_soilhgt .EQ. 1 ) .AND. &
825	( config_flags%sfcp_to_sfcp ) ) THEN
826	print *,'compute psfc from hi-res topography'
827	CALL sfcprs2(grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_psfc_gc, grid%ht, &
828	grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, &
829	ids , ide , jds , jde , 1 , num_metgrid_levels , &
830	ims , ime , jms , jme , 1 , num_metgrid_levels , &
831	its , ite , jts , jte , 1 , num_metgrid_levels )
832	!****MARS: here, in reality, grid%em_p_gc is not used
833
834	!****MARS: no sea-level pressure inputs possible
835	! ELSE
836	! CALL sfcprs (grid%em_t_gc, grid%em_qv_gc, grid%em_ght_gc, grid%em_pslv_gc, grid%ht, &
837	! grid%em_tavgsfc, grid%em_p_gc, grid%psfc, we_have_tavgsfc, &
838	! ids , ide , jds , jde , 1 , num_metgrid_levels , &
839	! ims , ime , jms , jme , 1 , num_metgrid_levels , &
840	! its , ite , jts , jte , 1 , num_metgrid_levels )
841	!****MARS: no sea-level pressure inputs possible
842
843
844	! If we have no input surface pressure, we'd better stick something in there.
845
846	IF ( flag_psfc .NE. 1 ) THEN
847	DO j = jts, MIN(jte,jde-1)
848	DO i = its, MIN(ite,ide-1)
849	grid%em_psfc_gc(i,j) = grid%psfc(i,j)
850	grid%em_p_gc(i,1,j) = grid%psfc(i,j)
851	END DO
852	END DO
853	END IF
854
855	END IF
856
857
858	!!!****MARS:
859	!!!****MARS: old stuff
860	!!! grid%em_p_gc is needed ... so it is computed from eta_gcm
861	!
862	!print *,'computing pressure levels for input data...'
863	!
864	! !! pressure is computed from eta_gcm and hi-res topography
865	! DO j = jts, MIN(jte,jde-1)
866	! DO i = its, MIN(ite,ide-1)
867	!!!psfc ou em_psfc_gc ? em_psfc_gc, sinon c'est faux et déclenche instabilités
868	!grid%em_p_gc(i,:,j)=grid%em_rh_gc(i,:,j)*(grid%em_psfc_gc(i,j)-grid%em_rh_gc(2,1,2))+grid%em_rh_gc(2,1,2)
869	!grid%em_p_gc(i,1,j)=grid%em_psfc_gc(i,j)
870	!
871	!
872	! END DO
873	! END DO
874	!!
875	!!****MARS:
876
877
878
879	!! Integrate the mixing ratio to get the vapor pressure.
880	!
881	!CALL integ_moist ( grid%em_qv_gc , grid%em_p_gc , grid%em_pd_gc , grid%em_t_gc , grid%em_ght_gc , grid%em_intq_gc , &
882	! ids , ide , jds , jde , 1 , num_metgrid_levels , &
883	! ims , ime , jms , jme , 1 , num_metgrid_levels , &
884	! its , ite , jts , jte , 1 , num_metgrid_levels )
885
886
887	!!****MARS
888	!!****MARS
889	!! and now, convert the GCM sigma levels into WRF sigma levels using hi-res surface pressure
890	!!DO j = jts , MIN ( jde-1 , jte )
891	!!DO i = its , MIN (ide-1 , ite )
892	!!
893	!! grid%em_pd_gc(i,:,j)=ap(i,:,j)+bp(i,:,j)*grid%psfc(i,j)
894	!!
895	!!END DO
896	!!END DO
897
898
899	!--get vertical size of the GCM input array and allocate new stuff
900	sizegcm=SIZE(grid%em_rh_gc(1,:,1))
901	ALLOCATE(sig(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
902	!ALLOCATE(ap(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
903	ALLOCATE(bp(MIN(ite,ide-1)-its+1,sizegcm, MIN(jte,jde-1)-jts+1))
904
905	DO j = jts , MIN ( jde-1 , jte )
906	DO i = its , MIN (ide-1 , ite )
907
908	!!! Define old sigma levels for each column
909	sig(i,:,j)=grid%em_p_gc(i,:,j)/grid%em_psfc_gc(i,j)
910
911	!!! Compute new sigma levels from old sigma levels with GCM (low-res) and WRF (hi-res) surface pressure
912	!!! (dimlevs,sigma_gcm, ps_gcm, ps_hr, sigma_hr)
913	CALL build_sigma_hr(sizegcm,sig(i,:,j),grid%em_psfc_gc(i,j),grid%psfc(i,j),bp(i,:,j))
914
915	!!! Calculate new pressure levels
916	grid%em_pd_gc(i,:,j)=bp(i,:,j)*grid%psfc(i,j)
917
918	END DO
919	END DO
920
921
922
923
924	!!****MARS
925	!grid%em_pd_gc=grid%em_p_gc
926	!!****MARS
927
928
929
930
931	!! Compute the difference between the dry, total surface pressure (input) and the
932	!! dry top pressure (constant).
933	!
934	!CALL p_dts ( grid%em_mu0 , grid%em_intq_gc , grid%psfc , grid%p_top , &
935	! ids , ide , jds , jde , 1 , num_metgrid_levels , &
936	! ims , ime , jms , jme , 1 , num_metgrid_levels , &
937	! its , ite , jts , jte , 1 , num_metgrid_levels )
938
939	!!****MARS
940	DO j = jts , MIN ( jde-1 , jte )
941	DO i = its , MIN (ide-1 , ite )
942
943	grid%em_mu0(i,j) = grid%psfc(i,j) - grid%p_top
944
945	END DO
946	END DO
947	!!****MARS
948
949
950	!! Compute the dry, hydrostatic surface pressure.
951	!
952	!CALL p_dhs ( grid%em_pdhs , grid%ht , p00 , t00 , a , &
953	! ids , ide , jds , jde , kds , kde , &
954	! ims , ime , jms , jme , kms , kme , &
955	! its , ite , jts , jte , kts , kte )
956	!!****MARS: voir remarques dans la routine
957	!!****MARS: dry hydrostatic pressure comes from the GCM ...
958	! DO j = jts , MIN ( jde-1 , jte )
959	! DO i = its , MIN (ide-1 , ite )
960	! grid%em_pdhs(i,j) = grid%psfc(i,j)
961	! END DO
962	! END DO
963	!!****MARS: em_pdhs ne sert qu'ici !
964
965
966	! Compute the eta levels if not defined already.
967
968	!!TODO: pb when ptop<1Pa
969
970	IF ( grid%em_znw(1) .NE. 1.0 ) THEN
971
972	eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
973	max_dz = model_config_rec%max_dz
974
975	!!****MARS
976	IF (grid%force_sfc_in_vinterp == 0) grid%force_sfc_in_vinterp = 8
977	!!default choice
978	!!****MARS
979
980	CALL compute_eta ( grid%em_znw , &
981	eta_levels , max_eta , max_dz , &
982	grid%force_sfc_in_vinterp, & !!ne sert pas par ailleurs
983	grid%p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , &
984	ids , ide , jds , jde , kds , kde , &
985	ims , ime , jms , jme , kms , kme , &
986	its , ite , jts , jte , kts , kte )
987	END IF
988
989	! The input field is temperature, we want potential temp.
990	!****MARS: here em_p_gc is really needed !
991	CALL t_to_theta ( grid%em_t_gc , grid%em_p_gc , p00 , &
992	ids , ide , jds , jde , 1 , num_metgrid_levels , &
993	ims , ime , jms , jme , 1 , num_metgrid_levels , &
994	its , ite , jts , jte , 1 , num_metgrid_levels )
995
996
997
998	! On the eta surfaces, compute the dry pressure = mu eta, stored in
999	! grid%em_pb, since it is a pressure, and we don't need another kms:kme 3d
1000	! array floating around. The grid%em_pb array is re-computed as the base pressure
1001	! later after the vertical interpolations are complete.
1002	CALL p_dry ( grid%em_mu0 , grid%em_znw , grid%p_top , grid%em_pb , &
1003	ids , ide , jds , jde , kds , kde , &
1004	ims , ime , jms , jme , kms , kme , &
1005	its , ite , jts , jte , kts , kte )
1006
1007
1008	!****MARS
1009	!****MARS: old stuff
1010	!****MARS
1011	!!! and now eta levels from the GCM are computed with the WRF ptop and GCM psfc
1012	!!! and em_pb is filled with WRF eta levels to prepare interpolation
1013	!print *,'computing eta levels for input data...'
1014	! DO j = jts, MIN(jte,jde-1)
1015	! DO i = its, MIN(ite,ide-1)
1016	!
1017	!!grid%em_psfc_gc: pb en haut!!!!
1018	!!!!valeurs plus grandes que 1 et extrapolation
1019	!!grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%psfc(i,j)-grid%p_top)
1020	!!!!utile si l'on est proche de la surface, mais pb plus haut !
1021	!grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%em_psfc_gc(i,j)-grid%p_top)
1022	!grid%em_pb(i,:,j)=grid%em_znw(:)
1023	!
1024	!!
1025	!!!!manage negative values
1026	!!DO k=1,num_metgrid_levels
1027	!! grid%em_p_gc(i,k,j)=MAX(0.,grid%em_p_gc(i,k,j))
1028	!!END DO
1029	!!
1030	!
1031	! END DO
1032	! END DO
1033	!!
1034	!!print *,'sample: eta GCM at its jts'
1035	!!print *,grid%em_p_gc(its,:,jts)
1036	!!print *,'sample: eta WRF at its jts'
1037	!!print *,grid%em_pb(its,:,jts)
1038	!!
1039	!!print *,grid%em_p_gc(:,2,:)
1040	!!print *, 'yeah yeah'
1041	!!grid%em_pd_gc(:,:,:)=grid%em_p_gc(:,:,:)
1042	!!
1043	!****MARS
1044	!****MARS: old stuff
1045	!****MARS
1046
1047
1048
1049
1050	! All of the vertical interpolations are done in dry-pressure space. The
1051	! input data has had the moisture removed (grid%em_pd_gc). The target levels (grid%em_pb)
1052	! had the vapor pressure removed from the surface pressure, then they were
1053	! scaled by the eta levels.
1054
1055	interp_type = grid%interp_type
1056	lagrange_order = grid%lagrange_order
1057	lowest_lev_from_sfc = grid%lowest_lev_from_sfc
1058	zap_close_levels = grid%zap_close_levels
1059	force_sfc_in_vinterp = grid%force_sfc_in_vinterp
1060
1061	!!****MARS: normalement c'est vert_interp
1062	!!****MARS: mais les résultats sont trop discontinus > retour à une
1063	!!****MARS: interpolation plus classique
1064	CALL vert_interp_old ( grid%em_qv_gc , grid%em_pd_gc , moist(:,:,:,P_QV) , grid%em_pb , &
1065	num_metgrid_levels , 'Q' , &
1066	interp_type , lagrange_order , lowest_lev_from_sfc , &
1067	zap_close_levels , force_sfc_in_vinterp , &
1068	ids , ide , jds , jde , kds , kde , &
1069	ims , ime , jms , jme , kms , kme , &
1070	its , ite , jts , jte , kts , kte )
1071
1072	!!****MARS: normalement c'est vert_interp
1073	CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2 , grid%em_pb , &
1074	num_metgrid_levels , 'T' , &
1075	interp_type , lagrange_order , lowest_lev_from_sfc , &
1076	zap_close_levels , force_sfc_in_vinterp , &
1077	ids , ide , jds , jde , kds , kde , &
1078	ims , ime , jms , jme , kms , kme , &
1079	its , ite , jts , jte , kts , kte )
1080
1081
1082	#if 0
1083	! Uncomment the Registry entries to activate these. This adds
1084	! noticeably to the allocated space for the model.
1085
1086	IF ( flag_qr .EQ. 1 ) THEN
1087	DO im = PARAM_FIRST_SCALAR, num_3d_m
1088	IF ( im .EQ. P_QR ) THEN
1089	CALL vert_interp_old ( qr_gc , grid%em_pd_gc , moist(:,:,:,P_QR) , grid%em_pb , &
1090	num_metgrid_levels , 'Q' , &
1091	interp_type , lagrange_order , lowest_lev_from_sfc , &
1092	zap_close_levels , force_sfc_in_vinterp , &
1093	ids , ide , jds , jde , kds , kde , &
1094	ims , ime , jms , jme , kms , kme , &
1095	its , ite , jts , jte , kts , kte )
1096	END IF
1097	END DO
1098	END IF
1099
1100	IF ( flag_qc .EQ. 1 ) THEN
1101	DO im = PARAM_FIRST_SCALAR, num_3d_m
1102	IF ( im .EQ. P_QC ) THEN
1103	CALL vert_interp_old ( qc_gc , grid%em_pd_gc , moist(:,:,:,P_QC) , grid%em_pb , &
1104	num_metgrid_levels , 'Q' , &
1105	interp_type , lagrange_order , lowest_lev_from_sfc , &
1106	zap_close_levels , force_sfc_in_vinterp , &
1107	ids , ide , jds , jde , kds , kde , &
1108	ims , ime , jms , jme , kms , kme , &
1109	its , ite , jts , jte , kts , kte )
1110	END IF
1111	END DO
1112	END IF
1113
1114	IF ( flag_qi .EQ. 1 ) THEN
1115	DO im = PARAM_FIRST_SCALAR, num_3d_m
1116	IF ( im .EQ. P_QI ) THEN
1117	CALL vert_interp_old ( qi_gc , grid%em_pd_gc , moist(:,:,:,P_QI) , grid%em_pb , &
1118	num_metgrid_levels , 'Q' , &
1119	interp_type , lagrange_order , lowest_lev_from_sfc , &
1120	zap_close_levels , force_sfc_in_vinterp , &
1121	ids , ide , jds , jde , kds , kde , &
1122	ims , ime , jms , jme , kms , kme , &
1123	its , ite , jts , jte , kts , kte )
1124	END IF
1125	END DO
1126	END IF
1127
1128	IF ( flag_qs .EQ. 1 ) THEN
1129	DO im = PARAM_FIRST_SCALAR, num_3d_m
1130	IF ( im .EQ. P_QS ) THEN
1131	CALL vert_interp_old ( qs_gc , grid%em_pd_gc , moist(:,:,:,P_QS) , grid%em_pb , &
1132	num_metgrid_levels , 'Q' , &
1133	interp_type , lagrange_order , lowest_lev_from_sfc , &
1134	zap_close_levels , force_sfc_in_vinterp , &
1135	ids , ide , jds , jde , kds , kde , &
1136	ims , ime , jms , jme , kms , kme , &
1137	its , ite , jts , jte , kts , kte )
1138	END IF
1139	END DO
1140	END IF
1141
1142	IF ( flag_qg .EQ. 1 ) THEN
1143	DO im = PARAM_FIRST_SCALAR, num_3d_m
1144	IF ( im .EQ. P_QG ) THEN
1145	CALL vert_interp_old ( qg_gc , grid%em_pd_gc , moist(:,:,:,P_QG) , grid%em_pb , &
1146	num_metgrid_levels , 'Q' , &
1147	interp_type , lagrange_order , lowest_lev_from_sfc , &
1148	zap_close_levels , force_sfc_in_vinterp , &
1149	ids , ide , jds , jde , kds , kde , &
1150	ims , ime , jms , jme , kms , kme , &
1151	its , ite , jts , jte , kts , kte )
1152	END IF
1153	END DO
1154	END IF
1155	#endif
1156
1157	#ifdef DM_PARALLEL
1158	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
1159
1160	! For the U and V vertical interpolation, we need the pressure defined
1161	! at both the locations for the horizontal momentum, which we get by
1162	! averaging two pressure values (i and i-1 for U, j and j-1 for V). The
1163	! pressure field on input (grid%em_pd_gc) and the pressure of the new coordinate
1164	! (grid%em_pb) are both communicated with an 8 stencil.
1165
1166	# include "HALO_EM_VINTERP_UV_1.inc"
1167	#endif
1168
1169	!!****MARS: normalement c'est vert_interp
1170	CALL vert_interp_old ( grid%em_u_gc , grid%em_pd_gc , grid%em_u_2, grid%em_pb , &
1171	num_metgrid_levels , 'U' , &
1172	interp_type , lagrange_order , lowest_lev_from_sfc , &
1173	zap_close_levels , force_sfc_in_vinterp , &
1174	ids , ide , jds , jde , kds , kde , &
1175	ims , ime , jms , jme , kms , kme , &
1176	its , ite , jts , jte , kts , kte )
1177	!!****MARS: normalement c'est vert_interp
1178	CALL vert_interp_old ( grid%em_v_gc , grid%em_pd_gc , grid%em_v_2, grid%em_pb , &
1179	num_metgrid_levels , 'V' , &
1180	interp_type , lagrange_order , lowest_lev_from_sfc , &
1181	zap_close_levels , force_sfc_in_vinterp , &
1182	ids , ide , jds , jde , kds , kde , &
1183	ims , ime , jms , jme , kms , kme , &
1184	its , ite , jts , jte , kts , kte )
1185
1186
1187	!!****MARS
1188	!!****MARS
1189	!!
1190	!! old obsolete method
1191	!! -------------------
1192	!!
1193	!!! and now eta levels from the GCM are computed with the WRF ptop and GCM psfc
1194	!!! and em_pb is filled with WRF eta levels to prepare interpolation
1195	!print *,'computing eta levels for input data...'
1196	!
1197	! DO j = jts, MIN(jte,jde-1)
1198	! DO i = its, MIN(ite,ide-1)
1199	!
1200	!! grid%em_psfc_gc: pb en haut!!!!
1201	!!!!valeurs plus grandes que 1 et extrapolation
1202	!! grid%em_p_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%psfc(i,j)-grid%p_top)
1203	!!!!utile si l'on est proche de la surface, mais pb plus haut !
1204	!grid%em_pd_gc(i,:,j)=(grid%em_p_gc(i,:,j)-grid%p_top)/(grid%em_psfc_gc(i,j)-grid%p_top)
1205	!grid%em_pb(i,:,j)=grid%em_znw(:)
1206	!
1207	!!
1208	!!!!manage negative values
1209	!!DO k=1,num_metgrid_levels
1210	!! grid%em_p_gc(i,k,j)=MAX(0.,grid%em_p_gc(i,k,j))
1211	!!END DO
1212	!!
1213	!
1214	! END DO
1215	! END DO
1216	!
1217	!print *,'sample: eta GCM at its jts'
1218	!print *,grid%em_pd_gc(its,:,jts)
1219	!print *,'sample: eta WRF at its jts'
1220	!print *,grid%em_pb(its,:,jts)
1221	!!
1222	!!!****MARS
1223	!
1224	!
1225	!
1226	!!!!****MARS
1227	!!!!
1228	!!!!grid%force_sfc_in_vinterp ne sert pas dans vert_interp_old :)
1229	!!!!peut donc servir pour préciser le nombre de niveaux
1230	!!!!pris à partir de l'interpolation eta
1231	!
1232	!IF (grid%force_sfc_in_vinterp .NE. 0) THEN
1233	!
1234	! !!!save in an array that is now unused
1235	! !!!the previously performed pressure interpolation
1236	! grid%em_qv_gc(:,:,:)=grid%em_t_2(:,:,:)
1237	!
1238	!
1239	! !!!perform interpolation on eta levels
1240	! print *, 'interpolate on eta levels for near-surface fields'
1241	! CALL vert_interp_old ( grid%em_t_gc , grid%em_pd_gc , grid%em_t_2, grid%em_pb , &
1242	! num_metgrid_levels , 'T' , &
1243	! interp_type , lagrange_order , lowest_lev_from_sfc ,&
1244	! zap_close_levels , force_sfc_in_vinterp , &
1245	! ids , ide , jds , jde , kds , kde , &
1246	! ims , ime , jms , jme , kms , kme , &
1247	! its , ite , jts , jte , kts , kte )
1248	!
1249	! !!!take the first layers from the eta interpolation
1250	! print *, 'the first ', &
1251	! grid%force_sfc_in_vinterp, &
1252	! 'layers will be taken from eta interpolation'
1253	! grid%em_qv_gc(:,1:grid%force_sfc_in_vinterp,:)=grid%em_t_2(:,1:grid%force_sfc_in_vinterp,:)
1254	!
1255	! !!!fix the possible little discontinuity at the limit
1256	! !!!between the two interpolation methods
1257	! grid%em_qv_gc(:,grid%force_sfc_in_vinterp+1,:)= &
1258	! 0.5*(grid%em_t_2(:,grid%force_sfc_in_vinterp,:) + & !!eta interpolation below
1259	! grid%em_qv_gc(:,grid%force_sfc_in_vinterp+2,:)) !!pressure interpolation above
1260	!
1261	!
1262	! !!!assign the final result in t_2
1263	! grid%em_t_2(:,:,:)=grid%em_qv_gc(:,:,:)
1264	! grid%em_qv_gc(:,:,:)=0.
1265	!
1266	!
1267	!ELSE
1268	!
1269	!
1270	!ENDIF
1271	!!****MARS
1272	!!****MARS
1273
1274
1275	END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
1276
1277
1278
1279	!****MARS: no need
1280	! ! Protect against bad grid%em_tsk values over water by supplying grid%sst (if it is
1281	! ! available, and if the grid%sst is reasonable).
1282	!
1283	! DO j = jts, MIN(jde-1,jte)
1284	! DO i = its, MIN(ide-1,ite)
1285	! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1286	! ( grid%sst(i,j) .GT. 200. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
1287	! grid%tsk(i,j) = grid%sst(i,j)
1288	! ENDIF
1289	! END DO
1290	! END DO
1291	!
1292	! ! Save the grid%em_tsk field for later use in the sea ice surface temperature
1293	! ! for the Noah LSM scheme.
1294	!
1295	! DO j = jts, MIN(jte,jde-1)
1296	! DO i = its, MIN(ite,ide-1)
1297	! grid%tsk_save(i,j) = grid%tsk(i,j)
1298	! END DO
1299	! END DO
1300	!
1301	!!****MARS: no need
1302	! ! Take the data from the input file and store it in the variables that
1303	! ! use the WRF naming and ordering conventions.
1304	!
1305	! DO j = jts, MIN(jte,jde-1)
1306	! DO i = its, MIN(ite,ide-1)
1307	! IF ( grid%snow(i,j) .GE. 10. ) then
1308	! grid%snowc(i,j) = 1.
1309	! ELSE
1310	! grid%snowc(i,j) = 0.0
1311	! END IF
1312	! END DO
1313	! END DO
1314	!
1315	! ! Set flag integers for presence of snowh and soilw fields
1316	!
1317	! grid%ifndsnowh = flag_snowh
1318	! IF (num_sw_levels_input .GE. 1) THEN
1319	! grid%ifndsoilw = 1
1320	! ELSE
1321	! grid%ifndsoilw = 0
1322	! END IF
1323	!
1324	!****MARS: no need
1325	! ! We require input data for the various LSM schemes.
1326	!
1327	! enough_data : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1328	!
1329	! CASE (LSMSCHEME)
1330	! IF ( num_st_levels_input .LT. 2 ) THEN
1331	! CALL wrf_error_fatal ( 'Not enough soil temperature data for Noah LSM scheme.')
1332	! END IF
1333	!
1334	! CASE (RUCLSMSCHEME)
1335	! IF ( num_st_levels_input .LT. 2 ) THEN
1336	! CALL wrf_error_fatal ( 'Not enough soil temperature data for RUC LSM scheme.')
1337	! END IF
1338	!
1339	! END SELECT enough_data
1340	!
1341	! ! For sf_surface_physics = 1, we want to use close to a 30 cm value
1342	! ! for the bottom level of the soil temps.
1343	!
1344	! fix_bottom_level_for_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1345	!
1346	! CASE (SLABSCHEME)
1347	! IF ( flag_tavgsfc .EQ. 1 ) THEN
1348	! DO j = jts , MIN(jde-1,jte)
1349	! DO i = its , MIN(ide-1,ite)
1350	! grid%tmn(i,j) = grid%em_tavgsfc(i,j)
1351	! END DO
1352	! END DO
1353	! ELSE IF ( flag_st010040 .EQ. 1 ) THEN
1354	! DO j = jts , MIN(jde-1,jte)
1355	! DO i = its , MIN(ide-1,ite)
1356	! grid%tmn(i,j) = grid%st010040(i,j)
1357	! END DO
1358	! END DO
1359	! ELSE IF ( flag_st000010 .EQ. 1 ) THEN
1360	! DO j = jts , MIN(jde-1,jte)
1361	! DO i = its , MIN(ide-1,ite)
1362	! grid%tmn(i,j) = grid%st000010(i,j)
1363	! END DO
1364	! END DO
1365	! ELSE IF ( flag_soilt020 .EQ. 1 ) THEN
1366	! DO j = jts , MIN(jde-1,jte)
1367	! DO i = its , MIN(ide-1,ite)
1368	! grid%tmn(i,j) = grid%soilt020(i,j)
1369	! END DO
1370	! END DO
1371	! ELSE IF ( flag_st007028 .EQ. 1 ) THEN
1372	! DO j = jts , MIN(jde-1,jte)
1373	! DO i = its , MIN(ide-1,ite)
1374	! grid%tmn(i,j) = grid%st007028(i,j)
1375	! END DO
1376	! END DO
1377	! ELSE
1378	! CALL wrf_debug ( 0 , 'No 10-40 cm, 0-10 cm, 7-28, or 20 cm soil temperature data for grid%em_tmn')
1379	! CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
1380	! END IF
1381	!
1382	! CASE (LSMSCHEME)
1383	!
1384	! CASE (RUCLSMSCHEME)
1385	!
1386	! END SELECT fix_bottom_level_for_temp
1387	!
1388	! ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is
1389	! ! is for the 5-layer scheme.
1390	!
1391	! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1392	! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1393	! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1394	! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1395	! CALL nl_get_isice ( grid%id , grid%isice )
1396	! CALL nl_get_iswater ( grid%id , grid%iswater )
1397	! CALL adjust_for_seaice_pre ( grid%xice , grid%landmask , grid%tsk , grid%ivgtyp , grid%vegcat , grid%lu_index , &
1398	! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , grid%soilctop , &
1399	! grid%soilcbot , grid%tmn , &
1400	! grid%seaice_threshold , &
1401	! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1402	! grid%iswater , grid%isice , &
1403	! model_config_rec%sf_surface_physics(grid%id) , &
1404	! ids , ide , jds , jde , kds , kde , &
1405	! ims , ime , jms , jme , kms , kme , &
1406	! its , ite , jts , jte , kts , kte )
1407	!
1408	! ! surface_input_source=1 => use data from static file (fractional category as input)
1409	! ! surface_input_source=2 => use data from grib file (dominant category as input)
1410	!
1411	! IF ( config_flags%surface_input_source .EQ. 1 ) THEN
1412	! grid%vegcat (its,jts) = 0
1413	! grid%soilcat(its,jts) = 0
1414	! END IF
1415	!
1416	! ! Generate the vegetation and soil category information from the fractional input
1417	! ! data, or use the existing dominant category fields if they exist.
1418	!
1419	! IF ( ( grid%soilcat(its,jts) .LT. 0.5 ) .AND. ( grid%vegcat(its,jts) .LT. 0.5 ) ) THEN
1420	!
1421	! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1422	! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1423	! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1424	!
1425	! CALL process_percent_cat_new ( grid%landmask , &
1426	! grid%landusef , grid%soilctop , grid%soilcbot , &
1427	! grid%isltyp , grid%ivgtyp , &
1428	! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1429	! ids , ide , jds , jde , kds , kde , &
1430	! ims , ime , jms , jme , kms , kme , &
1431	! its , ite , jts , jte , kts , kte , &
1432	! model_config_rec%iswater(grid%id) )
1433	!
1434	! ! Make all the veg/soil parms the same so as not to confuse the developer.
1435	!
1436	! DO j = jts , MIN(jde-1,jte)
1437	! DO i = its , MIN(ide-1,ite)
1438	! grid%vegcat(i,j) = grid%ivgtyp(i,j)
1439	! grid%soilcat(i,j) = grid%isltyp(i,j)
1440	! END DO
1441	! END DO
1442	!
1443	! ELSE
1444	!
1445	! ! Do we have dominant soil and veg data from the input already?
1446	!
1447	! IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
1448	! DO j = jts, MIN(jde-1,jte)
1449	! DO i = its, MIN(ide-1,ite)
1450	! grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
1451	! END DO
1452	! END DO
1453	! END IF
1454	! IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
1455	! DO j = jts, MIN(jde-1,jte)
1456	! DO i = its, MIN(ide-1,ite)
1457	! grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
1458	! END DO
1459	! END DO
1460	! END IF
1461	!
1462	! END IF
1463	!
1464	! ! Land use assignment.
1465	!
1466	! DO j = jts, MIN(jde-1,jte)
1467	! DO i = its, MIN(ide-1,ite)
1468	! grid%lu_index(i,j) = grid%ivgtyp(i,j)
1469	! IF ( grid%lu_index(i,j) .NE. model_config_rec%iswater(grid%id) ) THEN
1470	! grid%landmask(i,j) = 1
1471	! grid%xland(i,j) = 1
1472	! ELSE
1473	! grid%landmask(i,j) = 0
1474	! grid%xland(i,j) = 2
1475	! END IF
1476	! END DO
1477	! END DO
1478	!
1479	! ! Adjust the various soil temperature values depending on the difference in
1480	! ! in elevation between the current model's elevation and the incoming data's
1481	! ! orography.
1482	!
1483	! IF ( flag_soilhgt .EQ. 1 ) THEN
1484	! adjust_soil : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1485	!
1486	! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
1487	! CALL adjust_soil_temp_new ( grid%tmn , model_config_rec%sf_surface_physics(grid%id) , &
1488	! grid%tsk , grid%ht , grid%toposoil , grid%landmask , flag_soilhgt , &
1489	! grid%st000010 , grid%st010040 , grid%st040100 , grid%st100200 , grid%st010200 , &
1490	! flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
1491	! grid%st000007 , grid%st007028 , grid%st028100 , grid%st100255 , &
1492	! flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
1493	! grid%soilt000 , grid%soilt005 , grid%soilt020 , grid%soilt040 , grid%soilt160 , &
1494	! grid%soilt300 , &
1495	! flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , &
1496	! flag_soilt160 , flag_soilt300 , &
1497	! ids , ide , jds , jde , kds , kde , &
1498	! ims , ime , jms , jme , kms , kme , &
1499	! its , ite , jts , jte , kts , kte )
1500	!
1501	! END SELECT adjust_soil
1502	! END IF
1503	!
1504	! ! Fix grid%em_tmn and grid%em_tsk.
1505	!
1506	! fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1507	!
1508	! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
1509	! DO j = jts, MIN(jde-1,jte)
1510	! DO i = its, MIN(ide-1,ite)
1511	! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1512	! ( grid%sst(i,j) .GT. 240. ) .AND. ( grid%sst(i,j) .LT. 350. ) ) THEN
1513	! grid%tmn(i,j) = grid%sst(i,j)
1514	! grid%tsk(i,j) = grid%sst(i,j)
1515	! ELSE IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
1516	! grid%tmn(i,j) = grid%tsk(i,j)
1517	! END IF
1518	! END DO
1519	! END DO
1520	! END SELECT fix_tsk_tmn
1521	!
1522	! ! Is the grid%em_tsk reasonable?
1523	!
1524
1525
1526	!!**** MARS
1527	DO j = jts, MIN(jde-1,jte)
1528	DO i = its, MIN(ide-1,ite)
1529	!!grid%tsk(i,j)=200
1530	grid%tmn(i,j)=0
1531	grid%sst(i,j)=0 !!no use on Mars!!
1532	grid%tslb(i,:,j)=0 !!tslb is 3D field
1533	END DO
1534	END DO
1535	!!**** MARS
1536
1537	! IF ( internal_time_loop .NE. 1 ) THEN
1538	! DO j = jts, MIN(jde-1,jte)
1539	! DO i = its, MIN(ide-1,ite)
1540	! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1541	! grid%tsk(i,j) = grid%em_t_2(i,1,j)
1542	! END IF
1543	! END DO
1544	! END DO
1545	! ELSE
1546	! DO j = jts, MIN(jde-1,jte)
1547	! DO i = its, MIN(ide-1,ite)
1548	! IF ( grid%tsk(i,j) .LT. 170 .or. grid%tsk(i,j) .GT. 400. ) THEN
1549	! print *,'error in the grid%em_tsk'
1550	! print *,'i,j=',i,j
1551	! print *,'grid%landmask=',grid%landmask(i,j)
1552	! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1553	! if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1554	! grid%tsk(i,j)=grid%tmn(i,j)
1555	! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1556	! grid%tsk(i,j)=grid%sst(i,j)
1557	! else
1558	! CALL wrf_error_fatal ( 'grid%em_tsk unreasonable' )
1559	! end if
1560	! END IF
1561	! END DO
1562	! END DO
1563	! END IF
1564	!
1565	! ! Is the grid%em_tmn reasonable?
1566	!
1567	! DO j = jts, MIN(jde-1,jte)
1568	! DO i = its, MIN(ide-1,ite)
1569	! IF ( ( ( grid%tmn(i,j) .LT. 170. ) .OR. ( grid%tmn(i,j) .GT. 400. ) ) &
1570	! .AND. ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1571	! IF ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) THEN
1572	! print *,'error in the grid%em_tmn'
1573	! print *,'i,j=',i,j
1574	! print *,'grid%landmask=',grid%landmask(i,j)
1575	! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1576	! END IF
1577	!
1578	! if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1579	! grid%tmn(i,j)=grid%tsk(i,j)
1580	! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1581	! grid%tmn(i,j)=grid%sst(i,j)
1582	! else
1583	! CALL wrf_error_fatal ( 'grid%em_tmn unreasonable' )
1584	! endif
1585	! END IF
1586	! END DO
1587	! END DO
1588	!
1589	! interpolate_soil_tmw : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1590	!
1591	! CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
1592	! CALL process_soil_real ( grid%tsk , grid%tmn , &
1593	! grid%landmask , grid%sst , &
1594	! st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
1595	! grid%zs , grid%dzs , grid%tslb , grid%smois , grid%sh2o , &
1596	! flag_sst , flag_soilt000, flag_soilm000, &
1597	! ids , ide , jds , jde , kds , kde , &
1598	! ims , ime , jms , jme , kms , kme , &
1599	! its , ite , jts , jte , kts , kte , &
1600	! model_config_rec%sf_surface_physics(grid%id) , &
1601	! model_config_rec%num_soil_layers , &
1602	! model_config_rec%real_data_init_type , &
1603	! num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1604	! num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
1605	!
1606	! END SELECT interpolate_soil_tmw
1607	!
1608	! ! Minimum soil values, residual, from RUC LSM scheme. For input from Noah and using
1609	! ! RUC LSM scheme, this must be subtracted from the input total soil moisture. For
1610	! ! input RUC data and using the Noah LSM scheme, this value must be added to the soil
1611	! ! moisture input.
1612	!
1613	! lqmi(1:num_soil_top_cat) = &
1614	! (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1615	! 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1616	! 0.004, 0.065 /)
1617	!! 0.004, 0.065, 0.020, 0.004, 0.008 /) ! has extra levels for playa, lava, and white sand
1618	!
1619	! ! At the initial time we care about values of soil moisture and temperature, other times are
1620	! ! ignored by the model, so we ignore them, too.
1621	!
1622	! IF ( domain_ClockIsStartTime(grid) ) THEN
1623	! account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1624	!
1625	! CASE ( LSMSCHEME )
1626	! iicount = 0
1627	! IF ( FLAG_SM000010 .EQ. 1 ) THEN
1628	! DO j = jts, MIN(jde-1,jte)
1629	! DO i = its, MIN(ide-1,ite)
1630	! IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
1631	! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1632	! print *,'Noah -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1633	! iicount = iicount + 1
1634	! grid%smois(i,:,j) = 0.005
1635	! END IF
1636	! END DO
1637	! END DO
1638	! IF ( iicount .GT. 0 ) THEN
1639	! print *,'Noah -> Noah: total number of small soil moisture locations = ',iicount
1640	! END IF
1641	! ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1642	! DO j = jts, MIN(jde-1,jte)
1643	! DO i = its, MIN(ide-1,ite)
1644	! grid%smois(i,:,j) = grid%smois(i,:,j) + lqmi(grid%isltyp(i,j))
1645	! END DO
1646	! END DO
1647	! DO j = jts, MIN(jde-1,jte)
1648	! DO i = its, MIN(ide-1,ite)
1649	! IF ( (grid%landmask(i,j).gt.0.5) .and. ( grid%tslb(i,1,j) .gt. 200 ) .and. &
1650	! ( grid%tslb(i,1,j) .lt. 400 ) .and. ( grid%smois(i,1,j) .lt. 0.005 ) ) then
1651	! print *,'RUC -> Noah: bad soil moisture at i,j = ',i,j,grid%smois(i,:,j)
1652	! iicount = iicount + 1
1653	! grid%smois(i,:,j) = 0.005
1654	! END IF
1655	! END DO
1656	! END DO
1657	! IF ( iicount .GT. 0 ) THEN
1658	! print *,'RUC -> Noah: total number of small soil moisture locations = ',iicount
1659	! END IF
1660	! END IF
1661	!
1662	! CASE ( RUCLSMSCHEME )
1663	! iicount = 0
1664	! IF ( FLAG_SM000010 .EQ. 1 ) THEN
1665	! DO j = jts, MIN(jde-1,jte)
1666	! DO i = its, MIN(ide-1,ite)
1667	! grid%smois(i,:,j) = MAX ( grid%smois(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
1668	! END DO
1669	! END DO
1670	! ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1671	! ! no op
1672	! END IF
1673	!
1674	! END SELECT account_for_zero_soil_moisture
1675	! END IF
1676	!
1677	! ! Is the grid%tslb reasonable?
1678	!
1679	! IF ( internal_time_loop .NE. 1 ) THEN
1680	! DO j = jts, MIN(jde-1,jte)
1681	! DO ns = 1 , model_config_rec%num_soil_layers
1682	! DO i = its, MIN(ide-1,ite)
1683	! IF ( grid%tslb(i,ns,j) .LT. 170 .or. grid%tslb(i,ns,j) .GT. 400. ) THEN
1684	! grid%tslb(i,ns,j) = grid%em_t_2(i,1,j)
1685	! grid%smois(i,ns,j) = 0.3
1686	! END IF
1687	! END DO
1688	! END DO
1689	! END DO
1690	! ELSE
1691	! DO j = jts, MIN(jde-1,jte)
1692	! DO i = its, MIN(ide-1,ite)
1693	! IF ( ( ( grid%tslb(i,1,j) .LT. 170. ) .OR. ( grid%tslb(i,1,j) .GT. 400. ) ) .AND. &
1694	! ( grid%landmask(i,j) .GT. 0.5 ) ) THEN
1695	! IF ( ( model_config_rec%sf_surface_physics(grid%id) .NE. LSMSCHEME ) .AND. &
1696	! ( model_config_rec%sf_surface_physics(grid%id) .NE. RUCLSMSCHEME ) ) THEN
1697	! print *,'error in the grid%tslb'
1698	! print *,'i,j=',i,j
1699	! print *,'grid%landmask=',grid%landmask(i,j)
1700	! print *,'grid%tsk, grid%sst, grid%tmn=',grid%tsk(i,j),grid%sst(i,j),grid%tmn(i,j)
1701	! print *,'grid%tslb = ',grid%tslb(i,:,j)
1702	! print *,'old grid%smois = ',grid%smois(i,:,j)
1703	! grid%smois(i,1,j) = 0.3
1704	! grid%smois(i,2,j) = 0.3
1705	! grid%smois(i,3,j) = 0.3
1706	! grid%smois(i,4,j) = 0.3
1707	! END IF
1708	!
1709	! IF ( (grid%tsk(i,j).GT.170. .AND. grid%tsk(i,j).LT.400.) .AND. &
1710	! (grid%tmn(i,j).GT.170. .AND. grid%tmn(i,j).LT.400.) ) THEN
1711	! fake_soil_temp : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1712	! CASE ( SLABSCHEME )
1713	! DO ns = 1 , model_config_rec%num_soil_layers
1714	! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1715	! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1716	! END DO
1717	! CASE ( LSMSCHEME , RUCLSMSCHEME )
1718	! CALL wrf_error_fatal ( 'Assigning constant soil moisture, bad idea')
1719	! DO ns = 1 , model_config_rec%num_soil_layers
1720	! grid%tslb(i,ns,j) = ( grid%tsk(i,j)*(3.0 - grid%zs(ns)) + &
1721	! grid%tmn(i,j)*(0.0 - grid%zs(ns)) ) /(3.0 - 0.0)
1722	! END DO
1723	! END SELECT fake_soil_temp
1724	! else if(grid%tsk(i,j).gt.170. .and. grid%tsk(i,j).lt.400.)then
1725	! CALL wrf_error_fatal ( 'grid%tslb unreasonable 1' )
1726	! DO ns = 1 , model_config_rec%num_soil_layers
1727	! grid%tslb(i,ns,j)=grid%tsk(i,j)
1728	! END DO
1729	! else if(grid%sst(i,j).gt.170. .and. grid%sst(i,j).lt.400.)then
1730	! CALL wrf_error_fatal ( 'grid%tslb unreasonable 2' )
1731	! DO ns = 1 , model_config_rec%num_soil_layers
1732	! grid%tslb(i,ns,j)=grid%sst(i,j)
1733	! END DO
1734	! else if(grid%tmn(i,j).gt.170. .and. grid%tmn(i,j).lt.400.)then
1735	! CALL wrf_error_fatal ( 'grid%tslb unreasonable 3' )
1736	! DO ns = 1 , model_config_rec%num_soil_layers
1737	! grid%tslb(i,ns,j)=grid%tmn(i,j)
1738	! END DO
1739	! else
1740	! CALL wrf_error_fatal ( 'grid%tslb unreasonable 4' )
1741	! endif
1742	! END IF
1743	! END DO
1744	! END DO
1745	! END IF
1746	!
1747	! ! Adjustments for the seaice field AFTER the grid%tslb computations. This is
1748	! ! is for the Noah LSM scheme.
1749	!
1750	! num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1751	! num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1752	! num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1753	! CALL nl_get_seaice_threshold ( grid%id , grid%seaice_threshold )
1754	! CALL nl_get_isice ( grid%id , grid%isice )
1755	! CALL nl_get_iswater ( grid%id , grid%iswater )
1756	! CALL adjust_for_seaice_post ( grid%xice , grid%landmask , grid%tsk , grid%tsk_save , &
1757	! grid%ivgtyp , grid%vegcat , grid%lu_index , &
1758	! grid%xland , grid%landusef , grid%isltyp , grid%soilcat , &
1759	! grid%soilctop , &
1760	! grid%soilcbot , grid%tmn , grid%vegfra , &
1761	! grid%tslb , grid%smois , grid%sh2o , &
1762	! grid%seaice_threshold , &
1763	! num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1764	! model_config_rec%num_soil_layers , &
1765	! grid%iswater , grid%isice , &
1766	! model_config_rec%sf_surface_physics(grid%id) , &
1767	! ids , ide , jds , jde , kds , kde , &
1768	! ims , ime , jms , jme , kms , kme , &
1769	! its , ite , jts , jte , kts , kte )
1770	!
1771	! ! Let us make sure (again) that the grid%landmask and the veg/soil categories match.
1772	!
1773	!oops1=0
1774	!oops2=0
1775	! DO j = jts, MIN(jde-1,jte)
1776	! DO i = its, MIN(ide-1,ite)
1777	! IF ( ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. &
1778	! ( grid%ivgtyp(i,j) .NE. config_flags%iswater .OR. grid%isltyp(i,j) .NE. 14 ) ) .OR. &
1779	! ( ( grid%landmask(i,j) .GT. 0.5 ) .AND. &
1780	! ( grid%ivgtyp(i,j) .EQ. config_flags%iswater .OR. grid%isltyp(i,j) .EQ. 14 ) ) ) THEN
1781	! IF ( grid%tslb(i,1,j) .GT. 1. ) THEN
1782	!oops1=oops1+1
1783	! grid%ivgtyp(i,j) = 5
1784	! grid%isltyp(i,j) = 8
1785	! grid%landmask(i,j) = 1
1786	! grid%xland(i,j) = 1
1787	! ELSE IF ( grid%sst(i,j) .GT. 1. ) THEN
1788	!oops2=oops2+1
1789	! grid%ivgtyp(i,j) = config_flags%iswater
1790	! grid%isltyp(i,j) = 14
1791	! grid%landmask(i,j) = 0
1792	! grid%xland(i,j) = 2
1793	! ELSE
1794	! print *,'the grid%landmask and soil/veg cats do not match'
1795	! print *,'i,j=',i,j
1796	! print *,'grid%landmask=',grid%landmask(i,j)
1797	! print *,'grid%ivgtyp=',grid%ivgtyp(i,j)
1798	! print *,'grid%isltyp=',grid%isltyp(i,j)
1799	! print *,'iswater=', config_flags%iswater
1800	! print *,'grid%tslb=',grid%tslb(i,:,j)
1801	! print *,'grid%sst=',grid%sst(i,j)
1802	! CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
1803	! END IF
1804	! END IF
1805	! END DO
1806	! END DO
1807	!if (oops1.gt.0) then
1808	!print *,'points artificially set to land : ',oops1
1809	!endif
1810	!if(oops2.gt.0) then
1811	!print *,'points artificially set to water: ',oops2
1812	!endif
1813	!! fill grid%sst array with grid%em_tsk if missing in real input (needed for time-varying grid%sst in wrf)
1814	! DO j = jts, MIN(jde-1,jte)
1815	! DO i = its, MIN(ide-1,ite)
1816	! IF ( flag_sst .NE. 1 ) THEN
1817	! grid%sst(i,j) = grid%tsk(i,j)
1818	! ENDIF
1819	! END DO
1820	! END DO
1821
1822
1823	! From the full level data, we can get the half levels, reciprocals, and layer
1824	! thicknesses. These are all defined at half level locations, so one less level.
1825	! We allow the vertical coordinate to accidently come in upside down. We want
1826	! the first full level to be the ground surface.
1827
1828	! Check whether grid%em_znw (full level) data are truly full levels. If not, we need to adjust them
1829	! to be full levels.
1830	! in this test, we check if grid%em_znw(1) is neither 0 nor 1 (within a tolerance of 10**-5)
1831
1832	were_bad = .false.
1833	IF ( ( (grid%em_znw(1).LT.(1-1.E-5) ) .OR. ( grid%em_znw(1).GT.(1+1.E-5) ) ).AND. &
1834	( (grid%em_znw(1).LT.(0-1.E-5) ) .OR. ( grid%em_znw(1).GT.(0+1.E-5) ) ) ) THEN
1835	were_bad = .true.
1836	print *,'Your grid%em_znw input values are probably half-levels. '
1837	print *,grid%em_znw
1838	print *,'WRF expects grid%em_znw values to be full levels. '
1839	print *,'Adjusting now to full levels...'
1840	! We want to ignore the first value if it's negative
1841	IF (grid%em_znw(1).LT.0) THEN
1842	grid%em_znw(1)=0
1843	END IF
1844	DO k=2,kde
1845	grid%em_znw(k)=2*grid%em_znw(k)-grid%em_znw(k-1)
1846	END DO
1847	END IF
1848
1849	! Let's check our changes
1850
1851	IF ( ( ( grid%em_znw(1) .LT. (1-1.E-5) ) .OR. ( grid%em_znw(1) .GT. (1+1.E-5) ) ).AND. &
1852	( ( grid%em_znw(1) .LT. (0-1.E-5) ) .OR. ( grid%em_znw(1) .GT. (0+1.E-5) ) ) ) THEN
1853	print *,'The input grid%em_znw height values were half-levels or erroneous. '
1854	print *,'Attempts to treat the values as half-levels and change them '
1855	print *,'to valid full levels failed.'
1856	CALL wrf_error_fatal("bad grid%em_znw values from input files")
1857	ELSE IF ( were_bad ) THEN
1858	print *,'...adjusted. grid%em_znw array now contains full eta level values. '
1859	ENDIF
1860
1861	IF ( grid%em_znw(1) .LT. grid%em_znw(kde) ) THEN
1862	DO k=1, kde/2
1863	hold_znw = grid%em_znw(k)
1864	grid%em_znw(k)=grid%em_znw(kde+1-k)
1865	grid%em_znw(kde+1-k)=hold_znw
1866	END DO
1867	END IF
1868
1869	DO k=1, kde-1
1870	grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k)
1871	grid%em_rdnw(k) = 1./grid%em_dnw(k)
1872	grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k))
1873	END DO
1874
1875	! Now the same sort of computations with the half eta levels, even ANOTHER
1876	! level less than the one above.
1877
1878	DO k=2, kde-1
1879	grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1))
1880	grid%em_rdn(k) = 1./grid%em_dn(k)
1881	grid%em_fnp(k) = .5* grid%em_dnw(k )/grid%em_dn(k)
1882	grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k)
1883	END DO
1884
1885	! Scads of vertical coefficients.
1886
1887	cof1 = (2.grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))grid%em_dnw(1)/grid%em_dn(2)
1888	cof2 = grid%em_dn(2) /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3)
1889
1890	grid%cf1 = grid%em_fnp(2) + cof1
1891	grid%cf2 = grid%em_fnm(2) - cof1 - cof2
1892	grid%cf3 = cof2
1893
1894	grid%cfn = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1)
1895	grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1)
1896
1897	! Inverse grid distances.
1898
1899	grid%rdx = 1./config_flags%dx
1900	grid%rdy = 1./config_flags%dy
1901
1902	! Some of the many weird geopotential initializations that we'll see today: grid%em_ph0 is total,
1903	! and grid%em_ph_2 is a perturbation from the base state geopotential. We set the base geopotential
1904	! at the lowest level to terrain elevation * gravity.
1905
1906	DO j=jts,jte
1907	DO i=its,ite
1908	grid%em_ph0(i,1,j) = grid%ht(i,j) * g
1909	grid%em_ph_2(i,1,j) = 0.
1910	END DO
1911	END DO
1912
1913	! Base state potential temperature and inverse density (alpha = 1/rho) from
1914	! the half eta levels and the base-profile surface pressure. Compute 1/rho
1915	! from equation of state. The potential temperature is a perturbation from t0.
1916
1917	DO j = jts, MIN(jte,jde-1)
1918	DO i = its, MIN(ite,ide-1)
1919
1920
1921	!****MARS
1922	!TODO: etudier si une meilleure formule n'existe pas pour Mars
1923	!TODO: mais il s'agit juste d'un etat de base ...
1924	!****MARS
1925	! Base state pressure is a function of eta level and terrain, only, plus
1926	! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
1927	! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
1928
1929	!!****MARS
1930	!!ici il s'agit de definir un etat de base, de reference
1931	!!- on ne peut prendre le profil de temperature du modele
1932	!! qui conduit a des instabilites
1933	!! grid%em_t_init(i,k,j)=grid%em_t_2(i,k,j) - t0 est a eviter donc.
1934	!!- pour la pression de surface, aucune information
1935	!! sur un profil de temperature variable et non equilibre
1936	!! ne doit transparaitre
1937	!! p_surf = grid%psfc(i,j) pourquoi pas ... mais t y est utilisee ...
1938	!!
1939	!!>> l'etat de base ne doit dependre "geographiquement" que de la topographie
1940	!!
1941	!!****MARS
1942	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
1943
1944	DO k = 1, kte-1
1945	grid%em_php(i,k,j) = grid%em_znw(k)*(p_surf - grid%p_top) + grid%p_top ! temporary, full lev base pressure
1946	grid%em_pb(i,k,j) = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
1947	! temp = MAX ( 200., t00 + A*LOG(grid%em_pb(i,k,j)/p00) )
1948	temp = t00 + A*LOG(grid%em_pb(i,k,j)/p00)
1949	grid%em_t_init(i,k,j) = temp(p00/grid%em_pb(i,k,j))*(r_d/cp) - t0
1950	grid%em_alb(i,k,j) = (r_d/p1000mb)(grid%em_t_init(i,k,j)+t0)(grid%em_pb(i,k,j)/p1000mb)**cvpm
1951	END DO
1952
1953	! Base state mu is defined as base state surface pressure minus grid%p_top
1954
1955	grid%em_mub(i,j) = p_surf - grid%p_top
1956
1957	! Dry surface pressure is defined as the following (this mu is from the input file
1958	! computed from the dry pressure). Here the dry pressure is just reconstituted.
1959
1960	pd_surf = grid%em_mu0(i,j) + grid%p_top
1961
1962	! Integrate base geopotential, starting at terrain elevation. This assures that
1963	! the base state is in exact hydrostatic balance with respect to the model equations.
1964	! This field is on full levels.
1965
1966	grid%em_phb(i,1,j) = grid%ht(i,j) * g
1967	DO k = 2,kte
1968	grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)grid%em_mub(i,j)grid%em_alb(i,k-1,j)
1969	END DO
1970	END DO
1971	END DO
1972
1973	! Fill in the outer rows and columns to allow us to be sloppy.
1974
1975	IF ( ite .EQ. ide ) THEN
1976	i = ide
1977	DO j = jts, MIN(jde-1,jte)
1978	grid%em_mub(i,j) = grid%em_mub(i-1,j)
1979	grid%em_mu_2(i,j) = grid%em_mu_2(i-1,j)
1980	DO k = 1, kte-1
1981	grid%em_pb(i,k,j) = grid%em_pb(i-1,k,j)
1982	grid%em_t_init(i,k,j) = grid%em_t_init(i-1,k,j)
1983	grid%em_alb(i,k,j) = grid%em_alb(i-1,k,j)
1984	END DO
1985	DO k = 1, kte
1986	grid%em_phb(i,k,j) = grid%em_phb(i-1,k,j)
1987	END DO
1988	END DO
1989	END IF
1990
1991	IF ( jte .EQ. jde ) THEN
1992	j = jde
1993	DO i = its, ite
1994	grid%em_mub(i,j) = grid%em_mub(i,j-1)
1995	grid%em_mu_2(i,j) = grid%em_mu_2(i,j-1)
1996	DO k = 1, kte-1
1997	grid%em_pb(i,k,j) = grid%em_pb(i,k,j-1)
1998	grid%em_t_init(i,k,j) = grid%em_t_init(i,k,j-1)
1999	grid%em_alb(i,k,j) = grid%em_alb(i,k,j-1)
2000	END DO
2001	DO k = 1, kte
2002	grid%em_phb(i,k,j) = grid%em_phb(i,k,j-1)
2003	END DO
2004	END DO
2005	END IF
2006
2007	! Compute the perturbation dry pressure (grid%em_mub + grid%em_mu_2 + ptop = dry grid%em_psfc).
2008
2009	DO j = jts, min(jde-1,jte)
2010	DO i = its, min(ide-1,ite)
2011	grid%em_mu_2(i,j) = grid%em_mu0(i,j) - grid%em_mub(i,j)
2012	END DO
2013	END DO
2014
2015	! Fill in the outer rows and columns to allow us to be sloppy.
2016
2017	IF ( ite .EQ. ide ) THEN
2018	i = ide
2019	DO j = jts, MIN(jde-1,jte)
2020	grid%em_mu_2(i,j) = grid%em_mu_2(i-1,j)
2021	END DO
2022	END IF
2023
2024	IF ( jte .EQ. jde ) THEN
2025	j = jde
2026	DO i = its, ite
2027	grid%em_mu_2(i,j) = grid%em_mu_2(i,j-1)
2028	END DO
2029	END IF
2030
2031	lev500 = 0
2032	DO j = jts, min(jde-1,jte)
2033	DO i = its, min(ide-1,ite)
2034
2035	! Assign the potential temperature (perturbation from t0) and qv on all the mass
2036	! point locations.
2037
2038	DO k = 1 , kde-1
2039	grid%em_t_2(i,k,j) = grid%em_t_2(i,k,j) - t0
2040	END DO
2041
2042	!!---------------------------------------------------------------
2043	!!****MARS: no 500mb adjustment needed
2044	!!****MARS: must keep however the hydrostatic equation integration performed in this loop !
2045	!!****MARS: the DO WHILE loop is deactivated, since we will always be in the case
2046	!!****MARS: ... of "ELSE dpmu = 0."
2047	!!---------------------------------------------------------------
2048	! dpmu = 10001.
2049	! loop_count = 0
2050	!
2051	! DO WHILE ( ( ABS(dpmu) .GT. 10. ) .AND. &
2052	! ( loop_count .LT. 5 ) )
2053	!
2054	! loop_count = loop_count + 1
2055
2056	! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2057	! equation) down from the top to get the pressure perturbation. First get the pressure
2058	! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2059
2060	k = kte-1
2061
2062	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
2063	qvf2 = 1./(1.+qvf1)
2064	qvf1 = qvf1*qvf2
2065
2066	grid%em_p(i,k,j) = - 0.5(grid%em_mu_2(i,j)+qvf1grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2
2067	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2068	grid%em_alt(i,k,j) = (r_d/p1000mb)(grid%em_t_2(i,k,j)+t0)qvf&
2069	(((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)*cvpm)
2070	grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
2071
2072	! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2073	! inverse density fields (total and perturbation).
2074
2075	DO k=kte-2,1,-1
2076	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
2077	qvf2 = 1./(1.+qvf1)
2078	qvf1 = qvf1*qvf2
2079	grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_2(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1)
2080	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2081	grid%em_alt(i,k,j) = (r_d/p1000mb)(grid%em_t_2(i,k,j)+t0)qvf* &
2082	(((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
2083	grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
2084	END DO
2085
2086	! This is the hydrostatic equation used in the model after the small timesteps. In
2087	! the model, grid%em_al (inverse density) is computed from the geopotential.
2088
2089	DO k = 2,kte
2090	grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
2091	grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
2092	+ grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
2093	grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
2094	END DO
2095
2096	! ! Adjust the column pressure so that the computed 500 mb height is close to the
2097	! ! input value (of course, not when we are doing hybrid input).
2098	!
2099	! IF ( ( flag_metgrid .EQ. 1 ) .AND. ( i .EQ. its ) .AND. ( j .EQ. jts ) ) THEN
2100	! DO k = 1 , num_metgrid_levels
2101	! IF ( ABS ( grid%em_p_gc(i,k,j) - 50000. ) .LT. 1. ) THEN
2102	! lev500 = k
2103	! EXIT
2104	! END IF
2105	! END DO
2106	! END IF
2107	!
2108	! ! We only do the adjustment of height if we have the input data on pressure
2109	! ! surfaces, and folks have asked to do this option.
2110	!
2111	! IF ( ( flag_metgrid .EQ. 1 ) .AND. &
2112	! ( config_flags%adjust_heights ) .AND. &
2113	! ( lev500 .NE. 0 ) ) THEN
2114	!
2115	! DO k = 2 , kte-1
2116	!
2117	! ! Get the pressures on the full eta levels (grid%em_php is defined above as
2118	! ! the full-lev base pressure, an easy array to use for 3d space).
2119	!
2120	! pl = grid%em_php(i,k ,j) + &
2121	! ( grid%em_p(i,k-1 ,j) * ( grid%em_znw(k ) - grid%em_znu(k ) ) + &
2122	! grid%em_p(i,k ,j) * ( grid%em_znu(k-1 ) - grid%em_znw(k ) ) ) / &
2123	! ( grid%em_znu(k-1 ) - grid%em_znu(k ) )
2124	! pu = grid%em_php(i,k+1,j) + &
2125	! ( grid%em_p(i,k-1+1,j) * ( grid%em_znw(k +1) - grid%em_znu(k+1) ) + &
2126	! grid%em_p(i,k +1,j) * ( grid%em_znu(k-1+1) - grid%em_znw(k+1) ) ) / &
2127	! ( grid%em_znu(k-1+1) - grid%em_znu(k+1) )
2128	!
2129	! ! If these pressure levels trap 500 mb, use them to interpolate
2130	! ! to the 500 mb level of the computed height.
2131	!!**** PB on MARS .... ?
2132	! IF ( ( pl .GE. 50000. ) .AND. ( pu .LT. 50000. ) ) THEN
2133	! zl = ( grid%em_ph_2(i,k ,j) + grid%em_phb(i,k ,j) ) / g
2134	! zu = ( grid%em_ph_2(i,k+1,j) + grid%em_phb(i,k+1,j) ) / g
2135	!
2136	! z500 = ( zl * ( LOG(50000.) - LOG(pu ) ) + &
2137	! zu * ( LOG(pl ) - LOG(50000.) ) ) / &
2138	! ( LOG(pl) - LOG(pu) )
2139	!! z500 = ( zl * ( (50000.) - (pu ) ) + &
2140	!! zu * ( (pl ) - (50000.) ) ) / &
2141	!! ( (pl) - (pu) )
2142	!
2143	! ! Compute the difference of the 500 mb heights (computed minus input), and
2144	! ! then the change in grid%em_mu_2. The grid%em_php is still full-levels, base pressure.
2145	!
2146	! dz500 = z500 - grid%em_ght_gc(i,lev500,j)
2147	! 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)) &
2148	! (1.+0.6*moist(i,1,j,P_QV))
2149	! dpmu = ( grid%em_php(i,1,j) + grid%em_p(i,1,j) ) * EXP ( g * dz500 / ( r_d * tvsfc ) )
2150	! dpmu = dpmu - ( grid%em_php(i,1,j) + grid%em_p(i,1,j) )
2151	! grid%em_mu_2(i,j) = grid%em_mu_2(i,j) - dpmu
2152	! EXIT
2153	! END IF
2154	!
2155	! END DO
2156	! ELSE
2157	! dpmu = 0.
2158	! END IF
2159	!
2160	! END DO
2161
2162	END DO
2163	END DO
2164
2165	!!****MARS: we use WPS
2166	!
2167	! ! If this is data from the SI, then we probably do not have the original
2168	! ! surface data laying around. Note that these are all the lowest levels
2169	! ! of the respective 3d arrays. For surface pressure, we assume that the
2170	! ! vertical gradient of grid%em_p prime is zilch. This is not all that important.
2171	! ! These are filled in so that the various plotting routines have something
2172	! ! to play with at the initial time for the model.
2173	!
2174	! IF ( flag_metgrid .NE. 1 ) THEN
2175	! DO j = jts, min(jde-1,jte)
2176	! DO i = its, min(ide,ite)
2177	! grid%u10(i,j)=grid%em_u_2(i,1,j)
2178	! END DO
2179	! END DO
2180	!
2181	! DO j = jts, min(jde,jte)
2182	! DO i = its, min(ide-1,ite)
2183	! grid%v10(i,j)=grid%em_v_2(i,1,j)
2184	! END DO
2185	! END DO
2186	!
2187	! DO j = jts, min(jde-1,jte)
2188	! DO i = its, min(ide-1,ite)
2189	! p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
2190	! grid%psfc(i,j)=p_surf + grid%em_p(i,1,j)
2191	! grid%q2(i,j)=moist(i,1,j,P_QV)
2192	! grid%th2(i,j)=grid%em_t_2(i,1,j)+300.
2193	! grid%t2(i,j)=grid%th2(i,j)(((grid%em_p(i,1,j)+grid%em_pb(i,1,j))/p00)*(r_d/cp))
2194	! END DO
2195	! END DO
2196	!
2197	! ! If this data is from WPS, then we have previously assigned the surface
2198	! ! data for u, v, and t. If we have an input qv, welp, we assigned that one,
2199	! ! too. Now we pick up the left overs, and if RH came in - we assign the
2200	! ! mixing ratio.
2201	!
2202	! ELSE IF ( flag_metgrid .EQ. 1 ) THEN
2203	!
2204	!!****MARS: we use WPS
2205
2206	DO j = jts, min(jde-1,jte)
2207	DO i = its, min(ide-1,ite)
2208	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
2209	! recompute the value of surface pressure as calculated by sfcprs2
2210	grid%psfc(i,j)=p_surf + grid%em_p(i,1,j)
2211	grid%th2(i,j)=grid%t2(i,j)(p00/(grid%em_p(i,1,j)+grid%em_pb(i,1,j)))*(r_d/cp)
2212	END DO
2213	END DO
2214	IF ( flag_qv .NE. 1 ) THEN
2215	DO j = jts, min(jde-1,jte)
2216	DO i = its, min(ide-1,ite)
2217	grid%q2(i,j)=moist(i,1,j,P_QV)
2218	END DO
2219	END DO
2220	END IF
2221
2222	! END IF
2223
2224	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2225	#ifdef DM_PARALLEL
2226	# include "HALO_EM_INIT_1.inc"
2227	# include "HALO_EM_INIT_2.inc"
2228	# include "HALO_EM_INIT_3.inc"
2229	# include "HALO_EM_INIT_4.inc"
2230	# include "HALO_EM_INIT_5.inc"
2231	#endif
2232
2233	RETURN
2234
2235	END SUBROUTINE init_domain_rk
2236
2237	!---------------------------------------------------------------------
2238
2239	SUBROUTINE const_module_initialize ( p00 , t00 , a )
2240	USE module_configure
2241	IMPLICIT NONE
2242	! For the real-data-cases only.
2243	REAL , INTENT(OUT) :: p00 , t00 , a
2244	CALL nl_get_base_pres ( 1 , p00 )
2245	CALL nl_get_base_temp ( 1 , t00 )
2246	CALL nl_get_base_lapse ( 1 , a )
2247	END SUBROUTINE const_module_initialize
2248
2249	!-------------------------------------------------------------------
2250
2251	SUBROUTINE rebalance_driver ( grid )
2252
2253	IMPLICIT NONE
2254
2255	TYPE (domain) :: grid
2256
2257	CALL rebalance( grid &
2258	!
2259	#include "em_actual_new_args.inc"
2260	!
2261	)
2262
2263	END SUBROUTINE rebalance_driver
2264
2265	!---------------------------------------------------------------------
2266
2267	SUBROUTINE rebalance ( grid &
2268	!
2269	#include "em_dummy_new_args.inc"
2270	!
2271	)
2272	IMPLICIT NONE
2273
2274	TYPE (domain) :: grid
2275
2276	#include "em_dummy_new_decl.inc"
2277
2278	TYPE (grid_config_rec_type) :: config_flags
2279
2280	REAL :: p_surf , pd_surf, p_surf_int , pb_int , ht_hold
2281	REAL :: qvf , qvf1 , qvf2
2282	REAL :: p00 , t00 , a
2283	REAL , DIMENSION(:,:,:) , ALLOCATABLE :: t_init_int
2284
2285	! Local domain indices and counters.
2286
2287	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
2288
2289	INTEGER :: &
2290	ids, ide, jds, jde, kds, kde, &
2291	ims, ime, jms, jme, kms, kme, &
2292	its, ite, jts, jte, kts, kte, &
2293	ips, ipe, jps, jpe, kps, kpe, &
2294	i, j, k
2295
2296	#ifdef DM_PARALLEL
2297	# include "em_data_calls.inc"
2298	#endif
2299
2300	SELECT CASE ( model_data_order )
2301	CASE ( DATA_ORDER_ZXY )
2302	kds = grid%sd31 ; kde = grid%ed31 ;
2303	ids = grid%sd32 ; ide = grid%ed32 ;
2304	jds = grid%sd33 ; jde = grid%ed33 ;
2305
2306	kms = grid%sm31 ; kme = grid%em31 ;
2307	ims = grid%sm32 ; ime = grid%em32 ;
2308	jms = grid%sm33 ; jme = grid%em33 ;
2309
2310	kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
2311	its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
2312	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
2313
2314	CASE ( DATA_ORDER_XYZ )
2315	ids = grid%sd31 ; ide = grid%ed31 ;
2316	jds = grid%sd32 ; jde = grid%ed32 ;
2317	kds = grid%sd33 ; kde = grid%ed33 ;
2318
2319	ims = grid%sm31 ; ime = grid%em31 ;
2320	jms = grid%sm32 ; jme = grid%em32 ;
2321	kms = grid%sm33 ; kme = grid%em33 ;
2322
2323	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
2324	jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
2325	kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
2326
2327	CASE ( DATA_ORDER_XZY )
2328	ids = grid%sd31 ; ide = grid%ed31 ;
2329	kds = grid%sd32 ; kde = grid%ed32 ;
2330	jds = grid%sd33 ; jde = grid%ed33 ;
2331
2332	ims = grid%sm31 ; ime = grid%em31 ;
2333	kms = grid%sm32 ; kme = grid%em32 ;
2334	jms = grid%sm33 ; jme = grid%em33 ;
2335
2336	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
2337	kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
2338	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
2339
2340	END SELECT
2341
2342	ALLOCATE ( t_init_int(ims:ime,kms:kme,jms:jme) )
2343
2344	! Some of the many weird geopotential initializations that we'll see today: grid%em_ph0 is total,
2345	! and grid%em_ph_2 is a perturbation from the base state geopotential. We set the base geopotential
2346	! at the lowest level to terrain elevation * gravity.
2347
2348	DO j=jts,jte
2349	DO i=its,ite
2350	grid%em_ph0(i,1,j) = grid%ht_fine(i,j) * g
2351	grid%em_ph_2(i,1,j) = 0.
2352	END DO
2353	END DO
2354
2355	! To define the base state, we call a USER MODIFIED routine to set the three
2356	! necessary constants: p00 (sea level pressure, Pa), t00 (sea level temperature, K),
2357	! and A (temperature difference, from 1000 mb to 300 mb, K).
2358
2359	CALL const_module_initialize ( p00 , t00 , a )
2360
2361	! Base state potential temperature and inverse density (alpha = 1/rho) from
2362	! the half eta levels and the base-profile surface pressure. Compute 1/rho
2363	! from equation of state. The potential temperature is a perturbation from t0.
2364
2365	DO j = jts, MIN(jte,jde-1)
2366	DO i = its, MIN(ite,ide-1)
2367
2368	! Base state pressure is a function of eta level and terrain, only, plus
2369	! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
2370	! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
2371	! The fine grid terrain is ht_fine, the interpolated is grid%em_ht.
2372
2373	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht_fine(i,j)/a/r_d ) *0.5 )
2374	p_surf_int = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j) /a/r_d ) *0.5 )
2375
2376	DO k = 1, kte-1
2377	grid%em_pb(i,k,j) = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
2378	pb_int = grid%em_znu(k)*(p_surf_int - grid%p_top) + grid%p_top
2379	grid%em_t_init(i,k,j) = (t00 + ALOG(grid%em_pb(i,k,j)/p00))(p00/grid%em_pb(i,k,j))**(r_d/cp) - t0
2380	t_init_int(i,k,j)= (t00 + ALOG(pb_int /p00))(p00/pb_int )**(r_d/cp) - t0
2381	grid%em_alb(i,k,j) = (r_d/p1000mb)(grid%em_t_init(i,k,j)+t0)(grid%em_pb(i,k,j)/p1000mb)**cvpm
2382	END DO
2383
2384	! Base state mu is defined as base state surface pressure minus grid%p_top
2385
2386	grid%em_mub(i,j) = p_surf - grid%p_top
2387
2388	! Dry surface pressure is defined as the following (this mu is from the input file
2389	! computed from the dry pressure). Here the dry pressure is just reconstituted.
2390
2391	pd_surf = ( grid%em_mub(i,j) + grid%em_mu_2(i,j) ) + grid%p_top
2392
2393	! Integrate base geopotential, starting at terrain elevation. This assures that
2394	! the base state is in exact hydrostatic balance with respect to the model equations.
2395	! This field is on full levels.
2396
2397	grid%em_phb(i,1,j) = grid%ht_fine(i,j) * g
2398	DO k = 2,kte
2399	grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)grid%em_mub(i,j)grid%em_alb(i,k-1,j)
2400	END DO
2401	END DO
2402	END DO
2403
2404	! Replace interpolated terrain with fine grid values.
2405
2406	DO j = jts, MIN(jte,jde-1)
2407	DO i = its, MIN(ite,ide-1)
2408	grid%ht(i,j) = grid%ht_fine(i,j)
2409	END DO
2410	END DO
2411
2412	! Perturbation fields.
2413
2414	DO j = jts, min(jde-1,jte)
2415	DO i = its, min(ide-1,ite)
2416
2417	! The potential temperature is THETAnest = THETAinterp + ( TBARnest - TBARinterp)
2418
2419	DO k = 1 , kde-1
2420	grid%em_t_2(i,k,j) = grid%em_t_2(i,k,j) + ( grid%em_t_init(i,k,j) - t_init_int(i,k,j) )
2421	END DO
2422
2423	! Integrate the hydrostatic equation (from the RHS of the bigstep vertical momentum
2424	! equation) down from the top to get the pressure perturbation. First get the pressure
2425	! perturbation, moisture, and inverse density (total and perturbation) at the top-most level.
2426
2427	k = kte-1
2428
2429	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
2430	qvf2 = 1./(1.+qvf1)
2431	qvf1 = qvf1*qvf2
2432
2433	grid%em_p(i,k,j) = - 0.5(grid%em_mu_2(i,j)+qvf1grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2
2434	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2435	grid%em_alt(i,k,j) = (r_d/p1000mb)(grid%em_t_2(i,k,j)+t0)qvf* &
2436	(((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
2437	grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
2438
2439	! Now, integrate down the column to compute the pressure perturbation, and diagnose the two
2440	! inverse density fields (total and perturbation).
2441
2442	DO k=kte-2,1,-1
2443	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
2444	qvf2 = 1./(1.+qvf1)
2445	qvf1 = qvf1*qvf2
2446	grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_2(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1)
2447	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
2448	grid%em_alt(i,k,j) = (r_d/p1000mb)(grid%em_t_2(i,k,j)+t0)qvf* &
2449	(((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
2450	grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
2451	END DO
2452
2453	! This is the hydrostatic equation used in the model after the small timesteps. In
2454	! the model, grid%em_al (inverse density) is computed from the geopotential.
2455
2456	DO k = 2,kte
2457	grid%em_ph_2(i,k,j) = grid%em_ph_2(i,k-1,j) - &
2458	grid%em_dnw(k-1) * ( (grid%em_mub(i,j)+grid%em_mu_2(i,j))*grid%em_al(i,k-1,j) &
2459	+ grid%em_mu_2(i,j)*grid%em_alb(i,k-1,j) )
2460	grid%em_ph0(i,k,j) = grid%em_ph_2(i,k,j) + grid%em_phb(i,k,j)
2461	END DO
2462
2463	END DO
2464	END DO
2465
2466	DEALLOCATE ( t_init_int )
2467
2468	ips = its ; ipe = ite ; jps = jts ; jpe = jte ; kps = kts ; kpe = kte
2469	#ifdef DM_PARALLEL
2470	# include "HALO_EM_INIT_1.inc"
2471	# include "HALO_EM_INIT_2.inc"
2472	# include "HALO_EM_INIT_3.inc"
2473	# include "HALO_EM_INIT_4.inc"
2474	# include "HALO_EM_INIT_5.inc"
2475	#endif
2476	END SUBROUTINE rebalance
2477
2478	!---------------------------------------------------------------------
2479
2480	RECURSIVE SUBROUTINE find_my_parent ( grid_ptr_in , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2481
2482	USE module_domain
2483
2484	TYPE(domain) , POINTER :: grid_ptr_in , grid_ptr_out
2485	TYPE(domain) , POINTER :: grid_ptr_sibling
2486	INTEGER :: id_wanted , id_i_am
2487	LOGICAL :: found_the_id
2488
2489	found_the_id = .FALSE.
2490	grid_ptr_sibling => grid_ptr_in
2491	DO WHILE ( ASSOCIATED ( grid_ptr_sibling ) )
2492
2493	IF ( grid_ptr_sibling%grid_id .EQ. id_wanted ) THEN
2494	found_the_id = .TRUE.
2495	grid_ptr_out => grid_ptr_sibling
2496	RETURN
2497	ELSE IF ( grid_ptr_sibling%num_nests .GT. 0 ) THEN
2498	grid_ptr_sibling => grid_ptr_sibling%nests(1)%ptr
2499	CALL find_my_parent ( grid_ptr_sibling , grid_ptr_out , id_i_am , id_wanted , found_the_id )
2500	ELSE
2501	grid_ptr_sibling => grid_ptr_sibling%sibling
2502	END IF
2503
2504	END DO
2505
2506	END SUBROUTINE find_my_parent
2507
2508	#endif
2509
2510	!---------------------------------------------------------------------
2511
2512	#ifdef VERT_UNIT
2513
2514	!This is a main program for a small unit test for the vertical interpolation.
2515
2516	program vint
2517
2518	implicit none
2519
2520	integer , parameter :: ij = 3
2521	integer , parameter :: keta = 30
2522	integer , parameter :: kgen =20
2523
2524	integer :: ids , ide , jds , jde , kds , kde , &
2525	ims , ime , jms , jme , kms , kme , &
2526	its , ite , jts , jte , kts , kte
2527
2528	integer :: generic
2529
2530	real , dimension(1:ij,kgen,1:ij) :: fo , po
2531	real , dimension(1:ij,1:keta,1:ij) :: fn_calc , fn_interp , pn
2532
2533	integer, parameter :: interp_type = 1 ! 2
2534	! integer, parameter :: lagrange_order = 2 ! 1
2535	integer :: lagrange_order
2536	logical, parameter :: lowest_lev_from_sfc = .FALSE. ! .TRUE.
2537	real , parameter :: zap_close_levels = 500. ! 100.
2538	integer, parameter :: force_sfc_in_vinterp = 0 ! 6
2539
2540	integer :: k
2541
2542	ids = 1 ; ide = ij ; jds = 1 ; jde = ij ; kds = 1 ; kde = keta
2543	ims = 1 ; ime = ij ; jms = 1 ; jme = ij ; kms = 1 ; kme = keta
2544	its = 1 ; ite = ij ; jts = 1 ; jte = ij ; kts = 1 ; kte = keta
2545
2546	generic = kgen
2547
2548	print *,' '
2549	print *,'------------------------------------'
2550	print *,'UNIT TEST FOR VERTICAL INTERPOLATION'
2551	print *,'------------------------------------'
2552	print *,' '
2553	do lagrange_order = 1 , 2
2554	print *,' '
2555	print *,'------------------------------------'
2556	print *,'Lagrange Order = ',lagrange_order
2557	print *,'------------------------------------'
2558	print *,' '
2559	call fillitup ( fo , po , fn_calc , pn , &
2560	ids , ide , jds , jde , kds , kde , &
2561	ims , ime , jms , jme , kms , kme , &
2562	its , ite , jts , jte , kts , kte , &
2563	generic , lagrange_order )
2564
2565	print *,' '
2566	print *,'Level Pressure Field'
2567	print *,' (Pa) (generic)'
2568	print *,'------------------------------------'
2569	print *,' '
2570	do k = 1 , generic
2571	write (*,fmt='(i2,2x,f12.3,1x,g15.8)' ) &
2572	k,po(2,k,2),fo(2,k,2)
2573	end do
2574	print *,' '
2575
2576	call vert_interp ( fo , po , fn_interp , pn , &
2577	generic , 'T' , &
2578	interp_type , lagrange_order , lowest_lev_from_sfc , &
2579	zap_close_levels , force_sfc_in_vinterp , &
2580	ids , ide , jds , jde , kds , kde , &
2581	ims , ime , jms , jme , kms , kme , &
2582	its , ite , jts , jte , kts , kte )
2583
2584	print *,'Multi-Order Interpolator'
2585	print *,'------------------------------------'
2586	print *,' '
2587	print *,'Level Pressure Field Field Field'
2588	print *,' (Pa) Calc Interp Diff'
2589	print *,'------------------------------------'
2590	print *,' '
2591	do k = kts , kte-1
2592	write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2593	k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
2594	end do
2595
2596	call vert_interp_old ( fo , po , fn_interp , pn , &
2597	generic , 'T' , &
2598	interp_type , lagrange_order , lowest_lev_from_sfc , &
2599	zap_close_levels , force_sfc_in_vinterp , &
2600	ids , ide , jds , jde , kds , kde , &
2601	ims , ime , jms , jme , kms , kme , &
2602	its , ite , jts , jte , kts , kte )
2603
2604	print *,'Linear Interpolator'
2605	print *,'------------------------------------'
2606	print *,' '
2607	print *,'Level Pressure Field Field Field'
2608	print *,' (Pa) Calc Interp Diff'
2609	print *,'------------------------------------'
2610	print *,' '
2611	do k = kts , kte-1
2612	write (*,fmt='(i2,2x,f12.3,1x,3(g15.7))' ) &
2613	k,pn(2,k,2),fn_calc(2,k,2),fn_interp(2,k,2),fn_calc(2,k,2)-fn_interp(2,k,2)
2614	end do
2615	end do
2616
2617	end program vint
2618
2619	subroutine wrf_error_fatal (string)
2620	character (len=*) :: string
2621	print *,string
2622	stop
2623	end subroutine wrf_error_fatal
2624
2625	subroutine fillitup ( fo , po , fn , pn , &
2626	ids , ide , jds , jde , kds , kde , &
2627	ims , ime , jms , jme , kms , kme , &
2628	its , ite , jts , jte , kts , kte , &
2629	generic , lagrange_order )
2630
2631	implicit none
2632
2633	integer , intent(in) :: ids , ide , jds , jde , kds , kde , &
2634	ims , ime , jms , jme , kms , kme , &
2635	its , ite , jts , jte , kts , kte
2636
2637	integer , intent(in) :: generic , lagrange_order
2638
2639	real , dimension(ims:ime,generic,jms:jme) , intent(out) :: fo , po
2640	real , dimension(ims:ime,kms:kme,jms:jme) , intent(out) :: fn , pn
2641
2642	integer :: i , j , k
2643
2644	real , parameter :: piov2 = 3.14159265358 / 2.
2645
2646	k = 1
2647	do j = jts , jte
2648	do i = its , ite
2649	po(i,k,j) = 102000.
2650	end do
2651	end do
2652
2653	do k = 2 , generic
2654	do j = jts , jte
2655	do i = its , ite
2656	po(i,k,j) = ( 5000. * ( 1 - (k-1) ) + 100000. * ( (k-1) - (generic-1) ) ) / (1. - real(generic-1) )
2657	end do
2658	end do
2659	end do
2660
2661	if ( lagrange_order .eq. 1 ) then
2662	do k = 1 , generic
2663	do j = jts , jte
2664	do i = its , ite
2665	fo(i,k,j) = po(i,k,j)
2666	! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2667	end do
2668	end do
2669	end do
2670	else if ( lagrange_order .eq. 2 ) then
2671	do k = 1 , generic
2672	do j = jts , jte
2673	do i = its , ite
2674	fo(i,k,j) = (((po(i,k,j)-5000.)/102000.)((102000.-po(i,k,j))/102000.))102000.
2675	! fo(i,k,j) = sin(po(i,k,j) * piov2 / 102000. )
2676	end do
2677	end do
2678	end do
2679	end if
2680
2681	!!!!!!!!!!!!
2682
2683	do k = kts , kte
2684	do j = jts , jte
2685	do i = its , ite
2686	pn(i,k,j) = ( 5000. * ( 0 - (k-1) ) + 102000. * ( (k-1) - (kte-1) ) ) / (-1. * real(kte-1) )
2687	end do
2688	end do
2689	end do
2690
2691	do k = kts , kte-1
2692	do j = jts , jte
2693	do i = its , ite
2694	pn(i,k,j) = ( pn(i,k,j) + pn(i,k+1,j) ) /2.
2695	end do
2696	end do
2697	end do
2698
2699
2700	if ( lagrange_order .eq. 1 ) then
2701	do k = kts , kte-1
2702	do j = jts , jte
2703	do i = its , ite
2704	fn(i,k,j) = pn(i,k,j)
2705	! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2706	end do
2707	end do
2708	end do
2709	else if ( lagrange_order .eq. 2 ) then
2710	do k = kts , kte-1
2711	do j = jts , jte
2712	do i = its , ite
2713	fn(i,k,j) = (((pn(i,k,j)-5000.)/102000.)((102000.-pn(i,k,j))/102000.))102000.
2714	! fn(i,k,j) = sin(pn(i,k,j) * piov2 / 102000. )
2715	end do
2716	end do
2717	end do
2718	end if
2719
2720	end subroutine fillitup
2721
2722	#endif
2723
2724	!---------------------------------------------------------------------
2725
2726	SUBROUTINE vert_interp ( fo , po , fnew , pnu , &
2727	generic , var_type , &
2728	interp_type , lagrange_order , lowest_lev_from_sfc , &
2729	zap_close_levels , force_sfc_in_vinterp , &
2730	ids , ide , jds , jde , kds , kde , &
2731	ims , ime , jms , jme , kms , kme , &
2732	its , ite , jts , jte , kts , kte )
2733
2734	! Vertically interpolate the new field. The original field on the original
2735	! pressure levels is provided, and the new pressure surfaces to interpolate to.
2736
2737	IMPLICIT NONE
2738
2739	INTEGER , INTENT(IN) :: interp_type , lagrange_order
2740	LOGICAL , INTENT(IN) :: lowest_lev_from_sfc
2741	REAL , INTENT(IN) :: zap_close_levels
2742	INTEGER , INTENT(IN) :: force_sfc_in_vinterp
2743	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2744	ims , ime , jms , jme , kms , kme , &
2745	its , ite , jts , jte , kts , kte
2746	INTEGER , INTENT(IN) :: generic
2747
2748	CHARACTER (LEN=1) :: var_type
2749
2750	REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: fo , po
2751	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
2752	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
2753
2754	REAL , DIMENSION(ims:ime,generic,jms:jme) :: forig , porig
2755	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
2756
2757	! Local vars
2758
2759	INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2 , knext
2760	INTEGER :: istart , iend , jstart , jend , kstart , kend
2761	INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
2762	INTEGER , DIMENSION(ims:ime ) :: ks
2763	INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
2764	INTEGER :: count , zap , kst
2765
2766	LOGICAL :: any_below_ground
2767
2768	REAL :: p1 , p2 , pn, hold
2769	REAL , DIMENSION(1:generic) :: ordered_porig , ordered_forig
2770	REAL , DIMENSION(kts:kte) :: ordered_pnew , ordered_fnew
2771
2772	!****MARS
2773	!big problems ... discontinuity in the interpolated fields ...
2774	print *, '25/05/2007: decided to use simple linear interpolations'
2775	print *, 'use that one at your own risk'
2776	!stop
2777	!****MARS
2778
2779
2780	! Horiontal loop bounds for different variable types.
2781
2782	IF ( var_type .EQ. 'U' ) THEN
2783	istart = its
2784	iend = ite
2785	jstart = jts
2786	jend = MIN(jde-1,jte)
2787	kstart = kts
2788	kend = kte-1
2789	DO j = jstart,jend
2790	DO k = 1,generic
2791	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2792	porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
2793	END DO
2794	END DO
2795	IF ( ids .EQ. its ) THEN
2796	DO k = 1,generic
2797	porig(its,k,j) = po(its,k,j)
2798	END DO
2799	END IF
2800	IF ( ide .EQ. ite ) THEN
2801	DO k = 1,generic
2802	porig(ite,k,j) = po(ite-1,k,j)
2803	END DO
2804	END IF
2805
2806	DO k = kstart,kend
2807	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
2808	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
2809	END DO
2810	END DO
2811	IF ( ids .EQ. its ) THEN
2812	DO k = kstart,kend
2813	pnew(its,k,j) = pnu(its,k,j)
2814	END DO
2815	END IF
2816	IF ( ide .EQ. ite ) THEN
2817	DO k = kstart,kend
2818	pnew(ite,k,j) = pnu(ite-1,k,j)
2819	END DO
2820	END IF
2821	END DO
2822	ELSE IF ( var_type .EQ. 'V' ) THEN
2823	istart = its
2824	iend = MIN(ide-1,ite)
2825	jstart = jts
2826	jend = jte
2827	kstart = kts
2828	kend = kte-1
2829	DO i = istart,iend
2830	DO k = 1,generic
2831	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2832	porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
2833	END DO
2834	END DO
2835	IF ( jds .EQ. jts ) THEN
2836	DO k = 1,generic
2837	porig(i,k,jts) = po(i,k,jts)
2838	END DO
2839	END IF
2840	IF ( jde .EQ. jte ) THEN
2841	DO k = 1,generic
2842	porig(i,k,jte) = po(i,k,jte-1)
2843	END DO
2844	END IF
2845
2846	DO k = kstart,kend
2847	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
2848	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
2849	END DO
2850	END DO
2851	IF ( jds .EQ. jts ) THEN
2852	DO k = kstart,kend
2853	pnew(i,k,jts) = pnu(i,k,jts)
2854	END DO
2855	END IF
2856	IF ( jde .EQ. jte ) THEN
2857	DO k = kstart,kend
2858	pnew(i,k,jte) = pnu(i,k,jte-1)
2859	END DO
2860	END IF
2861	END DO
2862	ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
2863	istart = its
2864	iend = MIN(ide-1,ite)
2865	jstart = jts
2866	jend = MIN(jde-1,jte)
2867	kstart = kts
2868	kend = kte
2869	DO j = jstart,jend
2870	DO k = 1,generic
2871	DO i = istart,iend
2872	porig(i,k,j) = po(i,k,j)
2873	END DO
2874	END DO
2875
2876	DO k = kstart,kend
2877	DO i = istart,iend
2878	pnew(i,k,j) = pnu(i,k,j)
2879	END DO
2880	END DO
2881	END DO
2882	ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
2883	istart = its
2884	iend = MIN(ide-1,ite)
2885	jstart = jts
2886	jend = MIN(jde-1,jte)
2887	kstart = kts
2888	kend = kte-1
2889	DO j = jstart,jend
2890	DO k = 1,generic
2891	DO i = istart,iend
2892	porig(i,k,j) = po(i,k,j)
2893	END DO
2894	END DO
2895
2896	DO k = kstart,kend
2897	DO i = istart,iend
2898	pnew(i,k,j) = pnu(i,k,j)
2899	END DO
2900	END DO
2901	END DO
2902	ELSE
2903	istart = its
2904	iend = MIN(ide-1,ite)
2905	jstart = jts
2906	jend = MIN(jde-1,jte)
2907	kstart = kts
2908	kend = kte-1
2909	DO j = jstart,jend
2910	DO k = 1,generic
2911	DO i = istart,iend
2912	porig(i,k,j) = po(i,k,j)
2913	END DO
2914	END DO
2915
2916	DO k = kstart,kend
2917	DO i = istart,iend
2918	pnew(i,k,j) = pnu(i,k,j)
2919	END DO
2920	END DO
2921	END DO
2922	END IF
2923
2924	DO j = jstart , jend
2925
2926	! The lowest level is the surface. Levels 2 through "generic" are supposed to
2927	! be "bottom-up". Flip if they are not. This is based on the input pressure
2928	! array.
2929
2930	IF ( porig(its,2,j) .LT. porig(its,generic,j) ) THEN
2931	DO kn = 2 , ( generic + 1 ) / 2
2932	DO i = istart , iend
2933	hold = porig(i,kn,j)
2934	porig(i,kn,j) = porig(i,generic+2-kn,j)
2935	porig(i,generic+2-kn,j) = hold
2936	forig(i,kn,j) = fo (i,generic+2-kn,j)
2937	forig(i,generic+2-kn,j) = fo (i,kn,j)
2938	END DO
2939	DO i = istart , iend
2940	forig(i,1,j) = fo (i,1,j)
2941	END DO
2942	END DO
2943	ELSE
2944	DO kn = 1 , generic
2945	DO i = istart , iend
2946	forig(i,kn,j) = fo (i,kn,j)
2947	END DO
2948	END DO
2949	END IF
2950
2951	! Skip all of the levels below ground in the original data based upon the surface pressure.
2952	! The ko_above_sfc is the index in the pressure array that is above the surface. If there
2953	! are no levels underground, this is index = 2. The remaining levels are eligible for use
2954	! in the vertical interpolation.
2955
2956	DO i = istart , iend
2957	ko_above_sfc(i) = -1
2958	END DO
2959	DO ko = kstart+1 , kend
2960	DO i = istart , iend
2961	IF ( ko_above_sfc(i) .EQ. -1 ) THEN
2962	IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
2963	ko_above_sfc(i) = ko
2964	END IF
2965	END IF
2966	END DO
2967	END DO
2968
2969	! Piece together columns of the original input data. Pass the vertical columns to
2970	! the iterpolator.
2971
2972	DO i = istart , iend
2973
2974	! If the surface value is in the middle of the array, three steps: 1) do the
2975	! values below the ground (this is just to catch the occasional value that is
2976	! inconsistently below the surface based on input data), 2) do the surface level, then
2977	! 3) add in the levels that are above the surface. For the levels next to the surface,
2978	! we check to remove any levels that are "too close". When building the column of input
2979	! pressures, we also attend to the request for forcing the surface analysis to be used
2980	! in a few lower eta-levels.
2981
2982	! How many levels have we skipped in the input column.
2983
2984	zap = 0
2985
2986	! Fill in the column from up to the level just below the surface with the input
2987	! presssure and the input field (orig or old, which ever). For an isobaric input
2988	! file, this data is isobaric.
2989
2990	IF ( ko_above_sfc(i) .GT. 2 ) THEN
2991	count = 1
2992	DO ko = 2 , ko_above_sfc(i)-1
2993	ordered_porig(count) = porig(i,ko,j)
2994	ordered_forig(count) = forig(i,ko,j)
2995	count = count + 1
2996	END DO
2997
2998	! Make sure the pressure just below the surface is not "too close", this
2999	! will cause havoc with the higher order interpolators. In case of a "too close"
3000	! instance, we toss out the offending level (NOT the surface one) by simply
3001	! decrementing the accumulating loop counter.
3002
3003	IF ( ordered_porig(count-1) - porig(i,1,j) .LT. zap_close_levels ) THEN
3004	count = count -1
3005	zap = 1
3006	END IF
3007
3008	! Add in the surface values.
3009
3010	ordered_porig(count) = porig(i,1,j)
3011	ordered_forig(count) = forig(i,1,j)
3012	count = count + 1
3013
3014	! A usual way to do the vertical interpolation is to pay more attention to the
3015	! surface data. Why? Well it has about 20x the density as the upper air, so we
3016	! hope the analysis is better there. We more strongly use this data by artificially
3017	! tossing out levels above the surface that are beneath a certain number of prescribed
3018	! eta levels at this (i,j). The "zap" value is how many levels of input we are
3019	! removing, which is used to tell the interpolator how many valid values are in
3020	! the column. The "count" value is the increment to the index of levels, and is
3021	! only used for assignments.
3022
3023	IF ( force_sfc_in_vinterp .GT. 0 ) THEN
3024
3025	! Get the pressure at the eta level. We want to remove all input pressure levels
3026	! between the level above the surface to the pressure at this eta surface. That
3027	! forces the surface value to be used through the selected eta level. Keep track
3028	! of two things: the level to use above the eta levels, and how many levels we are
3029	! skipping.
3030
3031	knext = ko_above_sfc(i)
3032	find_level : DO ko = ko_above_sfc(i) , generic
3033	IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
3034	knext = ko
3035	exit find_level
3036	ELSE
3037	zap = zap + 1
3038	END IF
3039	END DO find_level
3040
3041	! No request for special interpolation, so we just assign the next level to use
3042	! above the surface as, ta da, the first level above the surface. I know, wow.
3043
3044	ELSE
3045	knext = ko_above_sfc(i)
3046	END IF
3047
3048	! One more time, make sure the pressure just above the surface is not "too close", this
3049	! will cause havoc with the higher order interpolators. In case of a "too close"
3050	! instance, we toss out the offending level above the surface (NOT the surface one) by simply
3051	! incrementing the loop counter. Here, count-1 is the surface level and knext is either
3052	! the next level up OR it is the level above the prescribed number of eta surfaces.
3053
3054	IF ( ordered_porig(count-1) - porig(i,knext,j) .LT. zap_close_levels ) THEN
3055	kst = knext+1
3056	zap = zap + 1
3057	ELSE
3058	kst = knext
3059	END IF
3060
3061	DO ko = kst , generic
3062	ordered_porig(count) = porig(i,ko,j)
3063	ordered_forig(count) = forig(i,ko,j)
3064	count = count + 1
3065	END DO
3066
3067	! This is easy, the surface is the lowest level, just stick them in, in this order. OK,
3068	! there are a couple of subtleties. We have to check for that special interpolation that
3069	! skips some input levels so that the surface is used for the lowest few eta levels. Also,
3070	! we must macke sure that we still do not have levels that are "too close" together.
3071
3072	ELSE
3073
3074	! Initialize no input levels have yet been removed from consideration.
3075
3076	zap = 0
3077
3078	! The surface is the lowest level, so it gets set right away to location 1.
3079
3080	ordered_porig(1) = porig(i,1,j)
3081	ordered_forig(1) = forig(i,1,j)
3082
3083	! We start filling in the array at loc 2, as in just above the level we just stored.
3084
3085	count = 2
3086
3087	! Are we forcing the interpolator to skip valid input levels so that the
3088	! surface data is used through more levels? Essentially as above.
3089
3090	IF ( force_sfc_in_vinterp .GT. 0 ) THEN
3091	knext = 2
3092	find_level2: DO ko = 2 , generic
3093	IF ( porig(i,ko,j) .LE. pnew(i,force_sfc_in_vinterp,j) ) THEN
3094	knext = ko
3095	exit find_level2
3096	ELSE
3097	zap = zap + 1
3098	END IF
3099	END DO find_level2
3100	ELSE
3101	knext = 2
3102	END IF
3103
3104	! Fill in the data above the surface. The "knext" index is either the one
3105	! just above the surface OR it is the index associated with the level that
3106	! is just above the pressure at this (i,j) of the top eta level that is to
3107	! be directly impacted with the surface level in interpolation.
3108
3109	DO ko = knext , generic
3110	IF ( ordered_porig(count-1) - porig(i,ko,j) .LT. zap_close_levels ) THEN
3111	zap = zap + 1
3112	CYCLE
3113	END IF
3114	ordered_porig(count) = porig(i,ko,j)
3115	ordered_forig(count) = forig(i,ko,j)
3116	count = count + 1
3117	END DO
3118
3119	END IF
3120
3121	! Now get the column of the "new" pressure data. So, this one is easy.
3122
3123	DO kn = kstart , kend
3124	ordered_pnew(kn) = pnew(i,kn,j)
3125	END DO
3126
3127	! The polynomials are either in pressure or LOG(pressure).
3128
3129	IF ( interp_type .EQ. 1 ) THEN
3130	CALL lagrange_setup ( var_type , &
3131	ordered_porig , ordered_forig , generic-zap , lagrange_order , &
3132	ordered_pnew , ordered_fnew , kend-kstart+1 ,i,j)
3133	ELSE
3134	CALL lagrange_setup ( var_type , &
3135	LOG(ordered_porig(1:generic-zap)) , ordered_forig , generic-zap , lagrange_order , &
3136	LOG(ordered_pnew(kstart:kend)) , ordered_fnew , kend-kstart+1 ,i,j)
3137	END IF
3138
3139	! Save the computed data.
3140
3141	DO kn = kstart , kend
3142	fnew(i,kn,j) = ordered_fnew(kn)
3143	END DO
3144
3145	! There may have been a request to have the surface data from the input field
3146	! to be assigned as to the lowest eta level. This assumes thin layers (usually
3147	! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
3148
3149	IF ( lowest_lev_from_sfc ) THEN
3150	fnew(i,1,j) = forig(i,ko_above_sfc(i)-1,j)
3151	END IF
3152
3153	END DO
3154
3155	END DO
3156
3157	END SUBROUTINE vert_interp
3158
3159	!---------------------------------------------------------------------
3160
3161	SUBROUTINE vert_interp_old ( forig , po , fnew , pnu , &
3162	generic , var_type , &
3163	interp_type , lagrange_order , lowest_lev_from_sfc , &
3164	zap_close_levels , force_sfc_in_vinterp , &
3165	ids , ide , jds , jde , kds , kde , &
3166	ims , ime , jms , jme , kms , kme , &
3167	its , ite , jts , jte , kts , kte )
3168
3169	! Vertically interpolate the new field. The original field on the original
3170	! pressure levels is provided, and the new pressure surfaces to interpolate to.
3171
3172	IMPLICIT NONE
3173
3174	INTEGER , INTENT(IN) :: interp_type , lagrange_order
3175	LOGICAL , INTENT(IN) :: lowest_lev_from_sfc
3176	REAL , INTENT(IN) :: zap_close_levels
3177	INTEGER , INTENT(IN) :: force_sfc_in_vinterp
3178	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3179	ims , ime , jms , jme , kms , kme , &
3180	its , ite , jts , jte , kts , kte
3181	INTEGER , INTENT(IN) :: generic
3182
3183	CHARACTER (LEN=1) :: var_type
3184
3185	! REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: forig , po
3186	!****MARS
3187	!error with g95 and warning with pgf90
3188	REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(IN) :: po
3189	REAL , DIMENSION(ims:ime,generic,jms:jme) , INTENT(INOUT) :: forig
3190	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: pnu
3191	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: fnew
3192
3193	REAL , DIMENSION(ims:ime,generic,jms:jme) :: porig
3194	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: pnew
3195
3196	! Local vars
3197
3198	INTEGER :: i , j , k , ko , kn , k1 , k2 , ko_1 , ko_2
3199	INTEGER :: istart , iend , jstart , jend , kstart , kend
3200	INTEGER , DIMENSION(ims:ime,kms:kme ) :: k_above , k_below
3201	INTEGER , DIMENSION(ims:ime ) :: ks
3202	INTEGER , DIMENSION(ims:ime ) :: ko_above_sfc
3203
3204	LOGICAL :: any_below_ground
3205
3206	REAL :: p1 , p2 , pn
3207	!****MARS
3208	integer vert_extrap
3209	integer kn_save
3210	vert_extrap = 0
3211	kn_save = 0
3212	!****MARS
3213
3214	! Horizontal loop bounds for different variable types.
3215
3216	IF ( var_type .EQ. 'U' ) THEN
3217	istart = its
3218	iend = ite
3219	jstart = jts
3220	jend = MIN(jde-1,jte)
3221	kstart = kts
3222	kend = kte-1
3223	DO j = jstart,jend
3224	DO k = 1,generic
3225	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
3226	porig(i,k,j) = ( po(i,k,j) + po(i-1,k,j) ) * 0.5
3227	END DO
3228	END DO
3229	IF ( ids .EQ. its ) THEN
3230	DO k = 1,generic
3231	porig(its,k,j) = po(its,k,j)
3232	END DO
3233	END IF
3234	IF ( ide .EQ. ite ) THEN
3235	DO k = 1,generic
3236	porig(ite,k,j) = po(ite-1,k,j)
3237	END DO
3238	END IF
3239
3240	DO k = kstart,kend
3241	DO i = MAX(ids+1,its) , MIN(ide-1,ite)
3242	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i-1,k,j) ) * 0.5
3243	END DO
3244	END DO
3245	IF ( ids .EQ. its ) THEN
3246	DO k = kstart,kend
3247	pnew(its,k,j) = pnu(its,k,j)
3248	END DO
3249	END IF
3250	IF ( ide .EQ. ite ) THEN
3251	DO k = kstart,kend
3252	pnew(ite,k,j) = pnu(ite-1,k,j)
3253	END DO
3254	END IF
3255	END DO
3256	ELSE IF ( var_type .EQ. 'V' ) THEN
3257	istart = its
3258	iend = MIN(ide-1,ite)
3259	jstart = jts
3260	jend = jte
3261	kstart = kts
3262	kend = kte-1
3263	DO i = istart,iend
3264	DO k = 1,generic
3265	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3266	porig(i,k,j) = ( po(i,k,j) + po(i,k,j-1) ) * 0.5
3267	END DO
3268	END DO
3269	IF ( jds .EQ. jts ) THEN
3270	DO k = 1,generic
3271	porig(i,k,jts) = po(i,k,jts)
3272	END DO
3273	END IF
3274	IF ( jde .EQ. jte ) THEN
3275	DO k = 1,generic
3276	porig(i,k,jte) = po(i,k,jte-1)
3277	END DO
3278	END IF
3279
3280	DO k = kstart,kend
3281	DO j = MAX(jds+1,jts) , MIN(jde-1,jte)
3282	pnew(i,k,j) = ( pnu(i,k,j) + pnu(i,k,j-1) ) * 0.5
3283	END DO
3284	END DO
3285	IF ( jds .EQ. jts ) THEN
3286	DO k = kstart,kend
3287	pnew(i,k,jts) = pnu(i,k,jts)
3288	END DO
3289	END IF
3290	IF ( jde .EQ. jte ) THEN
3291	DO k = kstart,kend
3292	pnew(i,k,jte) = pnu(i,k,jte-1)
3293	END DO
3294	END IF
3295	END DO
3296	ELSE IF ( ( var_type .EQ. 'W' ) .OR. ( var_type .EQ. 'Z' ) ) THEN
3297	istart = its
3298	iend = MIN(ide-1,ite)
3299	jstart = jts
3300	jend = MIN(jde-1,jte)
3301	kstart = kts
3302	kend = kte
3303	DO j = jstart,jend
3304	DO k = 1,generic
3305	DO i = istart,iend
3306	porig(i,k,j) = po(i,k,j)
3307	END DO
3308	END DO
3309
3310	DO k = kstart,kend
3311	DO i = istart,iend
3312	pnew(i,k,j) = pnu(i,k,j)
3313	END DO
3314	END DO
3315	END DO
3316	ELSE IF ( ( var_type .EQ. 'T' ) .OR. ( var_type .EQ. 'Q' ) ) THEN
3317	istart = its
3318	iend = MIN(ide-1,ite)
3319	jstart = jts
3320	jend = MIN(jde-1,jte)
3321	kstart = kts
3322	kend = kte-1
3323	DO j = jstart,jend
3324	DO k = 1,generic
3325	DO i = istart,iend
3326	porig(i,k,j) = po(i,k,j)
3327	END DO
3328	END DO
3329
3330	DO k = kstart,kend
3331	DO i = istart,iend
3332	pnew(i,k,j) = pnu(i,k,j)
3333	END DO
3334	END DO
3335	END DO
3336	ELSE
3337	istart = its
3338	iend = MIN(ide-1,ite)
3339	jstart = jts
3340	jend = MIN(jde-1,jte)
3341	kstart = kts
3342	kend = kte-1
3343	DO j = jstart,jend
3344	DO k = 1,generic
3345	DO i = istart,iend
3346	porig(i,k,j) = po(i,k,j)
3347	END DO
3348	END DO
3349
3350	DO k = kstart,kend
3351	DO i = istart,iend
3352	pnew(i,k,j) = pnu(i,k,j)
3353	END DO
3354	END DO
3355	END DO
3356	END IF
3357
3358
3359	DO j = jstart , jend
3360
3361	! Skip all of the levels below ground in the original data based upon the surface pressure.
3362	! The ko_above_sfc is the index in the pressure array that is above the surface. If there
3363	! are no levels underground, this is index = 2. The remaining levels are eligible for use
3364	! in the vertical interpolation.
3365
3366	DO i = istart , iend
3367	ko_above_sfc(i) = -1
3368	END DO
3369	DO ko = kstart+1 , kend
3370	DO i = istart , iend
3371
3372	IF ( ko_above_sfc(i) .EQ. -1 ) THEN
3373	IF ( porig(i,1,j) .GT. porig(i,ko,j) ) THEN
3374	ko_above_sfc(i) = ko
3375	!!****MARS
3376	!!old stuff
3377	!!
3378	!! Pressure level may be OK, however data from the diagfi is possibly missing
3379	!IF (forig(i,ko,j) .EQ. -1.0e+30) THEN
3380	! ko_above_sfc(i) = -1
3381	!END IF
3382	! !! Once the right start level is found, check that it is OK
3383	! !! >> first column should be 1e30 or so, second column should be a realistic value
3384	! !IF ( ko_above_sfc(i) .NE. -1 ) THEN
3385	! ! print *, 'verif', forig(i,ko-1,j), forig(i,ko,j), forig(i,ko+1,j), ko
3386	! !END IF
3387	!!
3388	!!****MARS
3389	END IF
3390	END IF
3391
3392	END DO
3393	END DO
3394
3395	! Initialize interpolation location. These are the levels in the original pressure
3396	! data that are physically below and above the targeted new pressure level.
3397
3398	DO kn = kts , kte
3399	DO i = its , ite
3400	k_above(i,kn) = -1
3401	k_below(i,kn) = -2
3402	END DO
3403	END DO
3404
3405	! Starting location is no lower than previous found location. This is for O(n logn)
3406	! and not O(n^2), where n is the number of vertical levels to search.
3407
3408	DO i = its , ite
3409	ks(i) = 1
3410	END DO
3411
3412	! Find trapping layer for interpolation. The kn index runs through all of the "new"
3413	! levels of data.
3414
3415	DO kn = kstart , kend
3416
3417	DO i = istart , iend
3418
3419	! For each "new" level (kn), we search to find the trapping levels in the "orig"
3420	! data. Most of the time, the "new" levels are the eta surfaces, and the "orig"
3421	! levels are the input pressure levels.
3422
3423	found_trap_above : DO ko = ks(i) , generic-1
3424
3425	! Because we can have levels in the interpolation that are not valid,
3426	! let's toss out any candidate orig pressure values that are below ground
3427	! based on the surface pressure. If the level =1, then this IS the surface
3428	! level, so we HAVE to keep that one, but maybe not the ones above. If the
3429	! level (ks) is NOT=1, then we have to just CYCLE our loop to find a legit
3430	! below-pressure value. If we are not below ground, then we choose two
3431	! neighboring levels to test whether they surround the new pressure level.
3432
3433	! The input trapping levels that we are trying is the surface and the first valid
3434	! level above the surface.
3435
3436	IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .EQ. 1 ) ) THEN
3437	ko_1 = ko
3438	ko_2 = ko_above_sfc(i)
3439	!!****MARS
3440	!!old remark: the possible issue is fixed later in the code ...
3441	!!****MARS
3442
3443	! The "below" level is underground, cycle until we get to a valid pressure
3444	! above ground.
3445
3446	ELSE IF ( ( ko .LT. ko_above_sfc(i) ) .AND. ( ko .NE. 1 ) ) THEN
3447	CYCLE found_trap_above
3448
3449	! The "below" level is above the surface, so we are in the clear to test these
3450	! two levels out.
3451
3452	ELSE
3453	ko_1 = ko
3454	ko_2 = ko+1
3455
3456	END IF
3457
3458
3459	! The test of the candidate levels: "below" has to have a larger pressure, and
3460	! "above" has to have a smaller pressure.
3461
3462	! OK, we found the correct two surrounding levels. The locations are saved for use in the
3463	! interpolation.
3464
3465	IF ( ( porig(i,ko_1,j) .GE. pnew(i,kn,j) ) .AND. &
3466	( porig(i,ko_2,j) .LT. pnew(i,kn,j) ) ) THEN
3467	k_above(i,kn) = ko_2
3468	k_below(i,kn) = ko_1
3469	ks(i) = ko_1
3470	EXIT found_trap_above
3471
3472	! What do we do is we need to extrapolate the data underground? This happens when the
3473	! lowest pressure that we have is physically "above" the new target pressure. Our
3474	! actions depend on the type of variable we are interpolating.
3475
3476	ELSE IF ( porig(i,1,j) .LT. pnew(i,kn,j) ) THEN
3477	!!****MARS
3478	!!old stuff
3479	!!check: values are usually quite close
3480	!print *,porig(i,1,j),pnew(i,kn,j)
3481	!!****MARS
3482
3483	! For horizontal winds and moisture, we keep a constant value under ground.
3484
3485	IF ( ( var_type .EQ. 'U' ) .OR. &
3486	( var_type .EQ. 'V' ) .OR. &
3487	( var_type .EQ. 'Q' ) ) THEN
3488	k_above(i,kn) = 1
3489	ks(i) = 1
3490
3491	! For temperature and height, we extrapolate the data. Hopefully, we are not
3492	! extrapolating too far. For pressure level input, the eta levels are always
3493	! contained within the surface to p_top levels, so no extrapolation is ever
3494	! required.
3495
3496	ELSE IF ( ( var_type .EQ. 'Z' ) .OR. &
3497	( var_type .EQ. 'T' ) ) THEN
3498	k_above(i,kn) = ko_above_sfc(i)
3499	k_below(i,kn) = 1
3500	ks(i) = 1
3501	!!!****MARS
3502	!!old stuff
3503	!k_above(i,kn) = 1
3504	!ks(i) = 1
3505	!!!"Hopefully, we are not extrapolating too far"
3506	!!!>> true on Mars ??
3507	!!!****MARS
3508
3509	! Just a catch all right now.
3510
3511	ELSE
3512	k_above(i,kn) = 1
3513	ks(i) = 1
3514	END IF
3515
3516	EXIT found_trap_above
3517
3518	! The other extrapolation that might be required is when we are going above the
3519	! top level of the input data. Usually this means we chose a P_PTOP value that
3520	! was inappropriate, and we should stop and let someone fix this mess.
3521
3522	ELSE IF ( porig(i,generic,j) .GT. pnew(i,kn,j) ) THEN
3523	print *,'data is too high, try a lower p_top'
3524	print *,'pnew=',pnew(i,kn,j),'i',i,'j',j,'kn',kn
3525	print *,'pnew=',pnew(i,:,j)
3526	print *,'porig=',porig(i,:,j)
3527	CALL wrf_error_fatal ('requested p_top is higher than input data, lower p_top')
3528
3529	END IF
3530	END DO found_trap_above
3531	END DO
3532	END DO
3533
3534	! Linear vertical interpolation.
3535
3536	DO kn = kstart , kend
3537	DO i = istart , iend
3538	IF ( k_above(i,kn) .EQ. 1 ) THEN
3539	!!!****MARS
3540	!!old stuff
3541	!!!ne doit pas arriver avec la temperature si l'on definit bien le champ au sol
3542	!IF (forig(i,1,j) .EQ. -1.0e+30) THEN
3543	! print *,'no data here - surface - var is ...',var_type,i,j,1
3544	! print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j)
3545	! forig(i,1,j) = forig(i,ko_above_sfc(i),j)
3546	! !IF ( ( var_type .EQ. 'U' ) .OR. &
3547	! ! ( var_type .EQ. 'V' ) .OR. &
3548	! ! ( var_type .EQ. 'Q' ) ) THEN
3549	! ! print *,'zero wind at the ground'
3550	! ! forig(i,1,j) = 0
3551	! !ENDIF
3552	! IF (forig(i,1,j) .EQ. -1.0e+30) THEN
3553	! print *,'well ... are you sure ?'
3554	! stop
3555	! ENDIF
3556	!END IF
3557	!!!****MARS
3558	fnew(i,kn,j) = forig(i,1,j)
3559	ELSE
3560	k2 = MAX ( k_above(i,kn) , 2)
3561	k1 = MAX ( k_below(i,kn) , 1)
3562	IF ( k1 .EQ. k2 ) THEN
3563	CALL wrf_error_fatal ( 'identical values in the interp, bad for divisions' )
3564	END IF
3565	!!!****MARS
3566	!!old stuff
3567	!IF (forig(i,k2,j) .EQ. -1.0e+30) THEN
3568	! print *,'no data here - level above - you_d better stop',i,j,k2
3569	! stop
3570	!END IF
3571	!IF (forig(i,k1,j) .EQ. -1.0e+30) THEN
3572	! print *,'no data here - level below - var is ...',var_type,i,j,k1
3573	! print *,'setting to first level with data...',ko_above_sfc(i),porig(i,ko_above_sfc(i),j)
3574	! forig(i,k1,j) = forig(i,ko_above_sfc(i),j)
3575	! !!!VERIFIER QUE LA TEMPERATURE AU SOL N'EST PAS CONCERNEE
3576	! !!!(montagnes=sources locales de chaleur)
3577	! !!!normalement, pas de souci, et lors de l'exécution rien ne s'affiche
3578	!END IF
3579	!!!****MARS
3580	IF ( interp_type .EQ. 1 ) THEN
3581	p1 = porig(i,k1,j)
3582	p2 = porig(i,k2,j)
3583	pn = pnew(i,kn,j)
3584	ELSE IF ( interp_type .EQ. 2 ) THEN
3585	p1 = ALOG(porig(i,k1,j))
3586	p2 = ALOG(porig(i,k2,j))
3587	pn = ALOG(pnew(i,kn,j))
3588	END IF
3589	IF ( ( p1-pn) * (p2-pn) > 0. ) THEN
3590	! CALL wrf_error_fatal ( 'both trapping pressures are on the same side of the new pressure' )
3591	! CALL wrf_debug ( 0 , 'both trapping pressures are on the same side of the new pressure' )
3592	!!!****MARS
3593	vert_extrap = vert_extrap + 1
3594	!print *, 'extrapolate', pnew(i,kn,j)-porig(i,k1,j), 'for WRF level', kn
3595	IF (kn_save < kn) kn_save=kn
3596	!!!****MARS
3597	END IF
3598	fnew(i,kn,j) = ( forig(i,k1,j) * ( p2 - pn ) + &
3599	forig(i,k2,j) * ( pn - p1 ) ) / &
3600	( p2 - p1 )
3601	END IF
3602	END DO
3603	END DO
3604
3605	search_below_ground : DO kn = kstart , kend
3606	any_below_ground = .FALSE.
3607	DO i = istart , iend
3608	IF ( k_above(i,kn) .EQ. 1 ) THEN
3609	fnew(i,kn,j) = forig(i,1,j)
3610	any_below_ground = .TRUE.
3611	END IF
3612	END DO
3613	IF ( .NOT. any_below_ground ) THEN
3614	EXIT search_below_ground
3615	END IF
3616	END DO search_below_ground
3617
3618	! There may have been a request to have the surface data from the input field
3619	! to be assigned as to the lowest eta level. This assumes thin layers (usually
3620	! the isobaric original field has the surface from 2-m T and RH, and 10-m U and V).
3621
3622
3623	DO i = istart , iend
3624	IF ( lowest_lev_from_sfc ) THEN
3625	fnew(i,1,j) = forig(i,ko_above_sfc(i),j)
3626	END IF
3627	END DO
3628
3629	END DO
3630	print *,'VERT EXTRAP = ', vert_extrap
3631	print *,'finished with ... ', var_type
3632	print *,'max WRF eta level where extrap. occurs: ',kn_save
3633
3634	END SUBROUTINE vert_interp_old
3635
3636	!---------------------------------------------------------------------
3637
3638	SUBROUTINE lagrange_setup ( var_type , all_x , all_y , all_dim , n , target_x , target_y , target_dim ,i,j)
3639
3640	! We call a Lagrange polynomial interpolator. The parallel concerns are put off as this
3641	! is initially set up for vertical use. The purpose is an input column of pressure (all_x),
3642	! and the associated pressure level data (all_y). These are assumed to be sorted (ascending
3643	! or descending, no matter). The locations to be interpolated to are the pressures in
3644	! target_x, probably the new vertical coordinate values. The field that is output is the
3645	! target_y, which is defined at the target_x location. Mostly we expect to be 2nd order
3646	! overlapping polynomials, with only a single 2nd order method near the top and bottom.
3647	! When n=1, this is linear; when n=2, this is a second order interpolator.
3648
3649	IMPLICIT NONE
3650
3651	CHARACTER (LEN=1) :: var_type
3652	INTEGER , INTENT(IN) :: all_dim , n , target_dim
3653	REAL, DIMENSION(all_dim) , INTENT(IN) :: all_x , all_y
3654	REAL , DIMENSION(target_dim) , INTENT(IN) :: target_x
3655	REAL , DIMENSION(target_dim) , INTENT(OUT) :: target_y
3656
3657	! Brought in for debug purposes, all of the computations are in a single column.
3658
3659	INTEGER , INTENT(IN) :: i,j
3660
3661	! Local vars
3662
3663	REAL , DIMENSION(n+1) :: x , y
3664	REAL :: target_y_1 , target_y_2
3665	LOGICAL :: found_loc
3666	INTEGER :: loop , loc_center_left , loc_center_right , ist , iend , target_loop
3667
3668
3669	IF ( all_dim .LT. n+1 ) THEN
3670	print *,'all_dim = ',all_dim
3671	print *,'order = ',n
3672	print *,'i,j = ',i,j
3673	print *,'p array = ',all_x
3674	print *,'f array = ',all_y
3675	print *,'p target= ',target_x
3676	CALL wrf_error_fatal ( 'troubles, the interpolating order is too large for this few input values' )
3677	END IF
3678
3679	IF ( n .LT. 1 ) THEN
3680	CALL wrf_error_fatal ( 'pal, linear is about as low as we go' )
3681	END IF
3682
3683	! Loop over the list of target x and y values.
3684
3685	DO target_loop = 1 , target_dim
3686
3687	! Find the two trapping x values, and keep the indices.
3688
3689	found_loc = .FALSE.
3690	find_trap : DO loop = 1 , all_dim -1
3691	IF ( ( target_x(target_loop) - all_x(loop) ) * ( target_x(target_loop) - all_x(loop+1) ) .LE. 0.0 ) THEN
3692	loc_center_left = loop
3693	loc_center_right = loop+1
3694	found_loc = .TRUE.
3695	!****MARS: check if no errors here
3696	!print *,'interpolating ... ',var_type
3697	! print *,'i,j = ',i,j
3698	! print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
3699	! DO loop = 1 , all_dim
3700	! print *,'column of pressure and value = ',all_x(loop),all_y(loop)
3701	! END DO
3702	!END IF
3703	!****MARS
3704	EXIT find_trap
3705	END IF
3706	END DO find_trap
3707
3708	IF ( ( .NOT. found_loc ) .AND. ( target_x(target_loop) .GT. all_x(1) ) ) THEN
3709	IF ( var_type .EQ. 'T' ) THEN
3710	write(6,fmt='(A,2i5,2f11.3)') &
3711	' --> extrapolating TEMPERATURE near sfc: i,j,psfc, p target = ',&
3712	i,j,all_x(1),target_x(target_loop)
3713	target_y(target_loop) = ( all_y(1) * ( target_x(target_loop) - all_x(2) ) + &
3714	all_y(2) * ( all_x(1) - target_x(target_loop) ) ) / &
3715	( all_x(1) - all_x(2) )
3716	ELSE
3717	!write(6,fmt='(A,2i5,2f11.3)') &
3718	!' --> extrapolating zero gradient near sfc: i,j,psfc, p target = ',&
3719	!i,j,all_x(1),target_x(target_loop)
3720	target_y(target_loop) = all_y(1)
3721	END IF
3722	CYCLE
3723	ELSE IF ( .NOT. found_loc ) THEN
3724	!****MARS: normally, no errors here (otherwise, keep this part commented ?)
3725	print *, var_type
3726	print *,'i,j = ',i,j
3727	print *,'target pressure and value = ',target_x(target_loop),target_y(target_loop)
3728	DO loop = 1 , all_dim
3729	print *,'column of pressure and value = ',all_x(loop),all_y(loop)
3730	END DO
3731	CALL wrf_error_fatal ( 'troubles, could not find trapping x locations' )
3732	!****MARS: end of 'keep this part commented'
3733	END IF
3734
3735	! Even or odd order? We can put the value in the middle if this is
3736	! an odd order interpolator. For the even guys, we'll do it twice
3737	! and shift the range one index, then get an average.
3738
3739	IF ( MOD(n,2) .NE. 0 ) THEN
3740	IF ( ( loc_center_left -(((n+1)/2)-1) .GE. 1 ) .AND. &
3741	( loc_center_right+(((n+1)/2)-1) .LE. all_dim ) ) THEN
3742	ist = loc_center_left -(((n+1)/2)-1)
3743	iend = iend + n
3744	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3745	ELSE
3746	IF ( .NOT. found_loc ) THEN
3747	CALL wrf_error_fatal ( 'I doubt this will happen, I will only do 2nd order for now' )
3748	END IF
3749	END IF
3750
3751	ELSE IF ( MOD(n,2) .EQ. 0 ) THEN
3752	IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
3753	( loc_center_right+(((n )/2) ) .LE. all_dim ) .AND. &
3754	( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
3755	( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
3756	ist = loc_center_left -(((n )/2)-1)
3757	iend = ist + n
3758	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_1 )
3759	ist = loc_center_left -(((n )/2) )
3760	iend = ist + n
3761	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y_2 )
3762	target_y(target_loop) = ( target_y_1 + target_y_2 ) * 0.5
3763
3764	ELSE IF ( ( loc_center_left -(((n )/2)-1) .GE. 1 ) .AND. &
3765	( loc_center_right+(((n )/2) ) .LE. all_dim ) ) THEN
3766	ist = loc_center_left -(((n )/2)-1)
3767	iend = ist + n
3768	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3769	ELSE IF ( ( loc_center_left -(((n )/2) ) .GE. 1 ) .AND. &
3770	( loc_center_right+(((n )/2)-1) .LE. all_dim ) ) THEN
3771	ist = loc_center_left -(((n )/2) )
3772	iend = ist + n
3773	CALL lagrange_interp ( all_x(ist:iend) , all_y(ist:iend) , n , target_x(target_loop) , target_y(target_loop) )
3774	ELSE
3775	CALL wrf_error_fatal ( 'unauthorized area, you should not be here' )
3776	END IF
3777
3778	END IF
3779
3780	END DO
3781
3782	END SUBROUTINE lagrange_setup
3783
3784	!---------------------------------------------------------------------
3785
3786	SUBROUTINE lagrange_interp ( x , y , n , target_x , target_y )
3787
3788	! Interpolation using Lagrange polynomials.
3789	! P(x) = f(x0)Ln0(x) + ... + f(xn)Lnn(x)
3790	! where Lnk(x) = (x -x0)(x -x1)...(x -xk-1)(x -xk+1)...(x -xn)
3791	! ---------------------------------------------
3792	! (xk-x0)(xk-x1)...(xk-xk-1)(xk-xk+1)...(xk-xn)
3793
3794	IMPLICIT NONE
3795
3796	INTEGER , INTENT(IN) :: n
3797	REAL , DIMENSION(0:n) , INTENT(IN) :: x , y
3798	REAL , INTENT(IN) :: target_x
3799
3800	REAL , INTENT(OUT) :: target_y
3801
3802	! Local vars
3803
3804	INTEGER :: i , k
3805	REAL :: numer , denom , Px
3806	REAL , DIMENSION(0:n) :: Ln
3807
3808	Px = 0.
3809	DO i = 0 , n
3810	numer = 1.
3811	denom = 1.
3812	DO k = 0 , n
3813	IF ( k .EQ. i ) CYCLE
3814	numer = numer * ( target_x - x(k) )
3815	denom = denom * ( x(i) - x(k) )
3816	END DO
3817	Ln(i) = y(i) * numer / denom
3818	Px = Px + Ln(i)
3819	END DO
3820	target_y = Px
3821
3822	END SUBROUTINE lagrange_interp
3823
3824	#ifndef VERT_UNIT
3825	!---------------------------------------------------------------------
3826
3827	SUBROUTINE p_dry ( mu0 , eta , pdht , pdry , &
3828	ids , ide , jds , jde , kds , kde , &
3829	ims , ime , jms , jme , kms , kme , &
3830	its , ite , jts , jte , kts , kte )
3831
3832	! Compute reference pressure and the reference mu.
3833
3834	IMPLICIT NONE
3835
3836	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3837	ims , ime , jms , jme , kms , kme , &
3838	its , ite , jts , jte , kts , kte
3839
3840	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: mu0
3841	REAL , DIMENSION( kms:kme ) , INTENT(IN) :: eta
3842	REAL :: pdht
3843	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pdry
3844
3845	! Local vars
3846
3847	INTEGER :: i , j , k
3848	REAL , DIMENSION( kms:kme ) :: eta_h
3849
3850	DO k = kts , kte-1
3851	eta_h(k) = ( eta(k) + eta(k+1) ) * 0.5
3852	END DO
3853
3854	DO j = jts , MIN ( jde-1 , jte )
3855	DO k = kts , kte-1
3856	DO i = its , MIN (ide-1 , ite )
3857	pdry(i,k,j) = eta_h(k) * mu0(i,j) + pdht
3858	END DO
3859	END DO
3860	END DO
3861
3862	END SUBROUTINE p_dry
3863
3864	!---------------------------------------------------------------------
3865
3866	SUBROUTINE p_dts ( pdts , intq , psfc , p_top , &
3867	ids , ide , jds , jde , kds , kde , &
3868	ims , ime , jms , jme , kms , kme , &
3869	its , ite , jts , jte , kts , kte )
3870
3871	! Compute difference between the dry, total surface pressure and the top pressure.
3872
3873	IMPLICIT NONE
3874
3875	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3876	ims , ime , jms , jme , kms , kme , &
3877	its , ite , jts , jte , kts , kte
3878
3879	REAL , INTENT(IN) :: p_top
3880	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: psfc
3881	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: intq
3882	REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: pdts
3883
3884	! Local vars
3885
3886	INTEGER :: i , j , k
3887
3888	DO j = jts , MIN ( jde-1 , jte )
3889	DO i = its , MIN (ide-1 , ite )
3890	pdts(i,j) = psfc(i,j) - intq(i,j) - p_top
3891	END DO
3892	END DO
3893
3894	END SUBROUTINE p_dts
3895
3896	!---------------------------------------------------------------------
3897
3898	SUBROUTINE p_dhs ( pdhs , ht , p0 , t0 , a , &
3899	ids , ide , jds , jde , kds , kde , &
3900	ims , ime , jms , jme , kms , kme , &
3901	its , ite , jts , jte , kts , kte )
3902
3903	! Compute dry, hydrostatic surface pressure.
3904
3905	IMPLICIT NONE
3906
3907	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3908	ims , ime , jms , jme , kms , kme , &
3909	its , ite , jts , jte , kts , kte
3910
3911	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(IN) :: ht
3912	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: pdhs
3913
3914	REAL , INTENT(IN) :: p0 , t0 , a
3915
3916	! Local vars
3917
3918	INTEGER :: i , j , k
3919	!****MARS ....
3920	REAL , PARAMETER :: Rd = 192.
3921	REAL , PARAMETER :: g = 3.72
3922	print *,'compute dry, hydrostatic surface pressure'
3923	!****MARS ....
3924
3925	DO j = jts , MIN ( jde-1 , jte )
3926	DO i = its , MIN (ide-1 , ite )
3927	pdhs(i,j) = p0 * EXP ( -t0/a + SQRT ( (t0/a)*2 - 2. g * ht(i,j)/(a * Rd) ) )
3928	END DO
3929	END DO
3930
3931	!****MARS
3932	!****MARS cette formule est-elle juste sur Mars ?
3933	!****MARS >> a premiere vue, ne donne pas de resultats absurdes
3934	!****TODO: il y a peut etre meilleur !
3935	!****MARS
3936
3937	!print *,pdhs
3938	!stop
3939
3940
3941	END SUBROUTINE p_dhs
3942
3943	!---------------------------------------------------------------------
3944
3945	SUBROUTINE find_p_top ( p , p_top , &
3946	ids , ide , jds , jde , kds , kde , &
3947	ims , ime , jms , jme , kms , kme , &
3948	its , ite , jts , jte , kts , kte )
3949
3950	! Find the largest pressure in the top level. This is our p_top. We are
3951	! assuming that the top level is the location where the pressure is a minimum
3952	! for each column. In cases where the top surface is not isobaric, a
3953	! communicated value must be shared in the calling routine. Also in cases
3954	! where the top surface is not isobaric, care must be taken that the new
3955	! maximum pressure is not greater than the previous value. This test is
3956	! also handled in the calling routine.
3957
3958	IMPLICIT NONE
3959
3960	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3961	ims , ime , jms , jme , kms , kme , &
3962	its , ite , jts , jte , kts , kte
3963
3964	REAL :: p_top
3965	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
3966
3967	! Local vars
3968
3969	INTEGER :: i , j , k, min_lev
3970
3971	i = its
3972	j = jts
3973	p_top = p(i,2,j)
3974	min_lev = 2
3975	DO k = 2 , kte
3976	IF ( p_top .GT. p(i,k,j) ) THEN
3977	p_top = p(i,k,j)
3978	min_lev = k
3979	END IF
3980	END DO
3981
3982	k = min_lev
3983	p_top = p(its,k,jts)
3984	DO j = jts , MIN ( jde-1 , jte )
3985	DO i = its , MIN (ide-1 , ite )
3986	p_top = MAX ( p_top , p(i,k,j) )
3987	END DO
3988	END DO
3989
3990	END SUBROUTINE find_p_top
3991
3992	!---------------------------------------------------------------------
3993
3994	SUBROUTINE t_to_theta ( t , p , p00 , &
3995	ids , ide , jds , jde , kds , kde , &
3996	ims , ime , jms , jme , kms , kme , &
3997	its , ite , jts , jte , kts , kte )
3998
3999	! Compute dry, hydrostatic surface pressure.
4000
4001	IMPLICIT NONE
4002
4003	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4004	ims , ime , jms , jme , kms , kme , &
4005	its , ite , jts , jte , kts , kte
4006
4007	REAL , INTENT(IN) :: p00
4008	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p
4009	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: t
4010
4011	! Local vars
4012
4013	INTEGER :: i , j , k
4014	!****MARS
4015	REAL , PARAMETER :: Rd = 192.
4016	REAL , PARAMETER :: Cp = 844.6
4017	!****MARS
4018
4019	DO j = jts , MIN ( jde-1 , jte )
4020	DO k = kts , kte
4021	DO i = its , MIN (ide-1 , ite )
4022	t(i,k,j) = t(i,k,j) * ( p00 / p(i,k,j) ) ** (Rd / Cp)
4023	END DO
4024	END DO
4025	END DO
4026
4027	END SUBROUTINE t_to_theta
4028
4029	!---------------------------------------------------------------------
4030
4031	SUBROUTINE integ_moist ( q_in , p_in , pd_out , t_in , ght_in , intq , &
4032	ids , ide , jds , jde , kds , kde , &
4033	ims , ime , jms , jme , kms , kme , &
4034	its , ite , jts , jte , kts , kte )
4035
4036	! Integrate the moisture field vertically. Mostly used to get the total
4037	! vapor pressure, which can be subtracted from the total pressure to get
4038	! the dry pressure.
4039
4040	IMPLICIT NONE
4041
4042	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4043	ims , ime , jms , jme , kms , kme , &
4044	its , ite , jts , jte , kts , kte
4045
4046	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: q_in , p_in , t_in , ght_in
4047	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: pd_out
4048	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: intq
4049
4050	! Local vars
4051
4052	INTEGER :: i , j , k
4053	INTEGER , DIMENSION(ims:ime) :: level_above_sfc
4054	REAL , DIMENSION(ims:ime,jms:jme) :: psfc , tsfc , qsfc, zsfc
4055	REAL , DIMENSION(ims:ime,kms:kme) :: q , p , t , ght, pd
4056
4057	REAL :: rhobar , qbar , dz
4058	REAL :: p1 , p2 , t1 , t2 , q1 , q2 , z1, z2
4059
4060	LOGICAL :: upside_down
4061
4062	!****MARS
4063	REAL , PARAMETER :: Rd = 192.
4064	REAL , PARAMETER :: g = 3.72
4065	!****MARS
4066
4067
4068	! Get a surface value, always the first level of a 3d field.
4069
4070	DO j = jts , MIN ( jde-1 , jte )
4071	DO i = its , MIN (ide-1 , ite )
4072	psfc(i,j) = p_in(i,kts,j)
4073	tsfc(i,j) = t_in(i,kts,j)
4074	qsfc(i,j) = q_in(i,kts,j)
4075	zsfc(i,j) = ght_in(i,kts,j)
4076	END DO
4077	END DO
4078
4079	IF ( p_in(its,kts+1,jts) .LT. p_in(its,kte,jts) ) THEN
4080	upside_down = .TRUE.
4081	ELSE
4082	upside_down = .FALSE.
4083	END IF
4084
4085	DO j = jts , MIN ( jde-1 , jte )
4086
4087	! Initialize the integrated quantity of moisture to zero.
4088
4089	DO i = its , MIN (ide-1 , ite )
4090	intq(i,j) = 0.
4091	END DO
4092
4093	IF ( upside_down ) THEN
4094	DO i = its , MIN (ide-1 , ite )
4095	p(i,kts) = p_in(i,kts,j)
4096	t(i,kts) = t_in(i,kts,j)
4097	q(i,kts) = q_in(i,kts,j)
4098	ght(i,kts) = ght_in(i,kts,j)
4099	DO k = kts+1,kte
4100	p(i,k) = p_in(i,kte+2-k,j)
4101	t(i,k) = t_in(i,kte+2-k,j)
4102	q(i,k) = q_in(i,kte+2-k,j)
4103	ght(i,k) = ght_in(i,kte+2-k,j)
4104	END DO
4105	END DO
4106	ELSE
4107	DO i = its , MIN (ide-1 , ite )
4108	DO k = kts,kte
4109	p(i,k) = p_in(i,k ,j)
4110	t(i,k) = t_in(i,k ,j)
4111	q(i,k) = q_in(i,k ,j)
4112	ght(i,k) = ght_in(i,k ,j)
4113	END DO
4114	END DO
4115	END IF
4116
4117	! Find the first level above the ground. If all of the levels are above ground, such as
4118	! a terrain following lower coordinate, then the first level above ground is index #2.
4119
4120	DO i = its , MIN (ide-1 , ite )
4121	level_above_sfc(i) = -1
4122	IF ( p(i,kts+1) .LT. psfc(i,j) ) THEN
4123	level_above_sfc(i) = kts+1
4124	ELSE
4125	find_k : DO k = kts+1,kte-1
4126	IF ( ( p(i,k )-psfc(i,j) .GE. 0. ) .AND. &
4127	( p(i,k+1)-psfc(i,j) .LT. 0. ) ) THEN
4128	level_above_sfc(i) = k+1
4129	EXIT find_k
4130	END IF
4131	END DO find_k
4132	IF ( level_above_sfc(i) .EQ. -1 ) THEN
4133	print *,'i,j = ',i,j
4134	print *,'p = ',p(i,:)
4135	print *,'p sfc = ',psfc(i,j)
4136	CALL wrf_error_fatal ( 'Could not find level above ground')
4137	END IF
4138	END IF
4139	END DO
4140
4141	DO i = its , MIN (ide-1 , ite )
4142
4143	! Account for the moisture above the ground.
4144
4145	pd(i,kte) = p(i,kte)
4146	DO k = kte-1,level_above_sfc(i),-1
4147	rhobar = ( p(i,k ) / ( Rd * t(i,k ) ) + &
4148	p(i,k+1) / ( Rd * t(i,k+1) ) ) * 0.5
4149	qbar = ( q(i,k ) + q(i,k+1) ) * 0.5
4150	dz = ght(i,k+1) - ght(i,k)
4151	intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
4152	pd(i,k) = p(i,k) - intq(i,j)
4153	END DO
4154
4155	! Account for the moisture between the surface and the first level up.
4156
4157	IF ( ( p(i,level_above_sfc(i)-1)-psfc(i,j) .GE. 0. ) .AND. &
4158	( p(i,level_above_sfc(i) )-psfc(i,j) .LT. 0. ) .AND. &
4159	( level_above_sfc(i) .GT. kts ) ) THEN
4160	p1 = psfc(i,j)
4161	p2 = p(i,level_above_sfc(i))
4162	t1 = tsfc(i,j)
4163	t2 = t(i,level_above_sfc(i))
4164	q1 = qsfc(i,j)
4165	q2 = q(i,level_above_sfc(i))
4166	z1 = zsfc(i,j)
4167	z2 = ght(i,level_above_sfc(i))
4168	rhobar = ( p1 / ( Rd * t1 ) + &
4169	p2 / ( Rd * t2 ) ) * 0.5
4170	qbar = ( q1 + q2 ) * 0.5
4171	dz = z2 - z1
4172	IF ( dz .GT. 0.1 ) THEN
4173	intq(i,j) = intq(i,j) + g * qbar * rhobar / (1. + qbar) * dz
4174	END IF
4175
4176	! Fix the underground values.
4177
4178	DO k = level_above_sfc(i)-1,kts+1,-1
4179	pd(i,k) = p(i,k) - intq(i,j)
4180	END DO
4181	END IF
4182	pd(i,kts) = psfc(i,j) - intq(i,j)
4183
4184	END DO
4185
4186	IF ( upside_down ) THEN
4187	DO i = its , MIN (ide-1 , ite )
4188	pd_out(i,kts,j) = pd(i,kts)
4189	DO k = kts+1,kte
4190	pd_out(i,kte+2-k,j) = pd(i,k)
4191	END DO
4192	END DO
4193	ELSE
4194	DO i = its , MIN (ide-1 , ite )
4195	DO k = kts,kte
4196	pd_out(i,k,j) = pd(i,k)
4197	END DO
4198	END DO
4199	END IF
4200
4201	END DO
4202
4203
4204	!!!****MARS: no water vapor pressure
4205	!! DO k = level_above_sfc(i)-1,kts+1,-1
4206	!! pd(i,k) = p(i,k)
4207	!! END DO
4208	!! pd(i,kts) = psfc(i,j)
4209	!!!****MARS
4210
4211
4212	END SUBROUTINE integ_moist
4213
4214	!---------------------------------------------------------------------
4215
4216	SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
4217	ids , ide , jds , jde , kds , kde , &
4218	ims , ime , jms , jme , kms , kme , &
4219	its , ite , jts , jte , kts , kte )
4220
4221	IMPLICIT NONE
4222
4223	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4224	ims , ime , jms , jme , kms , kme , &
4225	its , ite , jts , jte , kts , kte
4226
4227	LOGICAL , INTENT(IN) :: wrt_liquid
4228
4229	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
4230	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
4231	REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: q
4232
4233	! Local vars
4234
4235	INTEGER :: i , j , k
4236
4237	REAL :: ew , q1 , t1
4238	!****MARS .... regler si besoin ....
4239	!****MARS
4240	REAL, PARAMETER :: T_REF = 0.0
4241	REAL, PARAMETER :: MW_AIR = 28.966
4242	REAL, PARAMETER :: MW_VAP = 18.0152
4243
4244	REAL, PARAMETER :: A0 = 6.107799961
4245	REAL, PARAMETER :: A1 = 4.436518521e-01
4246	REAL, PARAMETER :: A2 = 1.428945805e-02
4247	REAL, PARAMETER :: A3 = 2.650648471e-04
4248	REAL, PARAMETER :: A4 = 3.031240396e-06
4249	REAL, PARAMETER :: A5 = 2.034080948e-08
4250	REAL, PARAMETER :: A6 = 6.136820929e-11
4251
4252	REAL, PARAMETER :: ES0 = 6.1121
4253
4254	REAL, PARAMETER :: C1 = 9.09718
4255	REAL, PARAMETER :: C2 = 3.56654
4256	REAL, PARAMETER :: C3 = 0.876793
4257	REAL, PARAMETER :: EIS = 6.1071
4258	REAL :: RHS
4259	REAL, PARAMETER :: TF = 273.16
4260	REAL :: TK
4261
4262	REAL :: ES
4263	REAL :: QS
4264	REAL, PARAMETER :: EPS = 0.622
4265	REAL, PARAMETER :: SVP1 = 0.6112
4266	REAL, PARAMETER :: SVP2 = 17.67
4267	REAL, PARAMETER :: SVP3 = 29.65
4268	REAL, PARAMETER :: SVPT0 = 273.15
4269	!****MARS
4270	!****MARS
4271
4272
4273	! This subroutine computes mixing ratio (q, kg/kg) from basic variables
4274	! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
4275	! The reference temperature (t_ref, C) is used to describe the temperature
4276	! at which the liquid and ice phase change occurs.
4277
4278	DO j = jts , MIN ( jde-1 , jte )
4279	DO k = kts , kte
4280	DO i = its , MIN (ide-1 , ite )
4281	rh(i,k,j) = MIN ( MAX ( rh(i,k,j) , 1. ) , 100. )
4282	END DO
4283	END DO
4284	END DO
4285
4286	IF ( wrt_liquid ) THEN
4287	DO j = jts , MIN ( jde-1 , jte )
4288	DO k = kts , kte
4289	DO i = its , MIN (ide-1 , ite )
4290	es=svp110.EXP(svp2*(t(i,k,j)-svpt0)/(t(i,k,j)-svp3))
4291	qs=eps*es/(p(i,k,j)/100.-es)
4292	q(i,k,j)=MAX(.01rh(i,k,j)qs,0.0)
4293	END DO
4294	END DO
4295	END DO
4296
4297	ELSE
4298	DO j = jts , MIN ( jde-1 , jte )
4299	DO k = kts , kte
4300	DO i = its , MIN (ide-1 , ite )
4301
4302	t1 = t(i,k,j) - 273.16
4303
4304	! Obviously dry.
4305
4306	IF ( t1 .lt. -200. ) THEN
4307	q(i,k,j) = 0
4308
4309	ELSE
4310
4311	! First compute the ambient vapor pressure of water
4312
4313	IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
4314	ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
4315
4316	ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
4317	ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
4318
4319	ELSE
4320	tk = t(i,k,j)
4321	rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
4322	c3 * (1. - tk / tf) + alog10(eis)
4323	ew = 10. ** rhs
4324
4325	END IF
4326
4327	! Now sat vap pres obtained compute local vapor pressure
4328
4329	ew = MAX ( ew , 0. ) * rh(i,k,j) * 0.01
4330
4331	! Now compute the specific humidity using the partial vapor
4332	! pressures of water vapor (ew) and dry air (p-ew). The
4333	! constants assume that the pressure is in hPa, so we divide
4334	! the pressures by 100.
4335
4336	q1 = mw_vap * ew
4337	q1 = q1 / (q1 + mw_air * (p(i,k,j)/100. - ew))
4338
4339	q(i,k,j) = q1 / (1. - q1 )
4340
4341	END IF
4342
4343	END DO
4344	END DO
4345	END DO
4346
4347	END IF
4348
4349	!!****MARS
4350	!!TODO: change once tracers are activated ?
4351	!q=0.
4352	!!****MARS
4353
4354	END SUBROUTINE rh_to_mxrat
4355
4356	!---------------------------------------------------------------------
4357
4358	SUBROUTINE compute_eta ( znw , &
4359	eta_levels , max_eta , max_dz , &
4360	fixedpbl, &
4361	p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0 , &
4362	ids , ide , jds , jde , kds , kde , &
4363	ims , ime , jms , jme , kms , kme , &
4364	its , ite , jts , jte , kts , kte )
4365
4366	! Compute eta levels, either using given values from the namelist (hardly
4367	! a computation, yep, I know), or assuming a constant dz above the PBL,
4368	! knowing p_top and the number of eta levels.
4369
4370	IMPLICIT NONE
4371
4372	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4373	ims , ime , jms , jme , kms , kme , &
4374	its , ite , jts , jte , kts , kte
4375	REAL , INTENT(IN) :: max_dz
4376	REAL , INTENT(IN) :: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0
4377	INTEGER , INTENT(IN) :: max_eta
4378	REAL , DIMENSION (max_eta) , INTENT(IN) :: eta_levels
4379
4380	REAL , DIMENSION (kts:kte) , INTENT(OUT) :: znw
4381
4382	! Local vars
4383
4384	INTEGER :: k
4385	REAL :: mub , t_init , p_surf , pb, ztop, ztop_pbl , dz , temp
4386	REAL , DIMENSION(kts:kte) :: dnw
4387
4388	INTEGER , PARAMETER :: prac_levels = 17
4389	INTEGER :: loop , loop1
4390	REAL , DIMENSION(prac_levels) :: znw_prac , znu_prac , dnw_prac
4391	REAL , DIMENSION(kts:kte) :: alb , phb
4392
4393
4394	!****MARS
4395	!****MARS
4396	INTEGER, fixedpbl ! usually, 8 first layers are fixed
4397	! change this parameter if the top is very
4398	! low
4399	print *, 'check Mars: p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0'
4400	print *, p_top , g , p00 , cvpm , a , r_d , cp , t00 , p1000mb , t0
4401	!-----solution alternative: definir dans la namelist les niveaux verticaux
4402	!****MARS
4403	!****MARS
4404
4405
4406	! Gee, do the eta levels come in from the namelist?
4407
4408	IF ( ABS(eta_levels(1)+1.) .GT. 0.0000001 ) THEN
4409
4410	IF ( ( ABS(eta_levels(1 )-1.) .LT. 0.0000001 ) .AND. &
4411	( ABS(eta_levels(kde)-0.) .LT. 0.0000001 ) ) THEN
4412	DO k = kds+1 , kde-1
4413	znw(k) = eta_levels(k)
4414	END DO
4415	znw( 1) = 1.
4416	znw(kde) = 0.
4417	ELSE
4418	CALL wrf_error_fatal ( 'First eta level should be 1.0 and the last 0.0 in namelist' )
4419	END IF
4420
4421	! Compute eta levels assuming a constant delta z above the PBL.
4422
4423	ELSE
4424
4425	! Compute top of the atmosphere with some silly levels. We just want to
4426	! integrate to get a reasonable value for ztop. We use the planned PBL-esque
4427	! levels, and then just coarse resolution above that. We know p_top, and we
4428	! have the base state vars.
4429
4430	p_surf = p00
4431
4432	! znw_prac = (/ 1.000 , 0.993 , 0.983 , 0.970 , 0.954 , 0.934 , 0.909 , &
4433	! 0.88 , 0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
4434
4435	!****MARS
4436	!****MARS
4437	! on Mars, this is important to correctly resolve the surface
4438	! -- levels were changed to get closer to the surface
4439	! -- values were chosen as done typically in LMD GCM simulations
4440	!TODO: better repartition ?
4441	znw_prac = (/ 1.000 , &
4442	0.9995 , &
4443	0.9980 , &
4444	0.9950 , &
4445	0.9850 , &
4446	0.9700 , &
4447	0.9400 , &
4448	0.9000 , &
4449	0.8 , 0.7 , 0.6 , 0.5 , 0.4 , 0.3 , 0.2 , 0.1 , 0.0 /)
4450	!****MARS
4451	!****MARS
4452
4453
4454	DO k = 1 , prac_levels - 1
4455	znu_prac(k) = ( znw_prac(k) + znw_prac(k+1) ) * 0.5
4456	dnw_prac(k) = znw_prac(k+1) - znw_prac(k)
4457	END DO
4458
4459	DO k = 1, prac_levels-1
4460	pb = znu_prac(k)*(p_surf - p_top) + p_top
4461	! temp = MAX ( 200., t00 + A*LOG(pb/p00) )
4462	temp = t00 + A*LOG(pb/p00)
4463	t_init = temp(p00/pb)*(r_d/cp) - t0
4464	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
4465	END DO
4466
4467	! Base state mu is defined as base state surface pressure minus p_top
4468
4469	mub = p_surf - p_top
4470
4471	! Integrate base geopotential, starting at terrain elevation.
4472
4473	phb(1) = 0.
4474	DO k = 2,prac_levels
4475	phb(k) = phb(k-1) - dnw_prac(k-1)mubalb(k-1)
4476	END DO
4477
4478	! So, now we know the model top in meters. Get the average depth above the PBL
4479	! of each of the remaining levels. We are going for a constant delta z thickness.
4480
4481	ztop = phb(prac_levels) / g
4482	ztop_pbl = phb(fixedpbl) / g
4483	dz = ( ztop - ztop_pbl ) / REAL ( kde - fixedpbl )
4484
4485	! Standard levels near the surface so no one gets in trouble.
4486	DO k = 1 , fixedpbl
4487	znw(k) = znw_prac(k)
4488	END DO
4489
4490	! Using d phb(k)/ d eta(k) = -mub * alb(k), eqn 2.9
4491	! Skamarock et al, NCAR TN 468. Use full levels, so
4492	! use twice the thickness.
4493
4494	DO k = fixedpbl, kte-1
4495	pb = znw(k) * (p_surf - p_top) + p_top
4496	! temp = MAX ( 200., t00 + A*LOG(pb/p00) )
4497	temp = t00 + A*LOG(pb/p00)
4498	t_init = temp(p00/pb)*(r_d/cp) - t0
4499	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
4500	znw(k+1) = znw(k) - dzg / ( mubalb(k) )
4501	END DO
4502	znw(kte) = 0.000
4503
4504	! There is some iteration. We want the top level, ztop, to be
4505	! consistent with the delta z, and we want the half level values
4506	! to be consistent with the eta levels. The inner loop to 10 gets
4507	! the eta levels very accurately, but has a residual at the top, due
4508	! to dz changing. We reset dz five times, and then things seem OK.
4509
4510
4511	DO loop1 = 1 , 5
4512	DO loop = 1 , 10
4513	DO k = fixedpbl, kte-1
4514	pb = (znw(k)+znw(k+1))0.5 (p_surf - p_top) + p_top
4515	! temp = MAX ( 200., t00 + A*LOG(pb/p00) )
4516	temp = t00 + A*LOG(pb/p00)
4517	t_init = temp(p00/pb)*(r_d/cp) - t0
4518	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
4519	znw(k+1) = znw(k) - dzg / ( mubalb(k) )
4520	!!****MARS
4521	!!attention 'base_lapse' ne doit pas etre trop grand
4522	!!sinon ... des NaN car temperatures negatives en haut
4523	!IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
4524	! IF (k .EQ. 8) THEN
4525	! print *, 'p,t,z,k'
4526	! END IF
4527	! print *, pb,temp,znw(k+1),k
4528	!END IF
4529	!****MARS
4530	END DO
4531	IF ( ( loop1 .EQ. 5 ) .AND. ( loop .EQ. 10 ) ) THEN
4532	print *,'Converged znw(kte) should be 0.0 = ',znw(kte)
4533	END IF
4534	znw(kte) = 0.000
4535	END DO
4536
4537	! Here is where we check the eta levels values we just computed.
4538
4539	DO k = 1, kde-1
4540	pb = (znw(k)+znw(k+1))0.5 (p_surf - p_top) + p_top
4541	! temp = MAX ( 200., t00 + A*LOG(pb/p00) )
4542	temp = t00 + A*LOG(pb/p00)
4543	t_init = temp(p00/pb)*(r_d/cp) - t0
4544	alb(k) = (r_d/p1000mb)(t_init+t0)(pb/p1000mb)**cvpm
4545	END DO
4546
4547	phb(1) = 0.
4548	DO k = 2,kde
4549	phb(k) = phb(k-1) - (znw(k)-znw(k-1)) * mub*alb(k-1)
4550	END DO
4551
4552	! Reset the model top and the dz, and iterate.
4553
4554	ztop = phb(kde)/g
4555	ztop_pbl = phb(fixedpbl)/g
4556	dz = ( ztop - ztop_pbl ) / REAL ( kde - fixedpbl )
4557	END DO
4558
4559
4560	! ****MARS
4561	! Display the computed levels
4562	print *,'WRF levels are:'
4563	print *,'z (m) = ',phb(1)/g
4564	do k = 2 ,kte
4565	print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
4566
4567
4568	!! little check of the repartition
4569	if (k>2) then
4570	if ((phb(k)-2.*phb(k-1)+phb(k-2))/g < -1.e-2) then
4571	print *, 'problem on the repartition'
4572	print *, '>> try to decrease force_sfc_in_vinterp (<8)'
4573	print *, '>> or increase model top (i.e. lower ptop)'
4574	print , (phb(k)-2.phb(k-1)+phb(k-2))/g
4575	stop
4576	endif
4577	endif
4578	end do
4579	! ****MARS
4580
4581
4582	IF ( dz .GT. max_dz ) THEN
4583	print *,'z (m) = ',phb(1)/g
4584	do k = 2 ,kte
4585	print *,'z (m) and dz (m) = ',phb(k)/g,(phb(k)-phb(k-1))/g
4586	end do
4587	print *,'dz (m) above fixed eta levels = ',dz
4588	print *,'namelist max_dz (m) = ',max_dz
4589	print *,'namelist p_top (Pa) = ',p_top
4590	CALL wrf_debug ( 0, 'You need one of three things:' )
4591	CALL wrf_debug ( 0, '1) More eta levels to reduce the dz: e_vert' )
4592	CALL wrf_debug ( 0, '2) A lower p_top so your total height is reduced: p_top_requested')
4593	CALL wrf_debug ( 0, '3) Increase the maximum allowable eta thickness: max_dz')
4594	CALL wrf_debug ( 0, 'All are namelist options')
4595	CALL wrf_error_fatal ( 'dz above fixed eta levels is too large')
4596	END IF
4597
4598	END IF
4599
4600	END SUBROUTINE compute_eta
4601
4602	!---------------------------------------------------------------------
4603
4604	SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
4605	ids , ide , jds , jde , kds , kde , &
4606	ims , ime , jms , jme , kms , kme , &
4607	its , ite , jts , jte , kts , kte )
4608
4609	! Plow through each month, find the max, min values for each i,j.
4610
4611	IMPLICIT NONE
4612
4613	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4614	ims , ime , jms , jme , kms , kme , &
4615	its , ite , jts , jte , kts , kte
4616
4617	REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
4618	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
4619
4620	! Local vars
4621
4622	INTEGER :: i , j , l
4623	REAL :: minner , maxxer
4624
4625	DO j = jts , MIN(jde-1,jte)
4626	DO i = its , MIN(ide-1,ite)
4627	minner = field_in(i,1,j)
4628	maxxer = field_in(i,1,j)
4629	DO l = 2 , 12
4630	IF ( field_in(i,l,j) .LT. minner ) THEN
4631	minner = field_in(i,l,j)
4632	END IF
4633	IF ( field_in(i,l,j) .GT. maxxer ) THEN
4634	maxxer = field_in(i,l,j)
4635	END IF
4636	END DO
4637	field_min(i,j) = minner
4638	field_max(i,j) = maxxer
4639	END DO
4640	END DO
4641
4642	END SUBROUTINE monthly_min_max
4643
4644	!---------------------------------------------------------------------
4645
4646	SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
4647	ids , ide , jds , jde , kds , kde , &
4648	ims , ime , jms , jme , kms , kme , &
4649	its , ite , jts , jte , kts , kte )
4650
4651	! Linrarly in time interpolate data to a current valid time. The data is
4652	! assumed to come in "monthly", valid at the 15th of every month.
4653
4654	IMPLICIT NONE
4655
4656	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4657	ims , ime , jms , jme , kms , kme , &
4658	its , ite , jts , jte , kts , kte
4659
4660	CHARACTER (LEN=24) , INTENT(IN) :: date_str
4661	REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in
4662	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
4663
4664	! Local vars
4665
4666	INTEGER :: i , j , l
4667	INTEGER , DIMENSION(0:13) :: middle
4668	INTEGER :: target_julyr , target_julday , target_date
4669	INTEGER :: julyr , julday , int_month , month1 , month2
4670	REAL :: gmt
4671	CHARACTER (LEN=4) :: yr
4672	CHARACTER (LEN=2) :: mon , day15
4673
4674
4675	WRITE(day15,FMT='(I2.2)') 15
4676	DO l = 1 , 12
4677	WRITE(mon,FMT='(I2.2)') l
4678	CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
4679	middle(l) = julyr*1000 + julday
4680	END DO
4681
4682	l = 0
4683	middle(l) = middle( 1) - 31
4684
4685	l = 13
4686	middle(l) = middle(12) + 31
4687
4688	CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
4689	target_date = target_julyr * 1000 + target_julday
4690	find_month : DO l = 0 , 12
4691	IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
4692	DO j = jts , MIN ( jde-1 , jte )
4693	DO i = its , MIN (ide-1 , ite )
4694	int_month = l
4695	IF ( ( int_month .EQ. 0 ) .OR. ( int_month .EQ. 12 ) ) THEN
4696	month1 = 12
4697	month2 = 1
4698	ELSE
4699	month1 = int_month
4700	month2 = month1 + 1
4701	END IF
4702	field_out(i,j) = ( field_in(i,month2,j) * ( target_date - middle(l) ) + &
4703	field_in(i,month1,j) * ( middle(l+1) - target_date ) ) / &
4704	( middle(l+1) - middle(l) )
4705	END DO
4706	END DO
4707	EXIT find_month
4708	END IF
4709	END DO find_month
4710
4711	END SUBROUTINE monthly_interp_to_date
4712
4713	!---------------------------------------------------------------------
4714
4715	SUBROUTINE sfcprs (t, q, height, pslv, ter, avgsfct, p, &
4716	psfc, ez_method, &
4717	ids , ide , jds , jde , kds , kde , &
4718	ims , ime , jms , jme , kms , kme , &
4719	its , ite , jts , jte , kts , kte )
4720
4721
4722	! Computes the surface pressure using the input height,
4723	! temperature and q (already computed from relative
4724	! humidity) on p surfaces. Sea level pressure is used
4725	! to extrapolate a first guess.
4726
4727	IMPLICIT NONE
4728
4729	!****MARS
4730	REAL , PARAMETER :: Rd = 192.
4731	REAL , PARAMETER :: Cp = 844.6
4732	REAL, PARAMETER :: g = 3.72
4733	REAL, PARAMETER :: pconst = 610.
4734	!****MARS
4735
4736	!****MARS .... to be changed if used
4737	REAL, PARAMETER :: gamma = 6.5E-3
4738	REAL, PARAMETER :: TC = 273.15 + 17.5
4739	REAL, PARAMETER :: gammarg = gamma * Rd / g
4740	REAL, PARAMETER :: rov2 = Rd / 2.
4741	!****MARS .... to be changed if used
4742
4743	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
4744	ims , ime , jms , jme , kms , kme , &
4745	its , ite , jts , jte , kts , kte
4746	LOGICAL , INTENT ( IN ) :: ez_method
4747
4748	REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
4749	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: pslv , ter, avgsfct
4750	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
4751
4752	INTEGER :: i
4753	INTEGER :: j
4754	INTEGER :: k
4755	INTEGER , DIMENSION (its:ite,jts:jte) :: k500 , k700 , k850
4756
4757	LOGICAL :: l1
4758	LOGICAL :: l2
4759	LOGICAL :: l3
4760	LOGICAL :: OK
4761
4762	REAL :: gamma78 ( its:ite,jts:jte )
4763	REAL :: gamma57 ( its:ite,jts:jte )
4764	REAL :: ht ( its:ite,jts:jte )
4765	REAL :: p1 ( its:ite,jts:jte )
4766	REAL :: t1 ( its:ite,jts:jte )
4767	REAL :: t500 ( its:ite,jts:jte )
4768	REAL :: t700 ( its:ite,jts:jte )
4769	REAL :: t850 ( its:ite,jts:jte )
4770	REAL :: tfixed ( its:ite,jts:jte )
4771	REAL :: tsfc ( its:ite,jts:jte )
4772	REAL :: tslv ( its:ite,jts:jte )
4773
4774	! We either compute the surface pressure from a time averaged surface temperature
4775	! (what we will call the "easy way"), or we try to remove the diurnal impact on the
4776	! surface temperature (what we will call the "other way"). Both are essentially
4777	! corrections to a sea level pressure with a high-resolution topography field.
4778
4779	!****MARS ....
4780	!****MARS .... the mean sea level method is abandoned
4781	print *, 'no sea level pressure on Mars, please'
4782	stop
4783	!****MARS ....
4784
4785	IF ( ez_method ) THEN
4786
4787	DO j = jts , MIN(jde-1,jte)
4788	DO i = its , MIN(ide-1,ite)
4789	psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / avgsfct(i,j) ) ** ( - g / ( Rd * gamma ) )
4790	END DO
4791	END DO
4792
4793	ELSE
4794
4795	! Find the locations of the 850, 700 and 500 mb levels.
4796
4797	k850 = 0 ! find k at: P=850
4798	k700 = 0 ! P=700
4799	k500 = 0 ! P=500
4800
4801	i = its
4802	j = jts
4803	DO k = kts+1 , kte
4804	IF (NINT(p(i,k,j)) .EQ. 85000) THEN
4805	k850(i,j) = k
4806	ELSE IF (NINT(p(i,k,j)) .EQ. 70000) THEN
4807	k700(i,j) = k
4808	ELSE IF (NINT(p(i,k,j)) .EQ. 50000) THEN
4809	k500(i,j) = k
4810	END IF
4811	END DO
4812
4813	IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
4814
4815	DO j = jts , MIN(jde-1,jte)
4816	DO i = its , MIN(ide-1,ite)
4817	psfc(i,j) = pslv(i,j) * ( 1.0 + gamma * ter(i,j) / t(i,1,j) ) ** ( - g / ( Rd * gamma ) )
4818	END DO
4819	END DO
4820
4821	RETURN
4822	#if 0
4823
4824	! Possibly it is just that we have a generalized vertical coord, so we do not
4825	! have the values exactly. Do a simple assignment to a close vertical level.
4826
4827	DO j = jts , MIN(jde-1,jte)
4828	DO i = its , MIN(ide-1,ite)
4829	DO k = kts+1 , kte-1
4830	IF ( ( p(i,k,j) - 85000. ) * ( p(i,k+1,j) - 85000. ) .LE. 0.0 ) THEN
4831	k850(i,j) = k
4832	END IF
4833	IF ( ( p(i,k,j) - 70000. ) * ( p(i,k+1,j) - 70000. ) .LE. 0.0 ) THEN
4834	k700(i,j) = k
4835	END IF
4836	IF ( ( p(i,k,j) - 50000. ) * ( p(i,k+1,j) - 50000. ) .LE. 0.0 ) THEN
4837	k500(i,j) = k
4838	END IF
4839	END DO
4840	END DO
4841	END DO
4842
4843	! If we still do not have the k levels, punt. I mean, we did try.
4844
4845	OK = .TRUE.
4846	DO j = jts , MIN(jde-1,jte)
4847	DO i = its , MIN(ide-1,ite)
4848	IF ( ( k850(i,j) .EQ. 0 ) .OR. ( k700(i,j) .EQ. 0 ) .OR. ( k500(i,j) .EQ. 0 ) ) THEN
4849	OK = .FALSE.
4850	PRINT '(A)','(i,j) = ',i,j,' Error in finding p level for 850, 700 or 500 hPa.'
4851	DO K = kts+1 , kte
4852	PRINT '(A,I3,A,F10.2,A)','K = ',k,' PRESSURE = ',p(i,k,j),' Pa'
4853	END DO
4854	PRINT '(A)','Expected 850, 700, and 500 mb values, at least.'
4855	END IF
4856	END DO
4857	END DO
4858	IF ( .NOT. OK ) THEN
4859	CALL wrf_error_fatal ( 'wrong pressure levels' )
4860	END IF
4861	#endif
4862
4863	! We are here if the data is isobaric and we found the levels for 850, 700,
4864	! and 500 mb right off the bat.
4865
4866	ELSE
4867	DO j = jts , MIN(jde-1,jte)
4868	DO i = its , MIN(ide-1,ite)
4869	k850(i,j) = k850(its,jts)
4870	k700(i,j) = k700(its,jts)
4871	k500(i,j) = k500(its,jts)
4872	END DO
4873	END DO
4874	END IF
4875
4876	! The 850 hPa level of geopotential height is called something special.
4877
4878	DO j = jts , MIN(jde-1,jte)
4879	DO i = its , MIN(ide-1,ite)
4880	ht(i,j) = height(i,k850(i,j),j)
4881	END DO
4882	END DO
4883
4884	! The variable ht is now -ter/ht(850 hPa). The plot thickens.
4885
4886	DO j = jts , MIN(jde-1,jte)
4887	DO i = its , MIN(ide-1,ite)
4888	ht(i,j) = -ter(i,j) / ht(i,j)
4889	END DO
4890	END DO
4891
4892	! Make an isothermal assumption to get a first guess at the surface
4893	! pressure. This is to tell us which levels to use for the lapse
4894	! rates in a bit.
4895
4896	DO j = jts , MIN(jde-1,jte)
4897	DO i = its , MIN(ide-1,ite)
4898	psfc(i,j) = pslv(i,j) * (pslv(i,j) / p(i,k850(i,j),j)) ** ht(i,j)
4899	END DO
4900	END DO
4901
4902	! Get a pressure more than pconst Pa above the surface - p1. The
4903	! p1 is the top of the level that we will use for our lapse rate
4904	! computations.
4905
4906	DO j = jts , MIN(jde-1,jte)
4907	DO i = its , MIN(ide-1,ite)
4908	IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
4909	p1(i,j) = 85000.
4910	ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0. ) THEN
4911	p1(i,j) = psfc(i,j) - pconst
4912	ELSE
4913	p1(i,j) = 50000.
4914	END IF
4915	END DO
4916	END DO
4917
4918	! Compute virtual temperatures for k850, k700, and k500 layers. Now
4919	! you see why we wanted Q on pressure levels, it all is beginning
4920	! to make sense.
4921
4922	DO j = jts , MIN(jde-1,jte)
4923	DO i = its , MIN(ide-1,ite)
4924	t850(i,j) = t(i,k850(i,j),j) * (1. + 0.608 * q(i,k850(i,j),j))
4925	t700(i,j) = t(i,k700(i,j),j) * (1. + 0.608 * q(i,k700(i,j),j))
4926	t500(i,j) = t(i,k500(i,j),j) * (1. + 0.608 * q(i,k500(i,j),j))
4927	END DO
4928	END DO
4929
4930	! Compute lapse rates between these three levels. These are
4931	! environmental values for each (i,j).
4932
4933	DO j = jts , MIN(jde-1,jte)
4934	DO i = its , MIN(ide-1,ite)
4935	gamma78(i,j) = ALOG(t850(i,j) / t700(i,j)) / ALOG (p(i,k850(i,j),j) / p(i,k700(i,j),j) )
4936	gamma57(i,j) = ALOG(t700(i,j) / t500(i,j)) / ALOG (p(i,k700(i,j),j) / p(i,k500(i,j),j) )
4937	END DO
4938	END DO
4939
4940	DO j = jts , MIN(jde-1,jte)
4941	DO i = its , MIN(ide-1,ite)
4942	IF ( ( psfc(i,j) - 95000. ) .GE. 0. ) THEN
4943	t1(i,j) = t850(i,j)
4944	ELSE IF ( ( psfc(i,j) - 85000. ) .GE. 0. ) THEN
4945	t1(i,j) = t700(i,j) * (p1(i,j) / (p(i,k700(i,j),j))) ** gamma78(i,j)
4946	ELSE IF ( ( psfc(i,j) - 70000. ) .GE. 0.) THEN
4947	t1(i,j) = t500(i,j) * (p1(i,j) / (p(i,k500(i,j),j))) ** gamma57(i,j)
4948	ELSE
4949	t1(i,j) = t500(i,j)
4950	ENDIF
4951	END DO
4952	END DO
4953
4954	! From our temperature way up in the air, we extrapolate down to
4955	! the sea level to get a guess at the sea level temperature.
4956
4957	DO j = jts , MIN(jde-1,jte)
4958	DO i = its , MIN(ide-1,ite)
4959	tslv(i,j) = t1(i,j) * (pslv(i,j) / p1(i,j)) ** gammarg
4960	END DO
4961	END DO
4962
4963	! The new surface temperature is computed from the with new sea level
4964	! temperature, just using the elevation and a lapse rate. This lapse
4965	! rate is -6.5 K/km.
4966
4967	DO j = jts , MIN(jde-1,jte)
4968	DO i = its , MIN(ide-1,ite)
4969	tsfc(i,j) = tslv(i,j) - gamma * ter(i,j)
4970	END DO
4971	END DO
4972
4973	! A small correction to the sea-level temperature, in case it is too warm.
4974
4975	DO j = jts , MIN(jde-1,jte)
4976	DO i = its , MIN(ide-1,ite)
4977	tfixed(i,j) = tc - 0.005 * (tsfc(i,j) - tc) ** 2
4978	END DO
4979	END DO
4980
4981	DO j = jts , MIN(jde-1,jte)
4982	DO i = its , MIN(ide-1,ite)
4983	l1 = tslv(i,j) .LT. tc
4984	l2 = tsfc(i,j) .LE. tc
4985	l3 = .NOT. l1
4986	IF ( l2 .AND. l3 ) THEN
4987	tslv(i,j) = tc
4988	ELSE IF ( ( .NOT. l2 ) .AND. l3 ) THEN
4989	tslv(i,j) = tfixed(i,j)
4990	END IF
4991	END DO
4992	END DO
4993
4994	! Finally, we can get to the surface pressure.
4995
4996	DO j = jts , MIN(jde-1,jte)
4997	DO i = its , MIN(ide-1,ite)
4998	p1(i,j) = - ter(i,j) * g / ( rov2 * ( tsfc(i,j) + tslv(i,j) ) )
4999	psfc(i,j) = pslv(i,j) * EXP ( p1(i,j) )
5000	END DO
5001	END DO
5002
5003	END IF
5004
5005	! Surface pressure and sea-level pressure are the same at sea level.
5006
5007	! DO j = jts , MIN(jde-1,jte)
5008	! DO i = its , MIN(ide-1,ite)
5009	! IF ( ABS ( ter(i,j) ) .LT. 0.1 ) THEN
5010	! psfc(i,j) = pslv(i,j)
5011	! END IF
5012	! END DO
5013	! END DO
5014
5015	END SUBROUTINE sfcprs
5016
5017	!---------------------------------------------------------------------
5018
5019	SUBROUTINE sfcprs2(t, q, height, psfc_in, ter, avgsfct, p, &
5020	psfc, ez_method, &
5021	ids , ide , jds , jde , kds , kde , &
5022	ims , ime , jms , jme , kms , kme , &
5023	its , ite , jts , jte , kts , kte )
5024
5025
5026	! Computes the surface pressure using the input height,
5027	! temperature and q (already computed from relative
5028	! humidity) on p surfaces. Sea level pressure is used
5029	! to extrapolate a first guess.
5030
5031	IMPLICIT NONE
5032
5033	!****MARS
5034	REAL , PARAMETER :: Rd = 192.
5035	REAL, PARAMETER :: g = 3.72
5036	!****MARS
5037
5038	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5039	ims , ime , jms , jme , kms , kme , &
5040	its , ite , jts , jte , kts , kte
5041	LOGICAL , INTENT ( IN ) :: ez_method
5042
5043	REAL , DIMENSION (ims:ime,kms:kme,jms:jme) , INTENT(IN ):: t, q, height, p
5044	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(IN ):: psfc_in , ter, avgsfct
5045	REAL , DIMENSION (ims:ime, jms:jme) , INTENT(OUT):: psfc
5046
5047	INTEGER :: i
5048	INTEGER :: j
5049	INTEGER :: k
5050
5051	REAL :: tv_sfc_avg , tv_sfc , del_z
5052
5053	! Compute the new surface pressure from the old surface pressure, and a
5054	! known change in elevation at the surface.
5055
5056
5057	!****MARS: as is done in MCD/pres0 with the MOLA topography :)
5058
5059	!!---------
5060	!! del_z = diff in surface topo, lo-res vs hi-res
5061	!grid%em_ght_gc - grid%ht
5062	!!---------
5063	!!* em_ght_gc: surface geopotential height from the GCM
5064	!!* ht: hi-res altimetry
5065	! psfc = psfc_in * exp ( g del_z / (Rd Tv_sfc ) )
5066	!!---------
5067
5068
5069	IF ( ez_method ) THEN
5070	!!
5071	!!****MARS: 'ez_method' is 'we_have_tavgsfc', hard-coded as false
5072	!!
5073	DO j = jts , MIN(jde-1,jte)
5074	DO i = its , MIN(ide-1,ite)
5075	tv_sfc_avg = avgsfct(i,j) * (1. + 0.608 * q(i,1,j))
5076	del_z = height(i,1,j) - ter(i,j)
5077	psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc_avg ) )
5078	END DO
5079	END DO
5080	ELSE
5081	!!
5082	!!****MARS .... here is what is done for Mars
5083	!!
5084	DO j = jts , MIN(jde-1,jte)
5085	DO i = its , MIN(ide-1,ite)
5086	! tv_sfc = t(i,1,j) * (1. + 0.608 * q(i,1,j))
5087	!!****MARS: 0.608 >> nonsense on Mars
5088	tv_sfc = t(i,1,j)
5089	!!****MARS .... changer pour t_1km - 7e couche GCM
5090	!!****MARS .... spiga et al. (2007)
5091	tv_sfc = t(i,8,j)
5092	del_z = height(i,1,j) - ter(i,j)
5093	psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * tv_sfc ) )
5094	!****MARS
5095	!****MARS .... which temperature is used in the Laplace formula ?
5096	!!****MARS: hardcoded as 220K (t0)
5097	!!****MARS: pas une enorme influence
5098	!psfc(i,j) = psfc_in(i,j) * EXP ( g * del_z / ( Rd * 220 ) )
5099
5100
5101	! !****MARS .... check of the altimetry differences
5102	! print *,del_z, tv_sfc
5103
5104	END DO
5105	END DO
5106	print *, '1 km temperatures - max'
5107	print *, MAXVAL(t(:,8,:))
5108	END IF
5109
5110	END SUBROUTINE sfcprs2
5111
5112	!---------------------------------------------------------------------
5113
5114	SUBROUTINE init_module_initialize
5115	END SUBROUTINE init_module_initialize
5116
5117	!---------------------------------------------------------------------
5118	SUBROUTINE constante3(field, field_custom, &
5119	ids , ide , jds , jde , kds , kde , &
5120	ims , ime , jms , jme , kms , kme , &
5121	its , ite , jts , jte , kts , kte )
5122
5123
5124	IMPLICIT NONE
5125
5126	REAL :: field_custom
5127	REAL, DIMENSION (ims:ime,kms:kme,jms:jme), INTENT(INOUT):: field
5128	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5129	ims , ime , jms , jme , kms , kme , &
5130	its , ite , jts , jte , kts , kte
5131
5132
5133	!!****MARS: set the 3D field to a constant value
5134	field(:,:,:)=field_custom
5135
5136	END SUBROUTINE constante3
5137	!---------------------------------------------------------------------
5138	SUBROUTINE constante2(field, field_custom, &
5139	ids , ide , jds , jde , kds , kde , &
5140	ims , ime , jms , jme , kms , kme , &
5141	its , ite , jts , jte , kts , kte )
5142
5143
5144	IMPLICIT NONE
5145
5146	REAL :: field_custom
5147	REAL, DIMENSION (ims:ime,jms:jme), INTENT(INOUT):: field
5148	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
5149	ims , ime , jms , jme , kms , kme , &
5150	its , ite , jts , jte , kts , kte
5151
5152
5153	!!****MARS: set the 3D field to a constant value
5154	field(:,:)=field_custom
5155
5156	END SUBROUTINE constante2
5157	!---------------------------------------------------------------------
5158
5159	subroutine build_sigma_hr(dimlevs,sigma_gcm,ps_gcm,ps_hr,sigma_hr) !,p_pgcm)
5160
5161	implicit none
5162	! include "constants_mcd.inc"
5163
5164	!---------------------------------------
5165	! written by E. Millour and F. Forget
5166	! Mars Climate Database v4.2
5167	! see DDD page 27 and following
5168	!---------------------------------------
5169
5170	INTEGER , INTENT(IN) :: dimlevs
5171
5172	! inputs
5173	real sigma_gcm(dimlevs) ! GCM sigma levels
5174	real ps_gcm ! GCM surface pressure
5175	real ps_hr ! High res surface pressure
5176	! outputs
5177	real sigma_hr(dimlevs) ! High res sigma levels
5178	! real p_pgcm(dimlevs) ! high res to GCM pressure ratios
5179
5180	! local variables
5181	integer l
5182	real x ! lower layer compression (-0.9<x<0) or dilatation (0.<x<0.9)
5183	real rp ! surface pressure ratio ps_hr/ps_gcm
5184	real deltaz ! corresponding pseudo-altitude difference (km)
5185	real f ! coefficient f= p_hr / p_gcm
5186	real z ! altitude of transition of p_hr toward p_gcm (km)
5187	real p_pgcm(dimlevs) ! high res to GCM pressure ratios
5188
5189	! 1. Coefficients
5190	rp=ps_hr/ps_gcm
5191	deltaz=-10.*log(rp)
5192	x = min(max(0.12*(abs(deltaz)-1.),0.),0.8)
5193	if(deltaz.gt.0) x=-x
5194	z=max(deltaz + 3.,3.)
5195
5196	do l=1,dimlevs
5197	f=rpsigma_gcm(l)*x
5198	! f=f+(1-f)0.5(1+tanh(6.(-10.log(sigma_gcm(l))-z)/z))
5199	! sigma_hr(l)=f*sigma_gcm(l)/rp
5200	p_pgcm(l)=f+(1-f)0.5(1+tanh(6.(-10.log(sigma_gcm(l))-z)/z))
5201	sigma_hr(l)=p_pgcm(l)*sigma_gcm(l)/rp
5202	enddo
5203
5204	end subroutine build_sigma_hr
5205
5206
5207
5208
5209
5210	END MODULE module_initialize
5211
5212	#endif

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

Download in other formats:

Original Format