Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real_method1_2interp.F @ 3531

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

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

Download in other formats:

Original Format