Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real.F @ 357

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

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

Download in other formats:

Original Format