Context Navigation

module_initialize_real.F @ 3026

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

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

Download in other formats:

Original Format