Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_real.F @ 2759

Last change on this file since 2759 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 150.0 KB

Line
1	!REAL:MODEL_LAYER:INITIALIZATION
2
3	! This MODULE holds the routines which are used to perform various initializations
4	! for the individual domains, specifically for the Eulerian, mass-based coordinate.
5
6	!-----------------------------------------------------------------------
7
8	MODULE module_initialize_real
9
10	USE module_bc
11	USE module_configure
12	USE module_domain
13	USE module_io_domain
14	USE module_model_constants
15	! USE module_si_io_nmm
16	USE module_state_description
17	USE module_timing
18	USE module_soil_pre
19	#ifdef DM_PARALLEL
20	USE module_dm
21	USE module_ext_internal
22	#endif
23
24	INTEGER :: internal_time_loop
25	INTEGER:: MPI_COMM_COMP,MYPE
26
27	CONTAINS
28
29	!-------------------------------------------------------------------
30
31	SUBROUTINE init_domain ( grid )
32
33	IMPLICIT NONE
34
35	! Input space and data. No gridded meteorological data has been stored, though.
36
37	! TYPE (domain), POINTER :: grid
38	TYPE (domain) :: grid
39
40	! Local data.
41
42	INTEGER :: idum1, idum2
43
44	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
45
46	CALL init_domain_nmm (grid &
47	!
48	#include <actual_args.inc>
49	!
50	)
51
52	END SUBROUTINE init_domain
53
54	!-------------------------------------------------------------------
55	!---------------------------------------------------------------------
56	SUBROUTINE init_domain_nmm ( grid &
57	!
58	# include <dummy_args.inc>
59	!
60	)
61
62	USE module_optional_input
63	IMPLICIT NONE
64
65	! Input space and data. No gridded meteorological data has been stored, though.
66
67	! TYPE (domain), POINTER :: grid
68	TYPE (domain) :: grid
69
70	# include <dummy_decl.inc>
71
72	TYPE (grid_config_rec_type) :: config_flags
73
74	! Local domain indices and counters.
75
76	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
77
78	INTEGER :: &
79	ids, ide, jds, jde, kds, kde, &
80	ims, ime, jms, jme, kms, kme, &
81	its, ite, jts, jte, kts, kte, &
82	ips, ipe, jps, jpe, kps, kpe, &
83	i, j, k, NNXP, NNYP
84
85	! Local data
86
87	CHARACTER(LEN=19):: start_date
88
89	#ifdef DM_PARALLEL
90
91	LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
92
93	! INTEGER :: DOMDESC
94	REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
95	INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
96	INTEGER, ALLOCATABLE:: ITEMP(:,:)
97	INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe
98	INTEGER :: istat,inpes,jnpes
99	#endif
100
101
102	CHARACTER (LEN=255) :: message
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, rmin,rmax
110
111	REAL :: p_top_requested , ptsgm
112	INTEGER :: num_metgrid_levels, ICOUNT
113	REAL , DIMENSION(max_eta) :: eta_levels
114
115	LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor
116
117	REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, &
118	PDVP,PSFC_OUTV
119
120	REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, &
121	QTMP,QTMP2
122
123	INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
124	KHLA,KHHA,KVLA,KVHA
125
126	! INTEGER, ALLOCATABLE, DIMENSION(:,:):: LU_INDEX
127
128	REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
129	FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
130	HDACJ,DDMPUJ,DDMPVJ
131
132	REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, &
133	SGML1,SGML2
134
135	!-- Carsel and Parrish [1988]
136	REAL , DIMENSION(100) :: lqmi
137	integer iicount
138
139	REAL:: TPH0D,TLM0D
140	REAL:: TPH0,WB,SB,TDLM,TDPH
141	REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
142	REAL:: TSPH,DTAD,DTCF
143	REAL:: ACDT,CDDAMP,DXP,FP
144	REAL:: WBD,SBD
145	REAL:: RSNOW,SNOFAC
146	REAL, PARAMETER:: SALP=2.60
147	REAL, PARAMETER:: SNUP=0.040
148	REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
149	REAL:: cur_smc, aposs_smc
150
151	INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300)
152	REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH
153
154	INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
155	INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook
156
157	REAL, PARAMETER:: DTR=0.01745329
158	REAL, PARAMETER:: W_NMM=0.08
159	REAL, PARAMETER:: COAC=1.6
160	REAL, PARAMETER:: CODAMP=6.4
161	REAL, PARAMETER:: TWOM=.00014584
162	REAL, PARAMETER:: CP=1004.6
163	REAL, PARAMETER:: DFC=1.0
164	REAL, PARAMETER:: DDFC=8.0
165	REAL, PARAMETER:: ROI=916.6
166	REAL, PARAMETER:: R=287.04
167	REAL, PARAMETER:: CI=2060.0
168	REAL, PARAMETER:: ROS=1500.
169	REAL, PARAMETER:: CS=1339.2
170	REAL, PARAMETER:: DS=0.050
171	REAL, PARAMETER:: AKS=.0000005
172	REAL, PARAMETER:: DZG=2.85
173	REAL, PARAMETER:: DI=.1000
174	REAL, PARAMETER:: AKI=0.000001075
175	REAL, PARAMETER:: DZI=2.0
176	REAL, PARAMETER:: THL=210.
177	REAL, PARAMETER:: PLQ=70000.
178	REAL, PARAMETER:: ERAD=6371200.
179	REAL, PARAMETER:: TG0=258.16
180	REAL, PARAMETER:: TGA=30.0
181
182	if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
183	if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
184
185	#define COPY_IN
186	#include <scalar_derefs.inc>
187	#ifdef DM_PARALLEL
188	# include <data_calls.inc>
189	#endif
190
191	SELECT CASE ( model_data_order )
192	CASE ( DATA_ORDER_ZXY )
193	kds = grid%sd31 ; kde = grid%ed31 ;
194	ids = grid%sd32 ; ide = grid%ed32 ;
195	jds = grid%sd33 ; jde = grid%ed33 ;
196
197	kms = grid%sm31 ; kme = grid%em31 ;
198	ims = grid%sm32 ; ime = grid%em32 ;
199	jms = grid%sm33 ; jme = grid%em33 ;
200
201	kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
202	its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
203	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
204
205	CASE ( DATA_ORDER_XYZ )
206	ids = grid%sd31 ; ide = grid%ed31 ;
207	jds = grid%sd32 ; jde = grid%ed32 ;
208	kds = grid%sd33 ; kde = grid%ed33 ;
209
210	ims = grid%sm31 ; ime = grid%em31 ;
211	jms = grid%sm32 ; jme = grid%em32 ;
212	kms = grid%sm33 ; kme = grid%em33 ;
213
214	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
215	jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
216	kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
217
218	CASE ( DATA_ORDER_XZY )
219	ids = grid%sd31 ; ide = grid%ed31 ;
220	kds = grid%sd32 ; kde = grid%ed32 ;
221	jds = grid%sd33 ; jde = grid%ed33 ;
222
223	ims = grid%sm31 ; ime = grid%em31 ;
224	kms = grid%sm32 ; kme = grid%em32 ;
225	jms = grid%sm33 ; jme = grid%em33 ;
226
227	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
228	kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
229	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
230
231	END SELECT
232
233	#ifdef DM_PARALLEL
234	CALL WRF_GET_MYPROC(MYPE)
235	CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
236	call wrf_get_nprocx(inpes)
237	call wrf_get_nprocy(jnpes)
238	!
239	allocate(itemp(inpes,jnpes),stat=istat)
240	npe=0
241	!
242	do j=1,jnpes
243	do i=1,inpes
244	itemp(i,j)=npe
245	if(npe==mype)then
246	myi=i
247	myj=j
248	endif
249	npe=npe+1
250	enddo
251	enddo
252	!
253	my_n=-1
254	if(myj+1<=jnpes)my_n=itemp(myi,myj+1)
255	!
256	my_e=-1
257	if(myi+1<=inpes)my_e=itemp(myi+1,myj)
258	!
259	my_s=-1
260	if(myj-1>=1)my_s=itemp(myi,myj-1)
261	!
262	my_w=-1
263	if(myi-1>=1)my_w=itemp(myi-1,myj)
264	!
265	my_ne=-1
266	if((myi+1<=inpes).and.(myj+1<=jnpes)) &
267	my_ne=itemp(myi+1,myj+1)
268	!
269	my_se=-1
270	if((myi+1<=inpes).and.(myj-1>=1)) &
271	my_se=itemp(myi+1,myj-1)
272	!
273	my_sw=-1
274	if((myi-1>=1).and.(myj-1>=1)) &
275	my_sw=itemp(myi-1,myj-1)
276	!
277	my_nw=-1
278	if((myi-1>=1).and.(myj+1<=jnpes)) &
279	my_nw=itemp(myi-1,myj+1)
280	!
281	! my_neb(1)=my_n
282	! my_neb(2)=my_e
283	! my_neb(3)=my_s
284	! my_neb(4)=my_w
285	! my_neb(5)=my_ne
286	! my_neb(6)=my_se
287	! my_neb(7)=my_sw
288	! my_neb(8)=my_nw
289	!
290	deallocate(itemp)
291	#endif
292
293	grid%DT=float(grid%TIME_STEP)
294
295	NNXP=min(ITE,IDE-1)
296	NNYP=min(JTE,JDE-1)
297
298	write(message,*) 'IDE, JDE: ', IDE,JDE
299	write(message,*) 'NNXP, NNYP: ', NNXP,NNYP
300	CALL wrf_message(message)
301
302	JAM=6+2*(JDE-JDS-10)
303
304	if (internal_time_loop .eq. 1) then
305	ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
306	ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
307	ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
308	ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
309	FADJ(JTS:NNYP))
310	ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
311	ALLOCATE(KHLA(JAM),KHHA(JAM))
312	ALLOCATE(KVLA(JAM),KVHA(JAM))
313	endif
314
315
316	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
317
318	write(message,*) 'cen_lat: ', config_flags%cen_lat
319	CALL wrf_debug(100,message)
320	write(message,*) 'cen_lon: ', config_flags%cen_lon
321	CALL wrf_debug(100,message)
322	write(message,*) 'dx: ', config_flags%dx
323	CALL wrf_debug(100,message)
324	write(message,*) 'dy: ', config_flags%dy
325	CALL wrf_debug(100,message)
326	write(message,*) 'config_flags%start_year: ', config_flags%start_year
327	CALL wrf_debug(100,message)
328	write(message,*) 'config_flags%start_month: ', config_flags%start_month
329	CALL wrf_debug(100,message)
330	write(message,*) 'config_flags%start_day: ', config_flags%start_day
331	CALL wrf_debug(100,message)
332	write(message,*) 'config_flags%start_hour: ', config_flags%start_hour
333	CALL wrf_debug(100,message)
334
335	write(start_date,435) config_flags%start_year, config_flags%start_month, &
336	config_flags%start_day, config_flags%start_hour
337	435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
338
339	dlmd=config_flags%dx
340	dphd=config_flags%dy
341	tph0d=config_flags%cen_lat
342	tlm0d=config_flags%cen_lon
343
344	!==========================================================================
345
346	!!
347
348	! Check to see if the boundary conditions are set
349	! properly in the namelist file.
350	! This checks for sufficiency and redundancy.
351
352	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
353
354	! Some sort of "this is the first time" initialization. Who knows.
355
356	grid%itimestep=0
357
358	! Pull in the info in the namelist to compare it to the input data.
359
360	grid%real_data_init_type = model_config_rec%real_data_init_type
361	write(message,*) 'what is flag_metgrid: ', flag_metgrid
362	CALL wrf_message(message)
363
364	IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
365
366	num_metgrid_levels = grid%num_metgrid_levels
367
368	IF (ght_gc(its,jts,num_metgrid_levels/2) .lt. ght_gc(its,jts,num_metgrid_levels/2+1)) THEN
369
370
371	write(message,*) 'normal ground up file order'
372	hyb_coor=.false.
373	CALL wrf_message(message)
374
375	ELSE
376
377	hyb_coor=.true.
378	write(message,*) 'reverse the order of coordinate'
379	CALL wrf_message(message)
380
381	CALL reverse_vert_coord(ght_gc, 2, num_metgrid_levels &
382	&, IDS,IDE,JDS,JDE,KDS,KDE &
383	&, IMS,IME,JMS,JME,KMS,KME &
384	&, ITS,ITE,JTS,JTE,KTS,KTE )
385
386	CALL reverse_vert_coord(p_gc, 2, num_metgrid_levels &
387	&, IDS,IDE,JDS,JDE,KDS,KDE &
388	&, IMS,IME,JMS,JME,KMS,KME &
389	&, ITS,ITE,JTS,JTE,KTS,KTE )
390
391	CALL reverse_vert_coord(t_gc, 2, num_metgrid_levels &
392	&, IDS,IDE,JDS,JDE,KDS,KDE &
393	&, IMS,IME,JMS,JME,KMS,KME &
394	&, ITS,ITE,JTS,JTE,KTS,KTE )
395
396	CALL reverse_vert_coord(u_gc, 2, num_metgrid_levels &
397	&, IDS,IDE,JDS,JDE,KDS,KDE &
398	&, IMS,IME,JMS,JME,KMS,KME &
399	&, ITS,ITE,JTS,JTE,KTS,KTE )
400
401	CALL reverse_vert_coord(v_gc, 2, num_metgrid_levels &
402	&, IDS,IDE,JDS,JDE,KDS,KDE &
403	&, IMS,IME,JMS,JME,KMS,KME &
404	&, ITS,ITE,JTS,JTE,KTS,KTE )
405
406	CALL reverse_vert_coord(rh_gc, 2, num_metgrid_levels &
407	&, IDS,IDE,JDS,JDE,KDS,KDE &
408	&, IMS,IME,JMS,JME,KMS,KME &
409	&, ITS,ITE,JTS,JTE,KTS,KTE )
410
411	endif
412
413
414	IF (hyb_coor) THEN
415	! limit extreme deviations from source model topography
416	! due to potential for nasty extrapolation/interpolation issues
417	!
418	write(message,*) 'min, max of ht_gc before adjust: ', minval(ht_gc), maxval(ht_gc)
419	CALL wrf_debug(100,message)
420	ICOUNT=0
421	DO J=JTS,min(JTE,JDE-1)
422	DO I=ITS,min(ITE,IDE-1)
423	IF ((ht_gc(I,J) - ght_gc(I,J,2)) .LT. -150.) THEN
424	ht_gc(I,J)=ght_gc(I,J,2)-150.
425	IF (ICOUNT .LT. 20) THEN
426	write(message,*) 'increasing NMM topo toward RUC ', I,J
427	CALL wrf_debug(100,message)
428	ICOUNT=ICOUNT+1
429	ENDIF
430	ELSEIF ((ht_gc(I,J) - ght_gc(I,J,2)) .GT. 150.) THEN
431	ht_gc(I,J)=ght_gc(I,J,2)+150.
432	IF (ICOUNT .LT. 20) THEN
433	write(message,*) 'decreasing NMM topo toward RUC ', I,J
434	CALL wrf_debug(100,message)
435	ICOUNT=ICOUNT+1
436	ENDIF
437	ENDIF
438	END DO
439	END DO
440
441	write(message,*) 'min, max of ht_gc after correct: ', minval(ht_gc), maxval(ht_gc)
442	CALL wrf_debug(100,message)
443	ENDIF
444
445	CALL boundary_smooth(ht_gc,landmask, grid, 12 , 12 &
446	&, IDS,IDE,JDS,JDE,KDS,KDE &
447	&, IMS,IME,JMS,JME,KMS,KME &
448	&, ITS,ITE,JTS,JTE,KTS,KTE )
449
450	DO j = jts, MIN(jte,jde-1)
451	DO i = its, MIN(ite,ide-1)
452	if (LANDMASK(I,J) .gt. 0.5) SM(I,J)=0.
453	if (LANDMASK(I,J) .le. 0.5) SM(I,J)=1.
454	if (tsk_gc(I,J) .gt. 0.) then
455	NMM_TSK(I,J)=tsk_gc(I,J)
456	else
457	NMM_TSK(I,J)=t_gc(I,J,1) ! stopgap measure
458	endif
459	!
460	GLAT(I,J)=hlat_gc(I,J)*DEGRAD
461	GLON(I,J)=hlon_gc(I,J)*DEGRAD
462	WEASD(I,J)=SNOW(I,J)
463	XICE(I,J)=XICE_gc(I,J)
464	ENDDO
465	ENDDO
466	! First item is to define the target vertical coordinate
467
468	num_metgrid_levels = grid%num_metgrid_levels
469	eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
470	ptsgm = model_config_rec%ptsgm
471	p_top_requested = grid%p_top_requested
472	pt=p_top_requested
473
474	if (internal_time_loop .eq. 1) then
475
476	if (eta_levels(1) .ne. 1.0) then
477	write(message,) '******************************************************************** '
478	CALL wrf_message(message)
479	write(message,) '* eta_levels appears not to be specified in the namelist'
480	CALL wrf_message(message)
481	write(message,) '* We will call compute_nmm_levels to define layer thicknesses.'
482	CALL wrf_message(message)
483	write(message,) '* These levels should be reasonable for running the model, '
484	CALL wrf_message(message)
485	write(message,) '* but may not be ideal for the simulation being made. Consider '
486	CALL wrf_message(message)
487	write(message,) '* defining your own levels by specifying eta_levels in the model '
488	CALL wrf_message(message)
489	write(message,) '* namelist. '
490	CALL wrf_message(message)
491	write(message,) '********************************************************************* '
492	CALL wrf_message(message)
493
494	CALL compute_nmm_levels(KDE,p_top_requested,eta_levels)
495
496	DO L=1,KDE
497	write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L)
498	CALL wrf_message(message)
499	ENDDO
500
501	endif
502
503	write(message,*) 'KDE-1: ', KDE-1
504	CALL wrf_debug(1,message)
505	allocate(SG1(1:KDE-1))
506	allocate(SG2(1:KDE-1))
507	allocate(DSG1(1:KDE-1))
508	allocate(DSG2(1:KDE-1))
509	allocate(SGML1(1:KDE))
510	allocate(SGML2(1:KDE))
511
512	CALL define_nmm_vertical_coord (kde-1, ptsgm, pt,pdtop, eta_levels, &
513	ETA1,DETA1,AETA1, &
514	ETA2,DETA2,AETA2, DFL, DFRLG )
515
516	DO L=KDS,KDE-1
517	DETA(L)=eta_levels(L)-eta_levels(L+1)
518	ENDDO
519	endif
520
521	if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME))
522	if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
523	if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME))
524	if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
525	if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels))
526
527	write(message,*) 'num_metgrid_levels: ', num_metgrid_levels
528	CALL wrf_message(message)
529
530	DO j = jts, MIN(jte,jde-1)
531	DO i = its, MIN(ite,ide-1)
532	FIS(I,J)=ht_gc(I,J)*g
533	!
534	! IF ( p_gc(I,J,1) .ne. 200100. .AND. (ht_gc(I,J) .eq. ght_gc(I,J,1)) .AND. ht_gc(I,J) .ne. 0) THEN
535	IF ( p_gc(I,J,1) .ne. 200100. .AND. (abs(ht_gc(I,J)-ght_gc(I,J,1)) .lt. 0.01) .AND. ht_gc(I,J) .ne. 0) THEN
536	IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN
537	write(message,*) 'ht_gc and toposoil to swap, flag_soilhgt ::: ', &
538	I,J, ht_gc(I,J),toposoil(I,J),flag_soilhgt
539	CALL wrf_debug(10,message)
540	ENDIF
541	IF ( ( flag_soilhgt.EQ. 1 ) ) THEN
542	ght_gc(I,J,1)=toposoil(I,J)
543	ENDIF
544	ENDIF
545
546	ENDDO
547	ENDDO
548
549	CALL compute_nmm_surfacep (ht_gc, ght_gc, p_gc , t_gc &
550	&, psfc_out, num_metgrid_levels &
551	&, IDS,IDE,JDS,JDE,KDS,KDE &
552	&, IMS,IME,JMS,JME,KMS,KME &
553	&, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points
554
555	CALL compute_3d_pressure (psfc_out,AETA1,AETA2 &
556	&, pdtop,pt,pd,p3d_out &
557	&, IDS,IDE,JDS,JDE,KDS,KDE &
558	&, IMS,IME,JMS,JME,KMS,KME &
559	&, ITS,ITE,JTS,JTE,KTS,KTE )
560
561	#ifdef DM_PARALLEL
562	ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
563	# include "HALO_NMM_MG2.inc"
564	#endif
565
566	#ifdef DM_PARALLEL
567	# include "HALO_NMM_MG3.inc"
568	#endif
569
570	do K=1,num_metgrid_levels
571	do J=JTS,min(JTE,JDE-1)
572	do I=ITS,min(ITE,IDE-1)
573
574	IF (K .eq. KTS) THEN
575	IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
576	PDVP(I,J)=0.5*(PD(I,J)+PD(I+1,J))
577	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J))
578	ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
579	PDVP(I,J)=0.5*(PD(I,J)+PD(I+1,J))
580	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J))
581	ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
582	PDVP(I,J)=0.5*(PD(I,J-1)+PD(I,J+1))
583	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J-1)+PSFC_OUT(I,J+1))
584	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
585	PDVP(I,J)=0.5*(PD(I,J-1)+PD(I,J+1))
586	PSFC_OUTV(I,J)=0.5*(PSFC_OUT(I,J-1)+PSFC_OUT(I,J+1))
587	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
588	PDVP(I,J)=PD(I,J)
589	PSFC_OUTV(I,J)=PSFC_OUT(I,J)
590	ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
591	PDVP(I,J)=0.25*(PD(I,J)+PD(I-1,J)+PD(I,J+1)+PD(I,J-1))
592	PSFC_OUTV(I,J)=0.25*(PSFC_OUT(I,J)+PSFC_OUT(I-1,J)+ &
593	PSFC_OUT(I,J+1)+PSFC_OUT(I,J-1))
594	ELSE ! interior odd row
595	PDVP(I,J)=0.25*(PD(I,J)+PD(I+1,J)+PD(I,J+1)+PD(I,J-1))
596	PSFC_OUTV(I,J)=0.25*(PSFC_OUT(I,J)+PSFC_OUT(I+1,J)+ &
597	PSFC_OUT(I,J+1)+PSFC_OUT(I,J-1))
598	ENDIF
599	ENDIF
600
601	IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
602	P3DV_IN(I,J,K)=0.5*(p_gc(I,J,K)+p_gc(I+1,J,K))
603	ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
604	P3DV_IN(I,J,K)=0.5*(p_gc(I,J,K)+p_gc(I+1,J,K))
605	ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
606	P3DV_IN(I,J,K)=0.5*(p_gc(I,J-1,K)+p_gc(I,J+1,K))
607	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
608	P3DV_IN(I,J,K)=0.5*(p_gc(I,J-1,K)+p_gc(I,J+1,K))
609	ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
610	P3DV_IN(I,J,K)=p_gc(I,J,K)
611	ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
612	P3DV_IN(I,J,K)=0.25*(p_gc(I,J,K)+p_gc(I-1,J,K) + &
613	p_gc(I,J+1,K)+p_gc(I,J-1,K))
614	ELSE ! interior odd row
615	P3DV_IN(I,J,K)=0.25*(p_gc(I,J,K)+p_gc(I+1,J,K) + &
616	p_gc(I,J+1,K)+p_gc(I,J-1,K))
617	ENDIF
618
619	enddo
620	enddo
621	enddo
622
623	CALL compute_3d_pressure (psfc_outv,AETA1,AETA2 &
624	&, pdtop,pt,pdvp,p3dv_out &
625	&, IDS,IDE,JDS,JDE,KDS,KDE &
626	&, IMS,IME,JMS,JME,KMS,KME &
627	&, ITS,ITE,JTS,JTE,KTS,KTE )
628
629	CALL interp_press2press_lin(p_gc, p3d_out &
630	&, t_gc, T,num_metgrid_levels &
631	&, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field
632	&, IDS,IDE,JDS,JDE,KDS,KDE &
633	&, IMS,IME,JMS,JME,KMS,KME &
634	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
635
636
637	CALL interp_press2press_lin(p3dv_in, p3dv_out &
638	&, u_gc, U,num_metgrid_levels &
639	&, .FALSE.,.TRUE.,.FALSE. &
640	&, IDS,IDE,JDS,JDE,KDS,KDE &
641	&, IMS,IME,JMS,JME,KMS,KME &
642	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
643
644	CALL interp_press2press_lin(p3dv_in, p3dv_out &
645	&, V_gc, V,num_metgrid_levels &
646	&, .FALSE.,.TRUE.,.FALSE. &
647	&, IDS,IDE,JDS,JDE,KDS,KDE &
648	&, IMS,IME,JMS,JME,KMS,KME &
649	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
650
651	IF (hyb_coor) THEN
652	CALL wind_adjust(p3dv_in,p3dv_out,U_gc,V_gc,U,V &
653	&, num_metgrid_levels,5000. &
654	&, IDS,IDE,JDS,JDE,KDS,KDE &
655	&, IMS,IME,JMS,JME,KMS,KME &
656	&, ITS,ITE,JTS,JTE,KTS,KTE )
657	ENDIF
658
659
660	ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels))
661	ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels))
662
663	CALL rh_to_mxrat (rh_gc, t_gc, p_gc, qtmp , .TRUE. , &
664	ids , ide , jds , jde , 1 , num_metgrid_levels , &
665	ims , ime , jms , jme , 1 , num_metgrid_levels , &
666	its , ite , jts , jte , 1 , num_metgrid_levels )
667
668	do K=1,num_metgrid_levels
669	do J=JTS,min(JTE,JDE-1)
670	do I=ITS,min(ITE,IDE-1)
671	QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K))
672	end do
673	end do
674	end do
675
676	CALL interp_press2press_log(p_gc, p3d_out &
677	&, QTMP2, Q,num_metgrid_levels &
678	&, .FALSE.,.TRUE. &
679	&, IDS,IDE,JDS,JDE,KDS,KDE &
680	&, IMS,IME,JMS,JME,KMS,KME &
681	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
682
683	IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP)
684	IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2)
685
686	! Get the monthly values interpolated to the current date
687	! for the traditional monthly
688	! fields of green-ness fraction and background albedo.
689
690	if (internal_time_loop .eq. 1) then
691
692	CALL monthly_interp_to_date ( greenfrac_gc , current_date , vegfra , &
693	ids , ide , jds , jde , kds , kde , &
694	ims , ime , jms , jme , kms , kme , &
695	its , ite , jts , jte , kts , kte )
696
697	CALL monthly_interp_to_date ( albedo12m_gc , current_date , albbck , &
698	ids , ide , jds , jde , kds , kde , &
699	ims , ime , jms , jme , kms , kme , &
700	its , ite , jts , jte , kts , kte )
701
702	! Get the min/max of each i,j for the monthly green-ness fraction.
703
704	CALL monthly_min_max ( greenfrac_gc , shdmin , shdmax , &
705	ids , ide , jds , jde , kds , kde , &
706	ims , ime , jms , jme , kms , kme , &
707	its , ite , jts , jte , kts , kte )
708
709	! The model expects the green-ness values in percent, not fraction.
710
711	DO j = jts, MIN(jte,jde-1)
712	DO i = its, MIN(ite,ide-1)
713	!! vegfra(i,j) = vegfra(i,j) * 100.
714	shdmax(i,j) = shdmax(i,j) * 100.
715	shdmin(i,j) = shdmin(i,j) * 100.
716	VEGFRC(I,J)=VEGFRA(I,J)
717	END DO
718	END DO
719
720	! The model expects the albedo fields as
721	! a fraction, not a percent. Set the water values to 8%.
722
723	DO j = jts, MIN(jte,jde-1)
724	DO i = its, MIN(ite,ide-1)
725	if (albbck(i,j) .lt. 5.) then
726	write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,albbck(I,J)
727	CALL wrf_debug(10,message)
728	albbck(I,J)=8.
729	endif
730	albbck(i,j) = albbck(i,j) / 100.
731	snoalb(i,j) = snoalb(i,j) / 100.
732	IF ( landmask(i,j) .LT. 0.5 ) THEN
733	albbck(i,j) = 0.08
734	snoalb(i,j) = 0.08
735	END IF
736	albase(i,j)=albbck(i,j)
737	mxsnal(i,j)=snoalb(i,j)
738	END DO
739	END DO
740
741	endif
742
743	! new deallocs
744	DEALLOCATE(p3d_out,p3dv_out,p3dv_in)
745
746	END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
747
748	if (internal_time_loop .eq. 1) then
749
750	!!! WEASD has "snow water equivalent" in mm
751
752	DO j = jts, MIN(jte,jde-1)
753	DO i = its, MIN(ite,ide-1)
754
755	IF(SM(I,J).GT.0.9) THEN
756
757	IF (XICE(I,J) .gt. 0) then
758	SI(I,J)=1.0
759	ENDIF
760
761	! SEA
762	EPSR(I,J)=.97
763	EMBCK(I,J)=.97
764	GFFC(I,J)=0.
765	ALBEDO(I,J)=.06
766	ALBASE(I,J)=.06
767	IF(SI (I,J).GT.0. ) THEN
768	! SEA-ICE
769	SM(I,J)=0.
770	SI(I,J)=0.
771	SICE(I,J)=1.
772	GFFC(I,J)=0. ! just leave zero as irrelevant
773	ALBEDO(I,J)=.60
774	ALBASE(I,J)=.60
775	ENDIF
776	ELSE
777
778	SI(I,J)=5.0*WEASD(I,J)/1000.
779	! LAND
780	EPSR(I,J)=1.0
781	EMBCK(I,J)=1.0
782	GFFC(I,J)=0.0 ! just leave zero as irrelevant
783	SICE(I,J)=0.
784	SNO(I,J)=SI(I,J)*.20
785	ENDIF
786	ENDDO
787	ENDDO
788
789	! DETERMINE ALBEDO OVER LAND
790	DO j = jts, MIN(jte,jde-1)
791	DO i = its, MIN(ite,ide-1)
792	IF(SM(I,J).LT.0.9.AND.SICE(I,J).LT.0.9) THEN
793	! SNOWFREE ALBEDO
794	IF ( (SNO(I,J) .EQ. 0.0) .OR. &
795	(ALBASE(I,J) .GE. MXSNAL(I,J) ) ) THEN
796	ALBEDO(I,J) = ALBASE(I,J)
797	ELSE
798	! MODIFY ALBEDO IF SNOWCOVER:
799	! BELOW SNOWDEPTH THRESHOLD...
800	IF (SNO(I,J) .LT. SNUP) THEN
801	RSNOW = SNO(I,J)/SNUP
802	SNOFAC = 1. - ( EXP(-SALPRSNOW) - RSNOWEXP(-SALP))
803	! ABOVE SNOWDEPTH THRESHOLD...
804	ELSE
805	SNOFAC = 1.0
806	ENDIF
807	! CALCULATE ALBEDO ACCOUNTING FOR SNOWDEPTH AND VGFRCK
808	ALBEDO(I,J) = ALBASE(I,J) &
809	+ (1.0-VEGFRA(I,J))SNOFAC(MXSNAL(I,J)-ALBASE(I,J))
810	ENDIF
811	END IF
812	SI(I,J)=5.0*WEASD(I,J)
813	SNO(I,J)=WEASD(I,J)
814
815	!! convert VEGFRA
816	VEGFRA(I,J)=VEGFRA(I,J)*100.
817	!
818	ENDDO
819	ENDDO
820
821	#ifdef DM_PARALLEL
822
823	ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
824
825	CALL WRF_PATCH_TO_GLOBAL_REAL( SICE(IMS,JMS) &
826	&, SICE_G,grid%DOMDESC &
827	&, 'z','xy' &
828	&, IDS,IDE-1,JDS,JDE-1,1,1 &
829	&, IMS,IME,JMS,JME,1,1 &
830	&, ITS,ITE,JTS,JTE,1,1 )
831
832	CALL WRF_PATCH_TO_GLOBAL_REAL( SM(IMS,JMS) &
833	&, SM_G,grid%DOMDESC &
834	&, 'z','xy' &
835	&, IDS,IDE-1,JDS,JDE-1,1,1 &
836	&, IMS,IME,JMS,JME,1,1 &
837	&, ITS,ITE,JTS,JTE,1,1 )
838
839
840	IF (WRF_DM_ON_MONITOR()) THEN
841
842	637 format(40(f3.0,1x))
843
844	allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
845	DO j = JDS, JDE-1
846	IHE_G(J)=MOD(J+1,2)
847	IHW_G(J)=IHE_G(J)-1
848	ENDDO
849
850	DO ITER=1,10
851	DO j = jds+1, (jde-1)-1
852	DO i = ids+1, (ide-1)-1
853
854	! any sea ice around point in question?
855
856	IF (SM_G(I,J) .ge. 0.9) THEN
857	SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
858	SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
859	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
860
861	IF ((SICE_G(I+IHE_G(J),J+1).lt.0.1 .and. SM_G(I+IHE_G(J),J+1).lt.0.1) .OR. &
862	(SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. &
863	(SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. &
864	(SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN
865
866	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
867
868	write(message,*) 'making seaice (1): ', I,J
869	CALL wrf_debug(100,message)
870	SICE_G(I,J)=1.0
871	SM_G(I,J)=0.
872
873	ENDIF
874
875	ELSEIF (SEAICESUM .ge. 3) THEN
876
877	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
878
879	write(message,*) 'making seaice (2): ', I,J
880	CALL wrf_debug(100,message)
881	SICE_G(I,J)=1.0
882	SM_G(I,J)=0.
883	ENDIF
884
885	ENDIF
886
887	ENDDO
888	ENDDO
889	ENDDO
890
891	ENDIF
892
893	CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, SICE &
894	&, grid%DOMDESC &
895	&, 'z','xy' &
896	&, IDS,IDE-1,JDS,JDE-1,1,1 &
897	&, IMS,IME,JMS,JME,1,1 &
898	&, ITS,ITE,JTS,JTE,1,1 )
899
900	CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,SM &
901	&, grid%DOMDESC &
902	&, 'z','xy' &
903	&, IDS,IDE-1,JDS,JDE-1,1,1 &
904	&, IMS,IME,JMS,JME,1,1 &
905	&, ITS,ITE,JTS,JTE,1,1 )
906
907	IF (WRF_DM_ON_MONITOR()) THEN
908
909	DEALLOCATE(SM_G,SICE_G)
910	DEALLOCATE(IHE_G,IHW_G)
911
912	ENDIF
913
914	! write(message,*) 'revised sea ice on patch'
915	! CALL wrf_debug(100,message)
916	! DO J=JTE,JTS,-(((JTE-JTS)/25)+1)
917	! write(message,637) (SICE(I,J),I=ITS,ITE,ITE/20)
918	! CALL wrf_debug(100,message)
919	! END DO
920
921	#else
922	! serial sea ice reprocessing
923
924	DO j = jts, MIN(jte,jde-1)
925	IHE(J)=MOD(J+1,2)
926	IHW(J)=IHE(J)-1
927	ENDDO
928
929	DO ITER=1,10
930	DO j = jts+1, MIN(jte,jde-1)-1
931	DO i = its+1, MIN(ite,ide-1)-1
932
933	! any sea ice around point in question?
934
935	IF (SM(I,J) .gt. 0.9) THEN
936	SEAICESUM=SICE(I+IHE(J),J+1)+SICE(I+IHW(J),J+1)+ &
937	SICE(I+IHE(J),J-1)+SICE(I+IHW(J),J-1)
938	IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
939	IF ((SICE(I+IHE(J),J+1).lt.0.1 .and. SM(I+IHE(J),J+1).lt.0.1) .OR. &
940	(SICE(I+IHW(J),J+1).lt.0.1 .and. SM(I+IHW(J),J+1).lt.0.1) .OR. &
941	(SICE(I+IHE(J),J-1).lt.0.1 .and. SM(I+IHE(J),J-1).lt.0.1) .OR. &
942	(SICE(I+IHW(J),J-1).lt.0.1 .and. SM(I+IHW(J),J-1).lt.0.1)) THEN
943
944	! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
945	SICE(I,J)=1.0
946	SM(I,J)=0.
947	ENDIF
948	ELSEIF (SEAICESUM .ge. 3) THEN
949	! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
950	SICE(I,J)=1.0
951	SM(I,J)=0.
952	ENDIF
953	ENDIF
954
955	ENDDO
956	ENDDO
957	ENDDO
958
959	#endif
960
961	! this block meant to guarantee land/sea agreement between SM and landmask
962
963	DO j = jts, MIN(jte,jde-1)
964	DO i = its, MIN(ite,ide-1)
965
966	IF (SM(I,J) .gt. 0.5) THEN
967	landmask(I,J)=0.0
968	ELSEIF (SM(I,J) .lt. 0.5 .and. SICE(I,J) .gt. 0.9) then
969	landmask(I,J)=0.0
970	ELSEIF (SM(I,J) .lt. 0.5 .and. SICE(I,J) .lt. 0.1) then
971	landmask(I,J)=1.0
972	ELSE
973	write(message,*) 'missed point in landmask definition ' , I,J
974	CALL wrf_message(message)
975	landmask(I,J)=0.0
976	ENDIF
977	!
978	IF (SICE(I,J) .gt. 0.5 .and. NMM_TSK(I,J) .lt. 0.1 .and. SST(I,J) .gt. 0.) THEN
979	write(message,*) 'set NMM_TSK to: ', SST(I,J)
980	CALL wrf_message(message)
981	NMM_TSK(I,J)=SST(I,J)
982	SST(I,J)=0.
983	endif
984
985	ENDDO
986	ENDDO
987
988	! For sf_surface_physics = 1, we want to use close to a 10 cm value
989	! for the bottom level of the soil temps.
990
991	IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
992	( flag_st000010 .EQ. 1 ) ) THEN
993	DO j = jts , MIN(jde-1,jte)
994	DO i = its , MIN(ide-1,ite)
995	soiltb(i,j) = st000010(i,j)
996	END DO
997	END DO
998	END IF
999
1000	! Adjust the various soil temperature values depending on the difference in
1001	! in elevation between the current model's elevation and the incoming data's
1002	! orography.
1003
1004	IF ( ( flag_toposoil .EQ. 1 ) ) THEN
1005
1006	ALLOCATE(HT(ims:ime,jms:jme))
1007
1008	DO J=jms,jme
1009	DO I=ims,ime
1010	HT(I,J)=FIS(I,J)/9.81
1011	END DO
1012	END DO
1013
1014	! if (maxval(toposoil) .gt. 100.) then
1015	!
1016	! Being avoided. Something to revisit eventually.
1017	!
1018	!1219 might be simply a matter of including TOPOSOIL
1019	!
1020	! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
1021	! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
1022	!
1023	! CALL adjust_soil_temp_new ( soiltb , 2 , &
1024	! nmm_tsk , ht , toposoil , landmask, flag_toposoil , &
1025	! st000010 , st010040 , st040100 , st100200 , st010200 , &
1026	! flag_st000010 , flag_st010040 , flag_st040100 , &
1027	! flag_st100200 , flag_st010200 , &
1028	! soilt000 , soilt005 , soilt020 , soilt040 , &
1029	! soilt160 , soilt300 , &
1030	! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
1031	! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
1032	! ids , ide , jds , jde , kds , kde , &
1033	! ims , ime , jms , jme , kms , kme , &
1034	! its , ite , jts , jte , kts , kte )
1035	! endif
1036
1037	END IF
1038
1039	! Process the LSM data.
1040
1041	! surface_input_source=1 => use data from static file
1042	! (fractional category as input)
1043	! surface_input_source=2 => use data from grib file
1044	! (dominant category as input)
1045
1046	IF ( config_flags%surface_input_source .EQ. 1 ) THEN
1047	vegcat (its,jts) = 0
1048	soilcat(its,jts) = 0
1049	END IF
1050
1051	! Generate the vegetation and soil category information
1052	! from the fractional input
1053	! data, or use the existing dominant category fields if they exist.
1054
1055	IF ((soilcat(its,jts) .LT. 0.5) .AND. (vegcat(its,jts) .LT. 0.5)) THEN
1056
1057	num_veg_cat = SIZE ( landusef_gc , DIM=3 )
1058	num_soil_top_cat = SIZE ( soilctop_gc , DIM=3 )
1059	num_soil_bot_cat = SIZE ( soilcbot_gc , DIM=3 )
1060
1061	do J=JMS,JME
1062	do K=1,num_veg_cat
1063	do I=IMS,IME
1064	landusef(I,K,J)=landusef_gc(I,J,K)
1065	enddo
1066	enddo
1067	enddo
1068
1069	do J=JMS,JME
1070	do K=1,num_soil_top_cat
1071	do I=IMS,IME
1072	soilctop(I,K,J)=soilctop_gc(I,J,K)
1073	enddo
1074	enddo
1075	enddo
1076
1077	do J=JMS,JME
1078	do K=1,num_soil_bot_cat
1079	do I=IMS,IME
1080	soilcbot(I,K,J)=soilcbot_gc(I,J,K)
1081	enddo
1082	enddo
1083	enddo
1084
1085	! sm (1=water, 0=land)
1086	! landmask(0=water, 1=land)
1087
1088
1089	write(message,*) 'landmask into process_percent_cat_new'
1090
1091	CALL wrf_debug(1,message)
1092	do J=JTE,JTS,-(((JTE-JTS)/20)+1)
1093	write(message,641) (landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1)
1094	CALL wrf_debug(1,message)
1095	enddo
1096	641 format(25(f3.0,1x))
1097
1098	CALL process_percent_cat_new ( landmask , &
1099	landusef , soilctop , soilcbot , &
1100	isltyp , ivgtyp , &
1101	num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1102	ids , ide , jds , jde , kds , kde , &
1103	ims , ime , jms , jme , kms , kme , &
1104	its , ite , jts , jte , kts , kte , &
1105	model_config_rec%iswater(grid%id) )
1106
1107	DO j = jts , MIN(jde-1,jte)
1108	DO i = its , MIN(ide-1,ite)
1109	vegcat(i,j) = ivgtyp(i,j)
1110	soilcat(i,j) = isltyp(i,j)
1111	END DO
1112	END DO
1113
1114	ELSE
1115
1116	! Do we have dominant soil and veg data from the input already?
1117
1118	IF ( soilcat(its,jts) .GT. 0.5 ) THEN
1119	DO j = jts, MIN(jde-1,jte)
1120	DO i = its, MIN(ide-1,ite)
1121	isltyp(i,j) = NINT( soilcat(i,j) )
1122	END DO
1123	END DO
1124	END IF
1125	IF ( vegcat(its,jts) .GT. 0.5 ) THEN
1126	DO j = jts, MIN(jde-1,jte)
1127	DO i = its, MIN(ide-1,ite)
1128	ivgtyp(i,j) = NINT( vegcat(i,j) )
1129	END DO
1130	END DO
1131	END IF
1132
1133	ENDIF
1134
1135	DO j = jts, MIN(jde-1,jte)
1136	DO i = its, MIN(ide-1,ite)
1137
1138	IF (SICE(I,J) .lt. 0.1) THEN
1139	IF (landmask(I,J) .gt. 0.5 .and. sm(I,J) .gt. 0.5) THEN
1140	write(message,*) 'land mask and SM both > 0.5: ', &
1141	I,J,landmask(I,J),sm(I,J)
1142	CALL wrf_message(message)
1143	SM(I,J)=0.
1144	ELSEIF (landmask(I,J) .lt. 0.5 .and. sm(I,J) .lt. 0.5) THEN
1145	write(message,*) 'land mask and SM both < 0.5: ', &
1146	I,J, landmask(I,J),sm(I,J)
1147	CALL wrf_message(message)
1148	SM(I,J)=1.
1149	ENDIF
1150	ELSE
1151	IF (landmask(I,J) .gt. 0.5 .and. SM(I,J)+SICE(I,J) .gt. 0.9) then
1152	write(message,*) 'landmask says LAND, SM/SICE say SEAICE: ', I,J
1153	ENDIF
1154	ENDIF
1155
1156	ENDDO
1157	ENDDO
1158
1159	DO j = jts, MIN(jde-1,jte)
1160	DO i = its, MIN(ide-1,ite)
1161
1162	if (SICE(I,J) .gt. 0.9) then
1163	ISLTYP(I,J)=16
1164	IVGTYP(I,J)=24
1165	endif
1166
1167	ENDDO
1168	ENDDO
1169
1170	DO j = jts, MIN(jde-1,jte)
1171	DO i = its, MIN(ide-1,ite)
1172
1173	if (SM(I,J) .lt. 0.5) then
1174	SST(I,J)=0.
1175	endif
1176
1177	if (SM(I,J) .gt. 0.5) then
1178	if (SST(I,J) .lt. 0.1) then
1179	SST(I,J)=NMM_TSK(I,J)
1180	endif
1181	NMM_TSK(I,J)=0.
1182	endif
1183
1184	IF ( (NMM_TSK(I,J)+SST(I,J)) .lt. 200. .or. &
1185	(NMM_TSK(I,J)+SST(I,J)) .gt. 350. ) THEN
1186	write(message,*) 'TSK, SST trouble at : ', I,J
1187	CALL wrf_message(message)
1188	write(message,*) 'SM, NMM_TSK,SST ', SM(I,J),NMM_TSK(I,J),SST(I,J)
1189	CALL wrf_message(message)
1190	ENDIF
1191
1192	ENDDO
1193	ENDDO
1194
1195	write(message,*) 'SM'
1196	CALL wrf_message(message)
1197
1198	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1199	write(message,635) (sm(i,J),I=its,ite,((ite-its)/10)+1)
1200	CALL wrf_message(message)
1201	END DO
1202
1203	write(message,*) 'SST/NMM_TSK'
1204	CALL wrf_debug(10,message)
1205	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1206	write(message,635) (SST(I,J)+NMM_TSK(I,J),I=ITS,min(ide-1,ite),((ite-its)/10)+1)
1207	CALL wrf_debug(10,message)
1208	END DO
1209
1210	635 format(20(f5.1,1x))
1211
1212	DO j = jts, MIN(jde-1,jte)
1213	DO i = its, MIN(ide-1,ite)
1214	IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
1215	soiltb(i,j) = sst(i,j)
1216	ELSE IF ( landmask(i,j) .GT. 0.5 ) THEN
1217	soiltb(i,j) = nmm_tsk(i,j)
1218	END IF
1219	END DO
1220	END DO
1221
1222	! END IF
1223
1224	! Land use categories, dominant soil and vegetation types (if available).
1225
1226	! allocate(lu_index(ims:ime,jms:jme))
1227
1228	DO j = jts, MIN(jde-1,jte)
1229	DO i = its, MIN(ide-1,ite)
1230	lu_index(i,j) = ivgtyp(i,j)
1231	END DO
1232	END DO
1233
1234	if (flag_sst .eq. 1) log_flag_sst=.true.
1235	if (flag_sst .eq. 0) log_flag_sst=.false.
1236
1237	write(message,*) 'st_input dimensions: ', size(st_input,dim=1), &
1238	size(st_input,dim=2),size(st_input,dim=3)
1239	CALL wrf_debug(100,message)
1240
1241	! write(message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
1242	! CALL wrf_message(message)
1243	! write(message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
1244	! CALL wrf_message(message)
1245	! write(message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
1246	! CALL wrf_message(message)
1247	! write(message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
1248	! CALL wrf_message(message)
1249
1250	! =============================================================
1251
1252	IF (.NOT. ALLOCATED(TG_ALT))ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme))
1253
1254	TPH0=TPH0D*DTR
1255	WBD=-(((ide-1)-1)*DLMD)
1256	WB= WBD*DTR
1257	SBD=-(((jde-1)/2)*DPHD)
1258	SB= SBD*DTR
1259	DLM=DLMD*DTR
1260	DPH=DPHD*DTR
1261	TDLM=DLM+DLM
1262	TDPH=DPH+DPH
1263	WBI=WB+TDLM
1264	SBI=SB+TDPH
1265	EBI=WB+(ide-2)*TDLM
1266	ANBI=SB+(jde-2)*DPH
1267	STPH0=SIN(TPH0)
1268	CTPH0=COS(TPH0)
1269	TSPH=3600./GRID%DT
1270	DO J=JTS,min(JTE,JDE-1)
1271	TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid
1272	TPH=SB+float(J-1)*DPH
1273	STPH=SIN(TPH)
1274	CTPH=COS(TPH)
1275	DO I=ITS,MIN(ITE,IDE-1)
1276
1277	if (I .eq. ITS) THEN
1278	TLM=TLM+TDLM*ITS
1279	else
1280	TLM=TLM+TDLM
1281	endif
1282
1283	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
1284	FP=TWOM*(TERM1)
1285	F(I,J)=0.5GRID%DTFP
1286	ENDDO
1287	ENDDO
1288	DO J=JTS,min(JTE,JDE-1)
1289	TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid
1290	TPH=SB+float(J-1)*DPH
1291	STPH=SIN(TPH)
1292	CTPH=COS(TPH)
1293	DO I=ITS,MIN(ITE,IDE-1)
1294
1295	if (I .eq. ITS) THEN
1296	TLM=TLM+TDLM*ITS
1297	else
1298	TLM=TLM+TDLM
1299	endif
1300
1301	TERM1=(STPH0CTPHCOS(TLM)+CTPH0*STPH)
1302	TERM1=MIN(TERM1,1.0D0)
1303	TERM1=MAX(TERM1,-1.0D0)
1304	APH=ASIN(TERM1)
1305	TG_ALT(I,J)=TG0+TGA*COS(APH)-FIS(I,J)/3333.
1306	ENDDO
1307	ENDDO
1308
1309	DO j = jts, MIN(jde-1,jte)
1310	DO i = its, MIN(ide-1,ite)
1311	! IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1312	! SICE(I,J) .eq. 0. ) THEN
1313	! TG(i,j) = sst(i,j)
1314	! ELSEIF (SICE(I,J) .eq. 1) THEN
1315	! TG(i,j) = 271.16
1316	! END IF
1317
1318	if (TG(I,J) .lt. 200.) then ! only use default TG_ALT definition if
1319	! not getting TGROUND from SI
1320	TG(I,J)=TG_ALT(I,J)
1321	endif
1322
1323	if (TG(I,J) .lt. 200. .or. TG(I,J) .gt. 320.) then
1324	write(message,*) 'problematic TG point at : ', I,J
1325	CALL wrf_message( message )
1326	endif
1327
1328	adum2d(i,j)=nmm_tsk(I,J)+sst(I,J)
1329
1330	END DO
1331	END DO
1332
1333	DEALLOCATE(TG_ALT)
1334
1335	write(message,*) 'call process_soil_real with num_st_levels_input: ', num_st_levels_input
1336	CALL wrf_message( message )
1337
1338	! =============================================================
1339
1340	CALL process_soil_real ( adum2d, TG , &
1341	landmask, sst, &
1342	st_input, sm_input, sw_input, &
1343	st_levels_input , sm_levels_input , &
1344	sw_levels_input , &
1345	sldpth , dzsoil , stc , smc , sh2o, &
1346	flag_sst , flag_soilt000, flag_soilm000, &
1347	ids , ide , jds , jde , kds , kde , &
1348	ims , ime , jms , jme , kms , kme , &
1349	its , ite , jts , jte , kts , kte , &
1350	model_config_rec%sf_surface_physics(grid%id) , &
1351	model_config_rec%num_soil_layers , &
1352	model_config_rec%real_data_init_type , &
1353	num_st_levels_input , num_sm_levels_input , &
1354	num_sw_levels_input , &
1355	num_st_levels_alloc , num_sm_levels_alloc , &
1356	num_sw_levels_alloc )
1357
1358	! =============================================================
1359
1360	! Minimum soil values, residual, from RUC LSM scheme.
1361	! For input from Noah and using
1362	! RUC LSM scheme, this must be subtracted from the input
1363	! total soil moisture. For input RUC data and using the Noah LSM scheme,
1364	! this value must be added to the soil moisture_input.
1365
1366	lqmi(1:num_soil_top_cat) = &
1367	(/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1368	0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1369	0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /)
1370
1371	! At the initial time we care about values of soil moisture and temperature,
1372	! other times are ignored by the model, so we ignore them, too.
1373
1374	account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1375
1376	CASE ( LSMSCHEME , NMMLSMSCHEME)
1377	iicount = 0
1378	IF ( FLAG_SM000010 .EQ. 1 ) THEN
1379	DO j = jts, MIN(jde-1,jte)
1380	DO i = its, MIN(ide-1,ite)
1381	IF ((landmask(i,j).gt.0.5) .and. (stc(i,1,j) .gt. 200) .and. &
1382	(stc(i,1,j) .lt. 400) .and. (smc(i,1,j) .lt. 0.005)) then
1383	write(message,*) 'Noah > Noah: bad soil moisture at i,j = ',i,j,smc(i,:,j)
1384	CALL wrf_message(message)
1385	iicount = iicount + 1
1386	smc(i,:,j) = 0.005
1387	END IF
1388	END DO
1389	END DO
1390	IF ( iicount .GT. 0 ) THEN
1391	write(message,*) 'Noah -> Noah: total number of small soil moisture locations= ',&
1392	iicount
1393	CALL wrf_message(message)
1394	END IF
1395	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1396	DO j = jts, MIN(jde-1,jte)
1397	DO i = its, MIN(ide-1,ite)
1398	smc(i,:,j) = smc(i,:,j) + lqmi(isltyp(i,j))
1399	END DO
1400	END DO
1401	DO j = jts, MIN(jde-1,jte)
1402	DO i = its, MIN(ide-1,ite)
1403	IF ((landmask(i,j).gt.0.5) .and. (stc(i,1,j) .gt. 200) .and. &
1404	(stc(i,1,j) .lt. 400) .and. (smc(i,1,j) .lt. 0.004)) then
1405	write(message,*) 'RUC -> Noah: bad soil moisture at i,j = ' &
1406	,i,j,smc(i,:,j)
1407	CALL wrf_message(message)
1408	iicount = iicount + 1
1409	smc(i,:,j) = 0.004
1410	END IF
1411	END DO
1412	END DO
1413	IF ( iicount .GT. 0 ) THEN
1414	write(message,*) 'RUC -> Noah: total number of small soil moisture locations = ',&
1415	iicount
1416	CALL wrf_message(message)
1417	END IF
1418	END IF
1419	CASE ( RUCLSMSCHEME )
1420	iicount = 0
1421	IF ( FLAG_SM000010 .EQ. 1 ) THEN
1422	DO j = jts, MIN(jde-1,jte)
1423	DO i = its, MIN(ide-1,ite)
1424	smc(i,:,j) = MAX ( smc(i,:,j) - lqmi(isltyp(i,j)) , 0. )
1425	END DO
1426	END DO
1427	ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1428	! no op
1429	END IF
1430
1431	END SELECT account_for_zero_soil_moisture
1432
1433	!!! zero out NMM_TSK at water points again
1434
1435	DO j = jts, MIN(jde-1,jte)
1436	DO i = its, MIN(ide-1,ite)
1437	if (SM(I,J) .gt. 0.5) then
1438	NMM_TSK(I,J)=0.
1439	endif
1440	END DO
1441	END DO
1442
1443	!! check on STC
1444
1445	DO j = jts, MIN(jde-1,jte)
1446	DO i = its, MIN(ide-1,ite)
1447
1448	IF (SICE(I,J) .gt. 0.9) then
1449	DO L = 1, grid%num_soil_layers
1450	STC(I,L,J)=271.16 ! TG value used by Eta/NMM
1451	END DO
1452	END IF
1453
1454	IF (SM(I,J) .gt. 0.9) then
1455	DO L = 1, grid%num_soil_layers
1456	STC(I,L,J)=273.16 ! TG value used by Eta/NMM
1457	END DO
1458	END IF
1459
1460	END DO
1461	END DO
1462
1463	DO j = jts, MIN(jde-1,jte)
1464	DO i = its, MIN(ide-1,ite)
1465
1466	if (SM(I,J) .lt. 0.1 .and. STC(I,1,J) .lt. 0.1) THEN
1467	write(message,*) 'troublesome SM,STC,SMC value: ', I,J,SM(I,J), stc(I,1,J),smc(I,1,J)
1468	CALL wrf_message(message)
1469	do JJ=J-1,J+1
1470	do L=1, grid%num_soil_layers
1471	do II=I-1,I+1
1472
1473	if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. &
1474	JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then
1475
1476	STC(I,L,J)=amax1(STC(I,L,J),STC(II,L,JJ))
1477	cur_smc=SMC(I,L,J)
1478
1479	if ( SMC(II,L,JJ) .gt. 0.005 .and. SMC(II,L,JJ) .lt. 1.0) then
1480	aposs_smc=SMC(II,L,JJ)
1481
1482	if ( cur_smc .eq. 0 ) then
1483	cur_smc=aposs_smc
1484	SMC(I,L,J)=cur_smc
1485	else
1486	cur_smc=amin1(cur_smc,aposs_smc)
1487	cur_smc=amin1(cur_smc,aposs_smc)
1488	SMC(I,L,J)=cur_smc
1489	endif
1490	endif
1491
1492	endif ! bounds check
1493
1494	enddo
1495	enddo
1496	enddo
1497	write(message,*) 'STC, SMC(1) now: ', stc(I,1,J),smc(I,1,J)
1498	CALL wrf_message(message)
1499	endif
1500
1501	if (STC(I,1,J) .lt. 0.1) then
1502	write(message,*) 'QUITTING DUE TO STILL troublesome STC value: ', I,J, stc(I,1,J),smc(I,1,J)
1503	call wrf_error_fatal(message)
1504	endif
1505
1506	ENDDO
1507	ENDDO
1508
1509	!hardwire soil stuff for time being
1510
1511	! RTDPTH=0.
1512	! RTDPTH(1)=0.1
1513	! RTDPTH(2)=0.3
1514	! RTDPTH(3)=0.6
1515
1516	! SLDPTH=0.
1517	! SLDPTH(1)=0.1
1518	! SLDPTH(2)=0.3
1519	! SLDPTH(3)=0.6
1520	! SLDPTH(4)=1.0
1521
1522	!!! main body of nmm_specific starts here
1523	!
1524	do J=jts,min(jte,jde-1)
1525	do I=its,min(ite,ide-1)
1526	RES(I,J)=1.
1527	enddo
1528	enddo
1529
1530	!! HBM2
1531
1532	HBM2=0.
1533
1534	do J=jts,min(jte,jde-1)
1535	do I=its,min(ite,ide-1)
1536
1537	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1538	(I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN
1539	HBM2(I,J)=1.
1540	ENDIF
1541	enddo
1542	enddo
1543
1544	!! HBM3
1545	HBM3=0.
1546
1547	!! LOOP OVER LOCAL DIMENSIONS
1548
1549	do J=jts,min(jte,jde-1)
1550	IHWG(J)=mod(J+1,2)-1
1551	IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN
1552	IHL=(ids+1)-IHWG(J)
1553	IHH=(ide-1)-2
1554	do I=its,min(ite,ide-1)
1555	IF (I .ge. IHL .and. I .le. IHH) HBM3(I,J)=1.
1556	enddo
1557	ENDIF
1558	enddo
1559
1560	!! VBM2
1561
1562	VBM2=0.
1563
1564	do J=jts,min(jte,jde-1)
1565	do I=its,min(ite,ide-1)
1566
1567	IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1568	(I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN
1569
1570	VBM2(I,J)=1.
1571
1572	ENDIF
1573
1574	enddo
1575	enddo
1576
1577	!! VBM3
1578
1579	VBM3=0.
1580
1581	do J=jts,min(jte,jde-1)
1582	do I=its,min(ite,ide-1)
1583
1584	IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. &
1585	(I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN
1586	VBM3(I,J)=1.
1587	ENDIF
1588
1589	enddo
1590	enddo
1591
1592	DTAD=1.0
1593	! IDTCF=DTCF, IDTCF=4
1594	DTCF=4.0 ! used?
1595
1596	DY_NMM=ERAD*DPH
1597	CPGFV=-GRID%DT/(48.*DY_NMM)
1598	EN= GRID%DT/( 4.DY_NMM)DTAD
1599	ENT=GRID%DT/(16.DY_NMM)DTAD
1600
1601	DO J=jts,nnyp
1602	KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
1603	KVL2(J)=(IDE-1)*(J-1)-J/2+2
1604	KHH2(J)=(IDE-1)*J-J/2-1
1605	KVH2(J)=(IDE-1)*J-(J+1)/2-1
1606	ENDDO
1607
1608	TPH=SB-DPH
1609
1610	DO J=jts,min(jte,jde-1)
1611	TPH=SB+float(J-1)*DPH
1612	DXP=ERADDLMCOS(TPH)
1613	DXJ(J)=DXP
1614	WPDARJ(J)=-W_NMM * &
1615	((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM2)/ &
1616	(GRID%DT32.DXP*DY_NMM)
1617
1618	CPGFUJ(J)=-GRID%DT/(48.*DXP)
1619	CURVJ(J)=.5GRID%DTTAN(TPH)/ERAD
1620	FCPJ(J)=GRID%DT/(CP192.DXP*DY_NMM)
1621	FDIVJ(J)=1./(12.DXPDY_NMM)
1622	! EMJ(J)= GRID%DT/( 4.DXP)DTAD
1623	! EMTJ(J)=GRID%DT/(16.DXP)DTAD
1624	FADJ(J)=-GRID%DT/(48.DXPDY_NMM)*DTAD
1625	ACDT=GRID%DTSQRT((ERADDLMAMIN1(COS(ANBI),COS(SBI)))2+DY_NMM*2)
1626	CDDAMP=CODAMP*ACDT
1627	HDACJ(J)=COACACDT/(4.DXP*DY_NMM)
1628	DDMPUJ(J)=CDDAMP/DXP
1629	DDMPVJ(J)=CDDAMP/DY_NMM
1630	ENDDO
1631
1632	DO J=JTS,min(JTE,JDE-1)
1633	TLM=WB-TDLM+MOD(J,2)*DLM
1634	TPH=SB+float(J-1)*DPH
1635	STPH=SIN(TPH)
1636	CTPH=COS(TPH)
1637	DO I=ITS,MIN(ITE,IDE-1)
1638
1639	if (I .eq. ITS) THEN
1640	TLM=TLM+TDLM*ITS
1641	else
1642	TLM=TLM+TDLM
1643	endif
1644
1645	FP=TWOM(CTPH0STPH+STPH0CTPHCOS(TLM))
1646	F(I,J)=0.5GRID%DTFP
1647
1648	ENDDO
1649	ENDDO
1650
1651	! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
1652
1653	EF4T=.5*GRID%DT/CP
1654	F4Q = -GRID%DT*DTAD
1655	F4D =-.5GRID%DTDTAD
1656
1657	DO L=KDS,KDE-1
1658	RDETA(L)=1./DETA(L)
1659	F4Q2(L)=-.25GRID%DTDTAD/DETA(L)
1660	ENDDO
1661
1662	DO J=JTS,min(JTE,JDE-1)
1663	DO I=ITS,min(ITE,IDE-1)
1664	DX_NMM(I,J)=DXJ(J)
1665	WPDAR(I,J)=WPDARJ(J)*HBM2(I,J)
1666	CPGFU(I,J)=CPGFUJ(J)*VBM2(I,J)
1667	CURV(I,J)=CURVJ(J)*VBM2(I,J)
1668	FCP(I,J)=FCPJ(J)*HBM2(I,J)
1669	FDIV(I,J)=FDIVJ(J)*HBM2(I,J)
1670	FAD(I,J)=FADJ(J)
1671	HDACV(I,J)=HDACJ(J)*VBM2(I,J)
1672	HDAC(I,J)=HDACJ(J)1.25HBM2(I,J)
1673	ENDDO
1674	ENDDO
1675
1676	DO J=JTS, MIN(JDE-1,JTE)
1677
1678	IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
1679
1680	KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
1681	DO I=ITS,MIN(IDE-1,ITE)
1682	IF (I .ge. 2 .and. I .le. KHH) THEN
1683	HDAC(I,J)=HDAC(I,J)* DFC
1684	ENDIF
1685	ENDDO
1686
1687	ELSE
1688
1689	KHH=2+MOD(J,2)
1690	DO I=ITS,MIN(IDE-1,ITE)
1691	IF (I .ge. 2 .and. I .le. KHH) THEN
1692	HDAC(I,J)=HDAC(I,J)* DFC
1693	ENDIF
1694	ENDDO
1695
1696	KHH=(IDE-1)-2+MOD(J,2)
1697
1698	DO I=ITS,MIN(IDE-1,ITE)
1699	IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
1700	HDAC(I,J)=HDAC(I,J)* DFC
1701	ENDIF
1702	ENDDO
1703	ENDIF
1704	ENDDO
1705
1706	DO J=JTS,min(JTE,JDE-1)
1707	DO I=ITS,min(ITE,IDE-1)
1708	DDMPU(I,J)=DDMPUJ(J)*VBM2(I,J)
1709	DDMPV(I,J)=DDMPVJ(J)*VBM2(I,J)
1710	HDACV(I,J)=HDACV(I,J)*VBM2(I,J)
1711	ENDDO
1712	ENDDO
1713	! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
1714
1715	DO J=JTS,MIN(JDE-1,JTE)
1716	IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
1717	KVH=(IDE-1)-1-MOD(J,2)
1718	DO I=ITS,min(IDE-1,ITE)
1719	IF (I .ge. 2 .and. I .le. KVH) THEN
1720	DDMPU(I,J)=DDMPU(I,J)*DDFC
1721	DDMPV(I,J)=DDMPV(I,J)*DDFC
1722	HDACV(I,J)=HDACV(I,J)* DFC
1723	ENDIF
1724	ENDDO
1725	ELSE
1726	KVH=3-MOD(J,2)
1727	DO I=ITS,min(IDE-1,ITE)
1728	IF (I .ge. 2 .and. I .le. KVH) THEN
1729	DDMPU(I,J)=DDMPU(I,J)*DDFC
1730	DDMPV(I,J)=DDMPV(I,J)*DDFC
1731	HDACV(I,J)=HDACV(I,J)* DFC
1732	ENDIF
1733	ENDDO
1734	KVH=(IDE-1)-1-MOD(J,2)
1735	DO I=ITS,min(IDE-1,ITE)
1736	IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
1737	DDMPU(I,J)=DDMPU(I,J)*DDFC
1738	DDMPV(I,J)=DDMPV(I,J)*DDFC
1739	HDACV(I,J)=HDACV(I,J)* DFC
1740	ENDIF
1741	ENDDO
1742	ENDIF
1743	ENDDO
1744
1745	write(message,*) 'STC(1)'
1746	CALL wrf_message(message)
1747	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1748	write(message,635) (stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
1749	CALL wrf_message(message)
1750	ENDDO
1751
1752	write(message,*) 'SMC(1)'
1753	CALL wrf_message(message)
1754	DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1755	write(message,635) (smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
1756	CALL wrf_message(message)
1757	ENDDO
1758
1759	DO j = jts, MIN(jde-1,jte)
1760	DO i= ITS, MIN(IDE-1,ITE)
1761
1762	if (SM(I,J) .lt. 0.1 .and. SMC(I,1,J) .gt. 0.5 .and. SICE(I,J) .lt. 0.1) then
1763	write(message,*) 'very moist on land point: ', I,J,SMC(I,1,J)
1764	CALL wrf_debug(10,message)
1765	endif
1766
1767	enddo
1768	enddo
1769
1770	!!! compute EMT, EM on global domain, and only on task 0.
1771
1772	#ifdef DM_PARALLEL
1773	IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
1774	#else
1775	IF (JDS .eq. JTS) THEN !! set unfailable condition for serial job
1776	#endif
1777
1778	ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
1779
1780	DO J=JDS,JDE-1
1781	TPH=SB+float(J-1)*DPH
1782	DXP=ERADDLMCOS(TPH)
1783	EMJ(J)= GRID%DT/( 4.DXP)DTAD
1784	EMTJ(J)=GRID%DT/(16.DXP)DTAD
1785	ENDDO
1786
1787	JA=0
1788	DO 161 J=3,5
1789	JA=JA+1
1790	KHLA(JA)=2
1791	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1792	161 EMT(JA)=EMTJ(J)
1793	DO 162 J=(JDE-1)-4,(JDE-1)-2
1794	JA=JA+1
1795	KHLA(JA)=2
1796	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1797	162 EMT(JA)=EMTJ(J)
1798	DO 163 J=6,(JDE-1)-5
1799	JA=JA+1
1800	KHLA(JA)=2
1801	KHHA(JA)=2+MOD(J,2)
1802	163 EMT(JA)=EMTJ(J)
1803	DO 164 J=6,(JDE-1)-5
1804	JA=JA+1
1805	KHLA(JA)=(IDE-1)-2
1806	KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1807	164 EMT(JA)=EMTJ(J)
1808
1809	! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
1810
1811	JA=0
1812	DO 171 J=3,5
1813	JA=JA+1
1814	KVLA(JA)=2
1815	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1816	171 EM(JA)=EMJ(J)
1817	DO 172 J=(JDE-1)-4,(JDE-1)-2
1818	JA=JA+1
1819	KVLA(JA)=2
1820	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1821	172 EM(JA)=EMJ(J)
1822	DO 173 J=6,(JDE-1)-5
1823	JA=JA+1
1824	KVLA(JA)=2
1825	KVHA(JA)=2+MOD(J+1,2)
1826	173 EM(JA)=EMJ(J)
1827	DO 174 J=6,(JDE-1)-5
1828	JA=JA+1
1829	KVLA(JA)=(IDE-1)-2
1830	KVHA(JA)=(IDE-1)-1-MOD(J,2)
1831	174 EM(JA)=EMJ(J)
1832
1833	696 continue
1834	ENDIF ! wrf_dm_on_monitor/serial job
1835
1836	call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,ISLTYP, &
1837	SM,SICE,STC,SMC,SH2O)
1838
1839	!! must be a better place to put this, but will eliminate "phantom"
1840	!! wind points here (no wind point on eastern boundary of odd numbered rows)
1841
1842	IF ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary
1843	write(message,*) 'zero phantom winds'
1844	CALL wrf_message(message)
1845	DO K=1,KDE-1
1846	DO J=JDS,JDE-1,2
1847	IF (J .ge. JTS .and. J .le. JTE) THEN
1848	u(IDE-1,J,K)=0.
1849	v(IDE-1,J,K)=0.
1850	ENDIF
1851	ENDDO
1852	ENDDO
1853	ENDIF
1854
1855	969 continue
1856
1857	DO j = jms, jme
1858	DO i = ims, ime
1859	fisx=max(fis(i,j),0.)
1860	Z0(I,J) =SM(I,J)Z0SEA+(1.-SM(I,J)) &
1861	& (0.Z0MAX+FISx FCM+Z0LAND)
1862	ENDDO
1863	ENDDO
1864
1865	write(message,*) 'Z0 over memory, leaving module_initialize_real'
1866	CALL wrf_message(message)
1867	DO J=JME,JMS,-((JME-JMS)/20+1)
1868	write(message,635) (Z0(I,J),I=IMS,IME,(IME-IMS)/14+1)
1869	CALL wrf_message(message)
1870	ENDDO
1871
1872
1873	endif ! on first_time check
1874
1875	write(message,*) 'leaving init_domain_nmm'
1876	CALL wrf_message( TRIM(message) )
1877	!
1878	write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, &
1879	& grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma
1880	CALL wrf_message( TRIM(message) )
1881	!==================================================================================
1882
1883	#define COPY_OUT
1884	#include <scalar_derefs.inc>
1885	RETURN
1886
1887	END SUBROUTINE init_domain_nmm
1888
1889	!------------------------------------------------------
1890
1891	SUBROUTINE define_nmm_vertical_coord ( LM, PTSGM, PT, PDTOP,HYBLEVS, &
1892	SG1,DSG1,SGML1, &
1893	SG2,DSG2,SGML2,DFL, DFRLG )
1894
1895	IMPLICIT NONE
1896
1897	! USE module_model_constants
1898
1899	!!! certain physical parameters here probably don't need to be defined, as defined
1900	!!! elsewhere within WRF. Done for initial testing purposes.
1901
1902	INTEGER :: LM, LPT2, L
1903	REAL :: PTSGM, PT, PL, PT2, PDTOP
1904	REAL :: RGOG, PSIG,PHYB,PHYBM
1905	REAL, PARAMETER :: Rd = 287.04 ! J deg{-1} kg{-1}
1906	REAL, PARAMETER :: CP=1004.6,GAMMA=.0065,PRF0=101325.,T0=288.
1907	REAL, PARAMETER :: g=9.81
1908
1909	REAL, DIMENSION(LM) :: DSG,DSG1,DSG2
1910	REAL, DIMENSION(LM) :: SGML1,SGML2
1911	REAL, DIMENSION(LM+1) :: SG1,SG2,HYBLEVS,DFL,DFRLG
1912
1913	CHARACTER(LEN=255) :: message
1914
1915	LPT2=LM+1
1916
1917	write(message,*) 'pt= ', pt
1918	CALL wrf_message(message)
1919
1920	DO L=LM+1,1,-1
1921	pl=HYBLEVS(L)*(101325.-pt)+pt
1922	if(pl.lt.ptSGm) LPT2=l
1923	ENDDO
1924
1925	IF(LPT2.lt.LM+1) THEN
1926	pt2=HYBLEVS(LPT2)*(101325.-pt)+pt
1927	ELSE
1928	pt2=pt
1929	ENDIF
1930
1931	write(message,) '** Sigma system starts at ',pt2,' Pa, from level ',LPT2
1932	CALL wrf_message(message)
1933
1934	pdtop=pt2-pt
1935
1936	write(message,*) 'allocating DSG,DSG1,DSG2 as ', LM
1937	CALL wrf_debug(10,message)
1938
1939	DSG=-99.
1940
1941	DO L=1,LM
1942	DSG(L)=HYBLEVS(L)- HYBLEVS(L+1)
1943	ENDDO
1944
1945	DSG1=0.
1946	DSG2=0.
1947
1948	DO L=LM,1,-1
1949
1950	IF(L.ge.LPT2) then
1951	DSG1(L)=DSG(L)
1952	ELSE
1953	DSG2(L)=DSG(L)
1954	ENDIF
1955
1956	ENDDO
1957
1958	SGML1=-99.
1959	SGML2=-99.
1960
1961	IF(LPT2.le.LM+1) THEN
1962
1963	DO L=LM+1,LPT2,-1
1964	SG2(L)=0.
1965	ENDDO
1966
1967	DO L=LPT2,2,-1
1968	SG2(L-1)=SG2(L)+DSG2(L-1)
1969	ENDDO
1970
1971	DO L=LPT2-1,1,-1
1972	SG2(L)=SG2(L)/SG2(1)
1973	ENDDO
1974	SG2(1)=1.
1975
1976	DO L=LPT2-1,1,-1
1977	DSG2(L)=SG2(L)-SG2(L+1)
1978	SGML2(l)=(SG2(l)+SG2(l+1))*0.5
1979	ENDDO
1980
1981	ENDIF
1982
1983	DO L=LM,LPT2,-1
1984	DSG2(L)=0.
1985	SGML2(L)=0.
1986	ENDDO
1987
1988	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1989
1990	SG1(LM+1)=0.
1991
1992	DO L=LM+1,LPT2,-1
1993	SG1(L-1)=SG1(L)+DSG1(L-1)
1994	ENDDO
1995
1996	DO L=LM,LPT2,-1
1997	SG1(L)=SG1(L)/SG1(LPT2-1)
1998	ENDDO
1999
2000	SG1(LPT2-1)=1.
2001
2002	do l=LPT2-2,1,-1
2003	SG1(l)=1.
2004	enddo
2005
2006
2007	DO L=LM,LPT2,-1
2008	DSG1(L)=SG1(L)-SG1(L+1)
2009	SGML1(L)=(SG1(L)+SG1(L+1))*0.5
2010	ENDDO
2011
2012	DO L=LPT2-1,1,-1
2013	DSG1(L)=0.
2014	SGML1(L)=1.
2015	ENDDO
2016
2017	1000 format('l,hyblevs,psig,SG1,SG2,phyb,phybm')
2018	1100 format(' ',i4,f7.4,f10.2,2f7.4,2f10.2)
2019
2020	write(message,1000)
2021	CALL wrf_debug(100,message)
2022
2023	do l=1,LM+1
2024	psig=HYBLEVS(L)*(101325.-pt)+pt
2025	phyb=SG1(l)pdtop+SG2(l)(101325.-pdtop-pt)+pt
2026	if(l.lt.LM+1) then
2027	phybm=SGML1(l)pdtop+SGML2(l)(101325.-pdtop-pt)+pt
2028	else
2029	phybm=-99.
2030	endif
2031
2032	write(message,1100) l,HYBLEVS(L),psig &
2033	,SG1(l),SG2(l),phyb,phybm
2034	CALL wrf_debug(100,message)
2035	enddo
2036
2037
2038	632 format(f9.6)
2039
2040	write(message,*) 'SG1'
2041	CALL wrf_debug(100,message)
2042	do L=LM+1,1,-1
2043	write(message,632) SG1(L)
2044	CALL wrf_debug(100,message)
2045	enddo
2046
2047	write(message,*) 'SG2'
2048	CALL wrf_debug(100,message)
2049	do L=LM+1,1,-1
2050	write(message,632) SG2(L)
2051	CALL wrf_debug(100,message)
2052	enddo
2053
2054	write(message,*) 'DSG1'
2055	CALL wrf_debug(100,message)
2056	do L=LM,1,-1
2057	write(message,632) DSG1(L)
2058	CALL wrf_debug(100,message)
2059	enddo
2060
2061	write(message,*) 'DSG2'
2062	CALL wrf_debug(100,message)
2063	do L=LM,1,-1
2064	write(message,632) DSG2(L)
2065	CALL wrf_debug(100,message)
2066	enddo
2067
2068	write(message,*) 'SGML1'
2069	CALL wrf_debug(100,message)
2070	do L=LM,1,-1
2071	write(message,632) SGML1(L)
2072	CALL wrf_debug(100,message)
2073	enddo
2074
2075	write(message,*) 'SGML2'
2076	CALL wrf_debug(100,message)
2077	do L=LM,1,-1
2078	write(message,632) SGML2(L)
2079	CALL wrf_debug(100,message)
2080	enddo
2081
2082	rgog=(rd*gamma)/g
2083	DO L=1,LM+1
2084	DFL(L)=gT0(1.-((pt+SG1(L)pdtop+SG2(L)(101325.-pt2)) &
2085	/101325.)**rgog)/gamma
2086	DFRLG(L)=DFL(L)/g
2087	write(message,*) 'L, DFL(L): ', L, DFL(L)
2088	CALL wrf_debug(100,message)
2089	ENDDO
2090
2091	END SUBROUTINE define_nmm_vertical_coord
2092
2093	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2094	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2095	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2096
2097	SUBROUTINE compute_nmm_surfacep ( TERRAIN_HGT_T, Z3D_IN, PRESS3D_IN, T3D_IN &
2098	&, psfc_out,generic &
2099	&, IDS,IDE,JDS,JDE,KDS,KDE &
2100	&, IMS,IME,JMS,JME,KMS,KME &
2101	&, ITS,ITE,JTS,JTE,KTS,KTE )
2102
2103
2104	IMPLICIT NONE
2105
2106	real, allocatable:: dum2d(:,:),DUM2DB(:,:)
2107
2108	integer :: IDS,IDE,JDS,JDE,KDS,KDE
2109	integer :: IMS,IME,JMS,JME,KMS,KME
2110	integer :: ITS,ITE,JTS,JTE,KTS,KTE,Ilook,Jlook
2111	integer :: I,J,II,generic,L,KINSERT,K,bot_lev,LL
2112	integer :: IHE(JMS:JME),IHW(JMS:JME), loopinc,iloopinc
2113
2114	real :: TERRAIN_HGT_T(IMS:IME,JMS:JME)
2115	real :: Z3D_IN(IMS:IME,JMS:JME,generic)
2116	real :: T3D_IN(IMS:IME,JMS:JME,generic)
2117	real :: PRESS3D_IN(IMS:IME,JMS:JME,generic)
2118	real :: PSFC_IN(IMS:IME,JMS:JME),TOPO_IN(IMS:IME,JMS:JME)
2119	real :: psfc_out(IMS:IME,JMS:JME),rincr(IMS:IME,JMS:JME)
2120	real :: dif1,dif2,dif3,dif4,dlnpdz,BOT_INPUT_HGT,BOT_INPUT_PRESS,dpdz,rhs
2121	real :: zin(generic),pin(generic)
2122
2123	character (len=255) :: message
2124
2125	logical :: DEFINED_PSFC(IMS:IME,JMS:JME), DEFINED_PSFCB(IMS:IME,JMS:JME)
2126
2127	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2128
2129	Ilook=25
2130	Jlook=25
2131
2132	DO j = JMS, JME
2133	IHE(J)=MOD(J+1,2)
2134	IHW(J)=IHE(J)-1
2135	ENDDO
2136
2137	DO J=JMS,JME
2138	DO I=IMS,IME
2139	DEFINED_PSFC(I,J)=.FALSE.
2140	DEFINED_PSFCB(I,J)=.FALSE.
2141	IF (PRESS3D_IN(I,J,1) .ne. 200100.) THEN
2142	PSFC_IN(I,J)=PRESS3D_IN(I,J,1)
2143	TOPO_IN(I,J)=Z3D_IN(I,J,1)
2144	ELSE
2145	PSFC_IN(I,J)=PRESS3D_IN(I,J,2)
2146	TOPO_IN(I,J)=Z3D_IN(I,J,2)
2147	ENDIF
2148	ENDDO
2149	ENDDO
2150
2151	! input surface pressure smoothing over the ocean - still needed for NAM?
2152
2153	II_loop: do II=1,8
2154
2155	CYCLE II_loop
2156
2157	do J=JTS+1,min(JTE,JDE-1)-1
2158	do I=ITS+1,min(ITE,IDE-1)-1
2159	rincr(I,J)=0.
2160
2161	if (PSFC_IN(I,J) .gt. 100000. .and. &
2162	PSFC_IN(I+IHE(J),J+1) .gt. 100000. .and. &
2163	PSFC_IN(I+IHE(J),J-1) .gt. 100000. .and. &
2164	PSFC_IN(I+IHW(J),J+1) .gt. 100000. .and. &
2165	PSFC_IN(I+IHW(J),J-1) .gt. 100000. ) then
2166
2167	dif1=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J+1))
2168	dif2=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J-1))
2169	dif3=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J+1))
2170	dif4=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J-1))
2171
2172	if (max(dif1,dif2,dif3,dif4) .lt. 200. .and. TOPO_IN(I,J).le. 0.5 .and. &
2173	TOPO_IN(I+IHE(J),J+1) .le. 0.5 .and. &
2174	TOPO_IN(I+IHW(J),J+1) .le. 0.5 .and. &
2175	TOPO_IN(I+IHE(J),J-1) .le. 0.5 .and. &
2176	TOPO_IN(I+IHW(J),J-1) .lt. 0.5) then
2177
2178	rincr(I,J)=0.125( 4.PSFC_IN(I,J)+ &
2179	PSFC_IN(I+IHE(J),J+1)+PSFC_IN(I+IHE(J),J-1)+ &
2180	PSFC_IN(I+IHW(J),J+1)+PSFC_IN(I+IHW(J),J-1) ) &
2181	- PSFC_IN(I,J)
2182
2183	! if (rincr(I,J) .ne. 0 .and. abs(rincr(I,J)) .gt. 20.) then
2184	! write(message,*) 'II, I,J,rincr: ', II, I,J,rincr(I,J)
2185	! CALL wrf_message(message)
2186	! endif
2187
2188	endif
2189	endif
2190
2191	ENDDO
2192	ENDDO
2193
2194	DO J=JTS+1,min(JTE,JDE-1)-1
2195	DO I=ITS+1,min(ITE,IDE-1)-1
2196	PSFC_IN(I,J)=PSFC_IN(I,J) + rincr(I,J)
2197	ENDDO
2198	ENDDO
2199
2200	! write(message,*) ' -------------------------------------------------- '
2201	! CALL wrf_message(message)
2202
2203	end do II_loop
2204
2205	ALLOCATE(DUM2D(IMS:IME,JMS:JME))
2206
2207	DO J=JMS,JME
2208	DO I=IMS,IME
2209	DUM2D(I,J)=-9.
2210	END DO
2211	END DO
2212
2213	DO J=JTS,min(JTE,JDE-1)
2214	I_loop: DO I=ITS,min(ITE,IDE-1)
2215
2216	IF (PSFC_IN(I,J) .lt. 0.1) THEN
2217	write(message,*) 'QUITTING BECAUSE I,J, PSFC_IN: ', I,J,PSFC_IN(I,J)
2218	call wrf_error_fatal(message)
2219	ENDIF
2220
2221	BOT_INPUT_PRESS=PSFC_IN(I,J)
2222	BOT_INPUT_HGT=TOPO_IN(I,J)
2223
2224	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2225
2226	write(message,*) ' TERRAIN_HGT_T: ', I,J, TERRAIN_HGT_T(I,J)
2227	CALL wrf_message(message)
2228	write(message,*) ' PSFC_IN, TOPO_IN: ', &
2229	I, J, PSFC_IN(I,J),TOPO_IN(I,J)
2230	CALL wrf_message(message)
2231
2232	DO L=1,generic
2233	write(message,*) ' L,PRESS3D_IN, Z3D_IN: ', &
2234	I,J,L, PRESS3D_IN(I,J,L),Z3D_IN(I,J,L)
2235	CALL wrf_debug(10,message)
2236	END DO
2237	ENDIF
2238
2239	DO L=2,generic-1
2240
2241	IF ( PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
2242	Z3D_IN(I,J,L) .lt. TERRAIN_HGT_T(I,J) .AND. &
2243	Z3D_IN(I,J,L+1) .gt. TERRAIN_HGT_T(I,J) ) THEN
2244
2245	BOT_INPUT_PRESS=PRESS3D_IN(i,j,L)
2246	BOT_INPUT_HGT=Z3D_IN(I,J,L)
2247
2248	! IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2249	! write(message,*) 'BOT_INPUT_PRESS, BOT_INPUT_HGT NOW : ', &
2250	! Ilook,Jlook, BOT_INPUT_PRESS, BOT_INPUT_HGT
2251	! CALL wrf_message(message)
2252	! ENDIF
2253
2254	ENDIF
2255	END DO
2256
2257	!!!!!!!!!!!!!!!!!!!!!! START HYDRO CHECK
2258
2259	IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
2260	(PSFC_IN(I,J) .gt. PRESS3D_IN(i,j,2) .OR. &
2261	TOPO_IN(I,J) .lt. Z3D_IN(I,J,2)) ) THEN ! extrapolate downward
2262
2263	IF (J .eq. JTS .AND. I .eq. ITS) THEN
2264	write(message,*) 'hydro check - should only be for isobaric input'
2265	CALL wrf_message(message)
2266	ENDIF
2267
2268	IF (Z3D_IN(I,J,2) .ne. TOPO_IN(I,J)) THEN
2269	dpdz=(PRESS3D_IN(i,j,2)-PSFC_IN(I,J))/(Z3D_IN(I,J,2)-TOPO_IN(I,J))
2270	rhs=-9.81((PRESS3D_IN(i,j,2)+ PSFC_IN(I,J))/2.)/(287.04 T3D_IN(I,J,2))
2271
2272	IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
2273	IF (dpdz .lt. 1.05rhs .OR. dpdz .gt. 0.95rhs) THEN
2274	write(message,*) 'I,J,P(2),Psfc,Z(2),Zsfc: ', &
2275	I,J,PRESS3D_IN(i,j,2),PSFC_IN(I,J),Z3D_IN(I,J,2),TOPO_IN(I,J)
2276	IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) CALL wrf_debug(50,message)
2277	CYCLE I_loop
2278	ENDIF
2279
2280	ENDIF
2281
2282	ELSE ! z(2) equals TOPO_IN
2283
2284	IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
2285	! write(message,*) 'all equal at I,J: ', I,J
2286	! CALL wrf_message(message)
2287	ELSE
2288	! write(message,*) 'heights equal, pressures not: ', &
2289	! PRESS3D_IN(i,j,2), PSFC_IN(I,J)
2290	! CALL wrf_message(message)
2291	CYCLE I_loop
2292	ENDIF
2293
2294	ENDIF
2295
2296	IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
2297	IF (PRESS3D_IN(i,j,2) .lt. PSFC_IN(I,J) .and. &
2298	Z3D_IN(I,J,2) .lt. TOPO_IN(I,J)) THEN
2299	! write(message,*) 'surface data mismatch(a) at I,J: ', I,J
2300	! CALL wrf_message(message)
2301	CYCLE I_loop
2302	ELSEIF (PRESS3D_IN(i,j,2) .gt. PSFC_IN(I,J) .AND. &
2303	Z3D_IN(I,J,2) .gt. TOPO_IN(I,J)) THEN
2304	! write(message,*) 'surface data mismatch(b) at I,J: ', I,J
2305	! CALL wrf_message(message)
2306	CYCLE I_loop
2307	ENDIF
2308	ENDIF
2309	ENDIF
2310
2311	!!!!!!! loop over a few more levels
2312
2313	DO L=3,6
2314	IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
2315	(((PSFC_IN(I,J)-PRESS3D_IN(i,j,L)) .lt. 400.) .OR. &
2316	TOPO_IN(I,J) .lt. Z3D_IN(I,J,L))) then
2317
2318	IF (Z3D_IN(I,J,L) .ne. TOPO_IN(I,J)) THEN
2319	dpdz=(PRESS3D_IN(i,j,L)-PSFC_IN(I,J))/ &
2320	(Z3D_IN(I,J,L)-TOPO_IN(I,J))
2321	rhs=-9.81*((PRESS3D_IN(i,j,L)+ PSFC_IN(I,J))/2.)/ &
2322	(287.04*T3D_IN(I,J,L))
2323	IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
2324	IF (dpdz .lt. 1.05rhs .or. dpdz .gt. 0.95rhs) THEN
2325	write(message,*) 'I,J,L,Piso,Psfc,Ziso,Zsfc: ', &
2326	I,J,L,PRESS3D_IN(i,j,L),PSFC_IN(I,J),&
2327	Z3D_IN(I,J,L),TOPO_IN(I,J)
2328	IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) &
2329	CALL wrf_debug(50,message)
2330	CYCLE I_loop
2331	ENDIF
2332	ENDIF
2333	ELSE
2334	IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
2335	! write(message,*) 'all equal at I,J: ', I,J
2336	! CALL wrf_message(message)
2337	ELSE
2338	CYCLE I_loop
2339	ENDIF
2340	ENDIF
2341	ENDIF
2342
2343	IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
2344	IF (PRESS3D_IN(i,j,L) .lt. PSFC_IN(I,J) .AND. &
2345	Z3D_IN(I,J,L) .lt. TOPO_IN(I,J)) THEN
2346	CYCLE I_loop
2347	ELSEIF (PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
2348	Z3D_IN(I,J,L) .gt. TOPO_IN(I,J)) THEN
2349	CYCLE I_loop
2350	ENDIF
2351	ENDIF
2352	END DO
2353	!!!!!!!!!!!!!!!!!!!!!! END HYDRO CHECK
2354
2355	IF (TERRAIN_HGT_T(I,J) .eq. BOT_INPUT_HGT ) THEN
2356	dum2d(I,J)=BOT_INPUT_PRESS
2357
2358	IF (BOT_INPUT_HGT .ne. 0. .and. (BOT_INPUT_HGT-INT(BOT_INPUT_HGT) .ne. 0.) ) THEN
2359	write(message,*) 'with BOT_INPUT_HGT: ', BOT_INPUT_HGT, &
2360	'set dum2d to bot_input_pres: ', I,J,dum2d(I,J)
2361	CALL wrf_message(message)
2362	ENDIF
2363
2364	IF (dum2d(I,J) .lt. 50000. .OR. dum2d(I,J) .gt. 109000.) THEN
2365	write(message,*) 'bad dum2d(a): ', I,J,DUM2D(I,J)
2366	CALL wrf_message(message)
2367	ENDIF
2368
2369	ELSEIF (TERRAIN_HGT_T(I,J) .lt. BOT_INPUT_HGT ) THEN
2370
2371	! target is below lowest possible input...extrapolate
2372
2373	IF ( BOT_INPUT_PRESS-PRESS3D_IN(I,J,2) .gt. 500. ) THEN
2374	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,2)) ) / &
2375	(BOT_INPUT_HGT-Z3D_IN(i,j,2))
2376	IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2377	write(message,*) 'I,J,dlnpdz(a): ', I,J,dlnpdz
2378	CALL wrf_message(message)
2379	ENDIF
2380
2381	ELSE
2382
2383	!! thin layer and/or just have lowest level - difference with 3rd level data
2384	IF ( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,3)) .gt. 290. ) THEN
2385
2386	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,3)) ) / &
2387	(BOT_INPUT_HGT-Z3D_IN(i,j,3))
2388
2389	IF (I .eq. Ilook .and. J .eq. Jlook) then
2390	write(message,*) 'p diff: ', BOT_INPUT_PRESS, PRESS3D_IN(i,j,3)
2391	CALL wrf_message(message)
2392	write(message,*) 'z diff: ', BOT_INPUT_HGT, Z3D_IN(i,j,3)
2393	CALL wrf_message(message)
2394	ENDIF
2395
2396	ELSE
2397
2398	!! Loop up to level 7 looking for a sufficiently thick layer
2399
2400	FIND_THICK: DO LL=4,7
2401	IF( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,LL)) .gt. 290.) THEN
2402	dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,LL)) ) / &
2403	(BOT_INPUT_HGT-Z3D_IN(i,j,LL))
2404	EXIT FIND_THICK
2405	ENDIF
2406	END DO FIND_THICK
2407
2408	ENDIF
2409
2410	ENDIF
2411
2412	dum2d(I,J)= exp(log(BOT_INPUT_PRESS) + dlnpdz * &
2413	(TERRAIN_HGT_T(I,J) - BOT_INPUT_HGT) )
2414
2415	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2416	write(message,*) 'bad dum2d(b): ', I,J,DUM2D(I,J)
2417	CALL wrf_message(message)
2418	write(message,*) 'BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T, BOT_INPUT_HGT: ', &
2419	BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T(I,J), BOT_INPUT_HGT
2420	CALL wrf_message(message)
2421	write(message,*) 'Z3D_IN: ', Z3D_IN(I,J,1:10)
2422	CALL wrf_message(message)
2423	write(message,*) 'PRESS3D_IN: ', PRESS3D_IN(I,J,1:10)
2424	CALL wrf_message(message)
2425	ENDIF
2426
2427	ELSE ! target level bounded by input levels
2428
2429	DO L=2,generic-1
2430	IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
2431	TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
2432	dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
2433	(Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
2434	dum2d(I,J)= log(PRESS3D_IN(i,j,l)) + &
2435	dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
2436	dum2d(i,j)=exp(dum2d(i,j))
2437	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2438	write(message,*) 'bad dum2d(c): ', I,J,DUM2D(I,J)
2439	CALL wrf_message(message)
2440	ENDIF
2441	ENDIF
2442	ENDDO
2443
2444	!!! account for situation where BOT_INPUT_HGT < TERRAIN_HGT_T < Z3D_IN(:,2,:)
2445	IF (dum2d(I,J) .eq. -9 .AND. BOT_INPUT_HGT .lt. TERRAIN_HGT_T(I,J) &
2446	.AND. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(I,J,2)) then
2447
2448	IF (mod(I,50) .eq. 0 .AND. mod(J,50) .eq. 0) THEN
2449	write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
2450	I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
2451	CALL wrf_message(message)
2452	ENDIF
2453
2454	dlnpdz= (log(PSFC_IN(i,j))-log(PRESS3D_IN(i,j,2)) ) / &
2455	(TOPO_IN(i,j)-Z3D_IN(i,j,2))
2456	dum2d(I,J)= log(PSFC_IN(i,j)) + &
2457	dlnpdz * (TERRAIN_HGT_T(I,J) - TOPO_IN(i,j) )
2458	dum2d(i,j)= exp(dum2d(i,j))
2459	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2460	write(message,*) 'bad dum2d(d): ', I,J,DUM2D(I,J)
2461	CALL wrf_message(message)
2462	ENDIF
2463	ENDIF
2464
2465	IF (dum2d(I,J) .eq. -9.) THEN
2466	write(message,*) 'must have flukey situation in new ', I,J
2467	CALL wrf_message(message)
2468	write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
2469	I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
2470	CALL wrf_message(message)
2471
2472	DO L=1,generic-1
2473	IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
2474	! problematic with HGT_M substitution for "input" surface height?
2475	dum2d(i,j)=PRESS3D_IN(I,J,L)
2476	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2477	write(message,*) 'bad dum2d(e): ', I,J,DUM2D(I,J)
2478	CALL wrf_message(message)
2479	ENDIF
2480	ENDIF
2481	ENDDO
2482
2483	IF ( TERRAIN_HGT_T(I,J) .eq. TOPO_IN(I,J)) THEN
2484	dum2d(I,J)=PSFC_IN(I,J)
2485	IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2486	write(message,*) 'bad dum2d(f): ', I,J,DUM2D(I,J)
2487	CALL wrf_message(message)
2488	ENDIF
2489	write(message,*) 'matched input topo, psfc: ', I,J,TOPO_IN(I,J),PSFC_IN(I,J)
2490	CALL wrf_message(message)
2491	ENDIF
2492
2493	IF (dum2d(I,J) .eq. -9.) THEN
2494	CALL wrf_error_fatal("quitting due to undefined surface pressure")
2495	ENDIF
2496	ENDIF
2497
2498	DEFINED_PSFC(I,J)=.TRUE.
2499
2500	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2501	write(message,*) 'newstyle psfc: ', I,J,dum2d(I,J)
2502	CALL wrf_message(message)
2503	ENDIF
2504
2505	ENDIF
2506
2507	ENDDO I_loop
2508	ENDDO
2509
2510	write(message,*) 'psfc points (new style)'
2511	CALL wrf_message(message)
2512	loopinc=max( (JTE-JTS)/20,1)
2513	iloopinc=max( (ITE-ITS)/10,1)
2514
2515	DO J=min(JTE,JDE-1),JTS,-loopinc
2516	write(message,633) (dum2d(I,J)/100.,I=ITS,min(ITE,IDE-1),iloopinc)
2517	END DO
2518
2519	633 format(35(f5.0,1x))
2520
2521	write(message,*) 'PSFC extremes (new style)'
2522	CALL wrf_message(message)
2523	write(message,*) minval(dum2d,MASK=DEFINED_PSFC),maxval(dum2d,MASK=DEFINED_PSFC)
2524	CALL wrf_message(message)
2525
2526	! IF (minval(dum2d,MASK=DEFINED_PSFC) .lt. 50000. .or. maxval(dum2d,MASK=DEFINED_PSFC) .gt. 108000.) THEN
2527
2528	DO J=JTS,min(JTE,JDE-1)
2529	DO I=ITS,min(ITE,IDE-1)
2530
2531	IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .lt. 50000. ) THEN
2532	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2533	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2534	CALL wrf_debug(2,message)
2535	ELSE
2536	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2537	ENDIF
2538	ENDIF
2539
2540	IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .gt. 108000. ) THEN
2541	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2542	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2543	CALL wrf_debug(2,message)
2544	ELSE
2545	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2546	ENDIF
2547	ENDIF
2548
2549	END DO
2550	END DO
2551
2552
2553
2554	!! "traditional" isobaric only approach ------------------------------------------------
2555
2556	ALLOCATE (DUM2DB(IMS:IME,JMS:JME))
2557	DO J=JMS,JME
2558	DO I=IMS,IME
2559	DUM2DB(I,J)=-9.
2560	END DO
2561	END DO
2562
2563	DO J=JTS,min(JTE,JDE-1)
2564	DO I=ITS,min(ITE,IDE-1)
2565
2566	IF (TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,2)) THEN ! targ below lowest
2567
2568	IF ( abs(PRESS3D_IN(i,j,2)-PRESS3D_IN(i,j,3)) .gt. 290.) THEN
2569	dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,3)) ) / &
2570	(Z3D_IN(i,j,2)-Z3D_IN(i,j,3))
2571	ELSE
2572	dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,4)) ) / &
2573	(Z3D_IN(i,j,2)-Z3D_IN(i,j,4))
2574	ENDIF
2575
2576	DUM2DB(I,J)= exp( log(PRESS3D_IN(i,j,2)) + dlnpdz * &
2577	(TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,2)) )
2578
2579	IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2580	write(message,*) 'I,K, trad: dlnpdz, press_in(2), terrain_t, Z3D_IN(2): ', I,J,dlnpdz, &
2581	PRESS3D_IN(i,j,2), TERRAIN_HGT_T(I,J), Z3D_IN(i,j,2)
2582	CALL wrf_message(message)
2583	ENDIF
2584
2585	DEFINED_PSFCB(i,j)=.true.
2586
2587	ELSEIF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,2)) THEN ! target level bounded by input levels
2588
2589	DO L=2,generic-1
2590	IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
2591	TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
2592
2593	dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
2594	(Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
2595
2596	DUM2DB(I,J)= log(PRESS3D_IN(i,j,l)) + &
2597	dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
2598	DUM2DB(i,j)=exp(DUM2DB(i,j))
2599
2600	DEFINED_PSFCB(i,j)=.true.
2601
2602	IF (DUM2DB(I,J) .lt. 13000.) THEN
2603	write(message,*) 'I,J,L,terrain,Z3d(L),z3d(L+1),p3d(L),p3d(l+1): ', I,J,L, &
2604	TERRAIN_HGT_T(I,J),Z3D_IN(I,J,L),Z3D_IN(I,J,L+1),PRESS3D_IN(I,J,L), &
2605	PRESS3D_IN(I,J,L+1)
2606	CALL wrf_error_fatal(message)
2607	ENDIF
2608	ENDIF
2609	ENDDO
2610
2611	ELSEIF (TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,2)) THEN
2612	DUM2DB(i,j)=PRESS3D_IN(I,J,2)
2613	DEFINED_PSFCB(i,j)=.true.
2614	ENDIF
2615
2616	IF (DUM2DB(I,J) .eq. -9.) THEN
2617	write(message,*) 'must have flukey situation in trad ', I,J
2618	CALL wrf_message(message)
2619	DO L=1,generic-1
2620	IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
2621	DUM2DB(i,j)=PRESS3D_IN(I,J,L)
2622	DEFINED_PSFCB(i,j)=.true.
2623	ENDIF
2624	ENDDO
2625	ENDIF
2626
2627	IF (DUM2DB(I,J) .eq. -9.) THEN
2628	write(message,*) 'HOPELESS PSFC, I QUIT'
2629	CALL wrf_error_fatal(message)
2630	ENDIF
2631
2632	if (I .eq. Ilook .and. J .eq. Jlook) THEN
2633	write(message,*) ' traditional psfc: ', I,J,DUM2DB(I,J)
2634	CALL wrf_message(message)
2635	ENDIF
2636
2637	ENDDO
2638	ENDDO
2639
2640	write(message,*) 'psfc points (traditional)'
2641	CALL wrf_message(message)
2642	DO J=min(JTE,JDE-1),JTS,-loopinc
2643	write(message,633) (DUM2DB(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
2644	CALL wrf_message(message)
2645	ENDDO
2646
2647	write(message,*) 'PSFC extremes (traditional)'
2648	CALL wrf_message(message)
2649	write(message,*) minval(DUM2DB,MASK=DEFINED_PSFCB),maxval(DUM2DB,MASK=DEFINED_PSFCB)
2650	CALL wrf_message(message)
2651
2652	DO J=JTS,min(JTE,JDE-1)
2653	DO I=ITS,min(ITE,IDE-1)
2654
2655	IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .lt. 50000. ) THEN
2656	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2657	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2658	CALL wrf_debug(2,message)
2659	ELSE
2660	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2661	ENDIF
2662	ENDIF
2663
2664	IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .gt. 108000. ) THEN
2665	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2666	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2667	CALL wrf_debug(2,message)
2668	ELSE
2669	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2670	ENDIF
2671	ENDIF
2672
2673	! IF (DEFINED_PSFCB(I,J) .AND. ( dum2db(I,J) .lt. 50000. .OR. dum2db(I,J) .gt. 108000. )) THEN
2674	! IF (TERRAIN_HGT_T(I,J) .gt. -10. .and. TERRAIN_HGT_T(I,J) .lt. 5000.) THEN
2675	! write(0,*) 'I,J, terrain_hgt_t, dum2db: ', I,J, terrain_hgt_t(I,J),dum2db(I,J)
2676	! CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2677	! ELSE
2678	! WRITE(message,*) 'surface pressure allowed because surface height is extreme value of: ', TERRAIN_HGT_T(I,J)
2679	! CALL wrf_debug(2,message)
2680	! ENDIF
2681	! ENDIF
2682
2683	ENDDO
2684	ENDDO
2685
2686	!!!!! end traditional
2687
2688	DO J=JTS,min(JTE,JDE-1)
2689	DO I=ITS,min(ITE,IDE-1)
2690	IF (DEFINED_PSFCB(I,J) .and. DEFINED_PSFC(I,J)) THEN
2691
2692	IF ( abs(dum2d(I,J)-DUM2DB(I,J)) .gt. 400.) THEN
2693	write(message,*) 'BIG DIFF I,J, dum2d, DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
2694	CALL wrf_message(message)
2695	ENDIF
2696
2697	!! do we have enough confidence in new style to give it more than 50% weight?
2698	psfc_out(I,J)=0.5*(dum2d(I,J)+DUM2DB(I,J))
2699
2700	ELSEIF (DEFINED_PSFC(I,J)) THEN
2701	psfc_out(I,J)=dum2d(I,J)
2702	ELSEIF (DEFINED_PSFCB(I,J)) THEN
2703	psfc_out(I,J)=DUM2DB(I,J)
2704	ELSE
2705	write(message,*) 'I,J,dum2d,DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
2706	CALL wrf_message(message)
2707	write(message,*) 'I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J): ', I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J)
2708	CALL wrf_message(message)
2709	call wrf_error_fatal("psfc_out completely undefined")
2710	ENDIF
2711
2712	IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2713	write(message,*) ' combined psfc: ', I,J,psfc_out(I,J)
2714	CALL wrf_message(message)
2715	ENDIF
2716
2717	IF (PSFC_OUT(I,J) .lt. 50000. ) THEN
2718	IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2719	WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2720	CALL wrf_debug(2,message)
2721	ELSE
2722	write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
2723	CALL wrf_debug(2,message)
2724	write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
2725	CALL wrf_debug(2,message)
2726	write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
2727	CALL wrf_debug(2,message)
2728	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2729	ENDIF
2730	ENDIF
2731
2732	IF (PSFC_OUT(I,J) .gt. 108000. ) THEN
2733	IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2734	WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2735	CALL wrf_debug(2,message)
2736	ELSE
2737	write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
2738	CALL wrf_debug(2,message)
2739	write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
2740	CALL wrf_debug(2,message)
2741	write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
2742	CALL wrf_debug(2,message)
2743	CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2744	ENDIF
2745	ENDIF
2746
2747	ENDDO
2748	ENDDO
2749
2750	deallocate(dum2d,dum2db)
2751
2752	END SUBROUTINE compute_nmm_surfacep
2753
2754	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2755	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2756
2757	SUBROUTINE compute_3d_pressure(psfc_out,SGML1,SGML2,pdtop,pt &
2758	&, pd,p3d_out &
2759	&, IDS,IDE,JDS,JDE,KDS,KDE &
2760	&, IMS,IME,JMS,JME,KMS,KME &
2761	&, ITS,ITE,JTS,JTE,KTS,KTE )
2762
2763
2764	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
2765	INTEGER :: IMS,IME,JMS,JME,KMS,KME
2766	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
2767
2768	REAL, INTENT(IN) :: psfc_out(IMS:IME,JMS:JME)
2769	REAL, INTENT(IN) :: SGML1(KDE),SGML2(KDE),pdtop,pt
2770
2771	REAL, INTENT(OUT):: p3d_out(IMS:IME,JMS:JME,KDS:KDE-1)
2772	REAL, INTENT(OUT):: PD(IMS:IME,JMS:JME)
2773
2774	CHARACTER (len=255) :: message
2775
2776	! write(message,*) 'pdtop, pt, psfc_out(1,1): ', pdtop, pt, psfc_out(1,1)
2777	! CALL wrf_message(message)
2778
2779	DO J=JTS,min(JTE,JDE-1)
2780	DO I=ITS,min(ITE,IDE-1)
2781	PD(I,J)=psfc_out(I,J)-PDTOP-PT
2782	ENDDO
2783	ENDDO
2784
2785	DO J=JTS,min(JTE,JDE-1)
2786	DO K=KDS,KDE-1
2787	DO I=ITS,min(ITE,IDE-1)
2788	p3d_out(I,J,K)=PD(I,J)SGML2(K)+PDTOPSGML1(K)+PT
2789
2790	IF (p3d_out(I,J,K) .ge. psfc_out(I,J) .or. p3d_out(I,J,K) .le. pt) THEN
2791	write(message,*) 'I,K,J,p3d_out: ', I,K,J,p3d_out(I,J,K)
2792	CALL wrf_error_fatal(message)
2793	ENDIF
2794
2795	ENDDO
2796	ENDDO
2797	ENDDO
2798
2799	END SUBROUTINE compute_3d_pressure
2800
2801	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2802	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2803	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2804
2805	SUBROUTINE interp_press2press_lin(press_in,press_out, &
2806	data_in, data_out,generic &
2807	&, extrapolate,ignore_lowest,TFIELD &
2808	&, IDS,IDE,JDS,JDE,KDS,KDE &
2809	&, IMS,IME,JMS,JME,KMS,KME &
2810	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
2811
2812	! Interpolates data from one set of pressure surfaces to
2813	! another set of pressures
2814
2815	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
2816	INTEGER :: IMS,IME,JMS,JME,KMS,KME
2817	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
2818	INTEGER :: internal_time
2819
2820	! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
2821	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
2822	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
2823	! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
2824	REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
2825	REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
2826	LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest, TFIELD
2827	LOGICAL :: col_smooth
2828
2829	INTEGER :: i,j
2830	INTEGER :: k,kk
2831	REAL :: desired_press
2832	REAL :: dvaldlnp,dlnp,tadiabat,tiso
2833
2834	REAL, PARAMETER :: ADIAFAC=9.81/1004.
2835	REAL, PARAMETER :: TSTEXTRAPFAC=.0065
2836
2837
2838
2839	DO K=KMS,KME
2840	DO J=JMS,JME
2841	DO I=IMS,IME
2842	DATA_OUT(I,J,K)=-99999.9
2843	ENDDO
2844	ENDDO
2845	ENDDO
2846
2847	IF (ignore_lowest) then
2848	LMIN=2
2849	ELSE
2850	LMIN=1
2851	ENDIF
2852
2853	DO j = JTS, min(JTE,JDE-1)
2854	test_i: DO i = ITS, min(ITE,IDE-1)
2855
2856	IF (internal_time_loop .gt. 1) THEN
2857	IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
2858	I .ne. IDS .and. I .ne. IDE-1 ) THEN
2859	!! not on boundary
2860	CYCLE test_i
2861	ENDIF
2862	ENDIF
2863
2864
2865	col_smooth=.false.
2866
2867	output_loop: DO k = KDS,KDE-1
2868
2869	desired_press = press_out(i,j,k)
2870
2871	if (K .gt. KDS) then
2872	if (TFIELD .and. col_smooth .and. desired_press .le. press_in(i,j,LMIN) &
2873	.and. press_out(i,j,k-1) .ge. press_in(i,j,LMIN)) then
2874	MAX_SMOOTH=K
2875	! write(message,*) 'I,J, MAX_SMOOTH: ', I,J, MAX_SMOOTH
2876	! CALL wrf_debug(100,message)
2877	endif
2878	endif
2879
2880	! keep track of where the extrapolation begins
2881
2882	IF (desired_press .GT. press_in(i,j,LMIN)) THEN
2883	IF (TFIELD .and. K .eq. 1 .and. (desired_press - press_in(i,j,LMIN)) .gt. 3000.) then
2884	col_smooth=.TRUE. ! due to large extrapolation distance
2885	ENDIF
2886
2887
2888	IF ((desired_press - press_in(i,j,LMIN)).LT. 50.) THEN ! 0.5 mb
2889	data_out(i,j,k) = data_in(i,j,LMIN)
2890	ELSE
2891	IF (extrapolate) THEN
2892	! Extrapolate downward because desired P level is below
2893	! the lowest level in our input data. Extrapolate using simple
2894	! 1st derivative of value with respect to ln P for the bottom 2
2895	! input layers.
2896
2897	! Add a check to make sure we are not using the gradient of
2898	! a very thin layer
2899
2900	if (TFIELD) then
2901	tiso=0.5*(data_in(i,j,1)+data_in(i,j,2))
2902	endif
2903
2904
2905	IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 500.) THEN ! likely isobaric data
2906	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
2907	dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+1)) / dlnp
2908	ELSE ! assume terrain following
2909	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+5))
2910	dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+5)) / dlnp
2911	ENDIF
2912	data_out(i,j,k) = data_in(i,j,LMIN) + dvaldlnp * &
2913	( log(desired_press)-log(press_in(i,j,LMIN)) )
2914
2915	if (TFIELD .and. data_out(i,j,k) .lt. tiso-0.2) then
2916
2917	! restrict slope to -1K/10 hPa
2918	dvaldlnp=max(dvaldlnp, -1.0/ &
2919	log( press_in(i,j,LMIN) / &
2920	( press_in(i,j,LMIN)-1000.) ))
2921
2922	data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
2923	( log(desired_press)-log(press_in(i,j,LMIN)) )
2924
2925	elseif (TFIELD .and. data_out(i,j,k) .gt. tiso+0.2) then
2926
2927	! restrict slope to +0.8K/10 hPa
2928	dvaldlnp=min(dvaldlnp, 0.8/ &
2929	log( press_in(i,j,LMIN) / &
2930	( press_in(i,j,LMIN)-1000.) ))
2931
2932	data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
2933	( log(desired_press)-log(press_in(i,j,LMIN)) )
2934
2935	endif
2936
2937	ELSE
2938	data_out(i,j,k) = data_in(i,j,LMIN)
2939	ENDIF
2940	ENDIF
2941	ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
2942	IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
2943	data_out(i,j,k) = data_in(i,j,generic)
2944	ELSE
2945	IF (extrapolate) THEN
2946	! Extrapolate upward
2947	IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
2948	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
2949	dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-1))/dlnp
2950	ELSE
2951	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
2952	dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-2))/dlnp
2953	ENDIF
2954	data_out(i,j,k) = data_in(i,j,generic) + &
2955	dvaldlnp * (log(desired_press)-log(press_in(i,j,generic)))
2956	ELSE
2957	data_out(i,j,k) = data_in(i,j,generic)
2958	ENDIF
2959	ENDIF
2960	ELSE
2961	! We can trap between two levels and linearly interpolate
2962
2963	input_loop: DO kk = LMIN, generic-1
2964	IF (desired_press .EQ. press_in(i,j,kk) )THEN
2965	data_out(i,j,k) = data_in(i,j,kk)
2966	EXIT input_loop
2967	ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
2968	(desired_press .GT. press_in(i,j,kk+1)) ) THEN
2969
2970	! do trapped in lnp
2971
2972	dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
2973	dvaldlnp = (data_in(i,j,kk)-data_in(i,j,kk+1))/dlnp
2974	data_out(i,j,k) = data_in(i,j,kk+1)+ &
2975	dvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1)))
2976
2977	EXIT input_loop
2978	ENDIF
2979
2980	ENDDO input_loop
2981	ENDIF
2982	ENDDO output_loop
2983
2984	if (col_smooth) then
2985	do K=max(KDS,MAX_SMOOTH-4),MAX_SMOOTH+4
2986	data_out(I,J,K)=0.5*(data_out(I,J,K)+data_out(I,J,K+1))
2987	enddo
2988	endif
2989
2990	ENDDO test_i
2991	ENDDO
2992	END SUBROUTINE interp_press2press_lin
2993
2994	SUBROUTINE wind_adjust(press_in,press_out, &
2995	U_in, V_in,U_out,V_out &
2996	&, generic,depth_replace &
2997	&, IDS,IDE,JDS,JDE,KDS,KDE &
2998	&, IMS,IME,JMS,JME,KMS,KME &
2999	&, ITS,ITE,JTS,JTE,KTS,KTE )
3000
3001	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3002	INTEGER :: IMS,IME,JMS,JME,KMS,KME
3003	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
3004	INTEGER :: MAXLIN,MAXLOUT
3005
3006	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
3007	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
3008	REAL, INTENT(IN) :: U_in(IMS:IME,JMS:JME,generic)
3009	REAL, INTENT(IN) :: V_in(IMS:IME,JMS:JME,generic)
3010	REAL, INTENT(INOUT) :: U_out(IMS:IME,KMS:KME,JMS:JME)
3011	REAL, INTENT(INOUT) :: V_out(IMS:IME,KMS:KME,JMS:JME)
3012	REAL :: p1d_in(generic)
3013	REAL :: p1d_out(KDS:KDE-1)
3014
3015
3016	DO j = JTS, min(JTE,JDE-1)
3017	DO i = ITS, min(ITE,IDE-1)
3018
3019	! IF (press_out(I,J,1) .lt. press_in(I,J,2)) then
3020	IF( (press_in(I,J,2)-press_out(I,J,1)) .gt. 200.) then
3021
3022	U_out(I,1,J)=U_in(I,J,2)
3023	V_out(I,1,J)=V_in(I,J,2)
3024
3025	INLOOP: DO L=2,generic
3026	p1d_in(L)=-9999.
3027	IF ( (press_in(I,J,2)-press_in(I,J,L)) .lt. depth_replace) THEN
3028	p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
3029	MAXLIN=L
3030	ELSE
3031	p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
3032	EXIT INLOOP
3033	ENDIF
3034	END DO INLOOP
3035
3036	OUTLOOP: DO L=KDS,KDE-1
3037	p1d_out(L)=-9999.
3038	IF ( (press_out(I,J,1)-press_out(I,J,L)) .lt. depth_replace) THEN
3039	p1d_out(L)=(press_out(I,J,1)-press_out(I,J,L))
3040	MAXLOUT=L
3041	ELSE
3042	EXIT OUTLOOP
3043	ENDIF
3044	END DO OUTLOOP
3045
3046	DO L=1,MAXLOUT
3047	ptarg=p1d_out(L)
3048
3049	FINDLOOP: DO LL=2,MAXLIN
3050
3051	if (p1d_in(LL) .lt. ptarg .and. p1d_in(LL+1) .gt. ptarg) then
3052
3053	dlnp=log(p1d_in(LL))-log(p1d_in(LL+1))
3054	dudlnp=(U_in(I,J,LL)-U_in(I,J,LL+1))/dlnp
3055	dvdlnp=(V_in(I,J,LL)-V_in(I,J,LL+1))/dlnp
3056	U_out(I,L,J)=U_in(I,J,LL)+dudlnp*(log(ptarg)-log(p1d_in(LL)))
3057	V_out(I,L,J)=V_in(I,J,LL)+dvdlnp*(log(ptarg)-log(p1d_in(LL)))
3058
3059	EXIT FINDLOOP
3060	endif
3061
3062	END DO FINDLOOP
3063	END DO ! MAXLOUT loop
3064
3065
3066	ENDIF
3067
3068	ENDDO
3069	ENDDO
3070
3071
3072
3073	END SUBROUTINE wind_adjust
3074	!--------------------------------------------------------------------
3075
3076	SUBROUTINE interp_press2press_log(press_in,press_out, &
3077	data_in, data_out, generic &
3078	&, extrapolate,ignore_lowest &
3079	&, IDS,IDE,JDS,JDE,KDS,KDE &
3080	&, IMS,IME,JMS,JME,KMS,KME &
3081	&, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
3082
3083	! Interpolates ln(data) from one set of pressure surfaces to
3084	! another set of pressures
3085
3086	INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3087	INTEGER :: IMS,IME,JMS,JME,KMS,KME
3088	INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
3089	INTEGER :: internal_time
3090
3091	! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
3092	REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
3093	REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
3094	! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
3095	! REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
3096	REAL :: data_in(IMS:IME,JMS:JME,generic)
3097	REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
3098	LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest
3099
3100	INTEGER :: i,j
3101	INTEGER :: k,kk
3102	REAL :: desired_press
3103	REAL :: dlnvaldlnp,dlnp
3104
3105
3106	DO K=1,generic
3107	DO j = JTS, min(JTE,JDE-1)
3108	DO i = ITS, min(ITE,IDE-1)
3109	DATA_IN(I,J,K)=max(DATA_in(I,J,K),1.e-13)
3110	ENDDO
3111	ENDDO
3112	ENDDO
3113
3114	DO K=KMS,KME
3115	DO J=JMS,JME
3116	DO I=IMS,IME
3117	DATA_OUT(I,J,K)=-99999.9
3118	ENDDO
3119	ENDDO
3120	ENDDO
3121
3122	IF (ignore_lowest) then
3123	LMIN=2
3124	ELSE
3125	LMIN=1
3126	ENDIF
3127
3128	DO j = JTS, min(JTE,JDE-1)
3129	test_i: DO i = ITS, min(ITE,IDE-1)
3130
3131	IF (internal_time .gt. 1) THEN
3132	IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
3133	I .ne. IDS .and. I .ne. IDE-1 ) THEN
3134	!! not on boundary
3135	CYCLE test_i
3136	ENDIF
3137	ENDIF
3138
3139
3140	output_loop: DO k = KDS,KDE-1
3141
3142	desired_press = press_out(i,j,k)
3143
3144	IF (desired_press .GT. press_in(i,j,LMIN)) THEN
3145
3146	IF ((desired_press - press_in(i,j,LMIN)).LT. 10.) THEN ! 0.1 mb
3147	data_out(i,j,k) = data_in(i,j,LMIN)
3148	ELSE
3149	IF (extrapolate) THEN
3150	! Extrapolate downward because desired P level is below
3151	! the lowest level in our input data. Extrapolate using simple
3152	! 1st derivative of value with respect to ln P for the bottom 2
3153	! input layers.
3154
3155	! Add a check to make sure we are not using the gradient of
3156	! a very thin layer
3157
3158	IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 100.) THEN
3159	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
3160	dlnvaldlnp = ( log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+1)) ) / dlnp
3161
3162	ELSE
3163
3164	dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+2))
3165	dlnvaldlnp = (log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+2))) / dlnp
3166
3167	ENDIF
3168
3169	data_out(i,j,k) = exp(log(data_in(i,j,LMIN)) + dlnvaldlnp * &
3170	( log(desired_press)-log(press_in(i,j,LMIN))))
3171	ELSE
3172	data_out(i,j,k) = data_in(i,j,LMIN)
3173	ENDIF
3174	ENDIF
3175	ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
3176	IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
3177	data_out(i,j,k) = data_in(i,j,generic)
3178	ELSE
3179	IF (extrapolate) THEN
3180	! Extrapolate upward
3181	IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
3182	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
3183	dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-1)))/dlnp
3184	ELSE
3185	dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
3186	dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-2)))/dlnp
3187	ENDIF
3188	data_out(i,j,k) = exp(log(data_in(i,j,generic)) + &
3189	dlnvaldlnp * (log(desired_press)-log(press_in(i,j,generic))))
3190	ELSE
3191	data_out(i,j,k) = data_in(i,j,generic)
3192	ENDIF
3193	ENDIF
3194	ELSE
3195	! We can trap between two levels and linearly interpolate
3196
3197	input_loop: DO kk = LMIN, generic-1
3198	IF (desired_press .EQ. press_in(i,j,kk) )THEN
3199	data_out(i,j,k) = data_in(i,j,kk)
3200	EXIT input_loop
3201	ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
3202	(desired_press .GT. press_in(i,j,kk+1)) ) THEN
3203
3204	! do trapped in lnp
3205
3206	dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
3207	dlnvaldlnp = (log(data_in(i,j,kk))-log(data_in(i,j,kk+1)))/dlnp
3208	data_out(i,j,k) = exp(log(data_in(i,j,kk+1))+ &
3209	dlnvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1))))
3210
3211	EXIT input_loop
3212
3213	ENDIF
3214
3215	ENDDO input_loop
3216	ENDIF
3217	ENDDO output_loop
3218	ENDDO test_i
3219	ENDDO
3220	END SUBROUTINE interp_press2press_log
3221
3222	!-------------------------------------------------------------------
3223	SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
3224	ids , ide , jds , jde , kds , kde , &
3225	ims , ime , jms , jme , kms , kme , &
3226	its , ite , jts , jte , kts , kte )
3227
3228	IMPLICIT NONE
3229
3230	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3231	ims , ime , jms , jme , kms , kme , &
3232	its , ite , jts , jte , kts , kte
3233
3234	LOGICAL , INTENT(IN) :: wrt_liquid
3235
3236	! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
3237	! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
3238	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: p , t
3239	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(INOUT) :: rh
3240	REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(OUT) :: q
3241
3242	! Local vars
3243
3244	INTEGER :: i , j , k
3245
3246	REAL :: ew , q1 , t1
3247
3248	REAL, PARAMETER :: T_REF = 0.0
3249	REAL, PARAMETER :: MW_AIR = 28.966
3250	REAL, PARAMETER :: MW_VAP = 18.0152
3251
3252	REAL, PARAMETER :: A0 = 6.107799961
3253	REAL, PARAMETER :: A1 = 4.436518521e-01
3254	REAL, PARAMETER :: A2 = 1.428945805e-02
3255	REAL, PARAMETER :: A3 = 2.650648471e-04
3256	REAL, PARAMETER :: A4 = 3.031240396e-06
3257	REAL, PARAMETER :: A5 = 2.034080948e-08
3258	REAL, PARAMETER :: A6 = 6.136820929e-11
3259
3260	REAL, PARAMETER :: ES0 = 6.1121
3261
3262	REAL, PARAMETER :: C1 = 9.09718
3263	REAL, PARAMETER :: C2 = 3.56654
3264	REAL, PARAMETER :: C3 = 0.876793
3265	REAL, PARAMETER :: EIS = 6.1071
3266	REAL :: RHS
3267	REAL, PARAMETER :: TF = 273.16
3268	REAL :: TK
3269
3270	REAL :: ES
3271	REAL :: QS
3272	REAL, PARAMETER :: EPS = 0.622
3273	REAL, PARAMETER :: SVP1 = 0.6112
3274	REAL, PARAMETER :: SVP2 = 17.67
3275	REAL, PARAMETER :: SVP3 = 29.65
3276	REAL, PARAMETER :: SVPT0 = 273.15
3277
3278	! This subroutine computes mixing ratio (q, kg/kg) from basic variables
3279	! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
3280	! The reference temperature (t_ref, C) is used to describe the temperature
3281	! at which the liquid and ice phase change occurs.
3282
3283	DO k = kts , kte
3284	DO j = jts , MIN ( jde-1 , jte )
3285	DO i = its , MIN (ide-1 , ite )
3286	rh(i,j,k) = MIN ( MAX ( rh(i,j,k) , 1. ) , 100. )
3287	END DO
3288	END DO
3289	END DO
3290
3291	IF ( wrt_liquid ) THEN
3292	DO k = kts , kte
3293	DO j = jts , MIN ( jde-1 , jte )
3294	DO i = its , MIN (ide-1 , ite )
3295	es=svp110.EXP(svp2*(t(i,j,k)-svpt0)/(t(i,j,k)-svp3))
3296	qs=eps*es/(p(i,j,k)/100.-es)
3297	q(i,j,k)=MAX(.01rh(i,j,k)qs,0.0)
3298	END DO
3299	END DO
3300	END DO
3301
3302	ELSE
3303	DO k = kts , kte
3304	DO j = jts , MIN ( jde-1 , jte )
3305	DO i = its , MIN (ide-1 , ite )
3306
3307	t1 = t(i,j,k) - 273.16
3308
3309	! Obviously dry.
3310
3311	IF ( t1 .lt. -200. ) THEN
3312	q(i,j,k) = 0
3313
3314	ELSE
3315
3316	! First compute the ambient vapor pressure of water
3317
3318	IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
3319	ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
3320
3321	ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
3322	ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
3323
3324	ELSE
3325	tk = t(i,j,k)
3326	rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
3327	c3 * (1. - tk / tf) + alog10(eis)
3328	ew = 10. ** rhs
3329
3330	END IF
3331
3332	! Now sat vap pres obtained compute local vapor pressure
3333
3334	ew = MAX ( ew , 0. ) * rh(i,j,k) * 0.01
3335
3336	! Now compute the specific humidity using the partial vapor
3337	! pressures of water vapor (ew) and dry air (p-ew). The
3338	! constants assume that the pressure is in hPa, so we divide
3339	! the pressures by 100.
3340
3341	q1 = mw_vap * ew
3342	q1 = q1 / (q1 + mw_air * (p(i,j,k)/100. - ew))
3343
3344	q(i,j,k) = q1 / (1. - q1 )
3345
3346	END IF
3347
3348	END DO
3349	END DO
3350	END DO
3351
3352	END IF
3353
3354	END SUBROUTINE rh_to_mxrat
3355
3356	!--=------------------------------------------------------------------
3357
3358	SUBROUTINE boundary_smooth(h, landmask, grid, nsmth , nrow &
3359	&, IDS,IDE,JDS,JDE,KDS,KDE &
3360	&, IMS,IME,JMS,JME,KMS,KME &
3361	&, ITS,ITE,JTS,JTE,KTS,KTE )
3362
3363	implicit none
3364
3365	TYPE (domain) :: grid
3366
3367	integer :: IDS,IDE,JDS,JDE,KDS,KDE
3368	integer :: IMS,IME,JMS,JME,KMS,KME
3369	integer :: ITS,ITE,JTS,JTE,KTS,KTE
3370	integer:: ihw(JDS:JDE-1),ihe(JDS:JDE-1),nsmth,nrow
3371	real:: h(IMS:IME,JMS:JME),landmask(IMS:IME,JMS:JME)
3372	real :: h_old(IMS:IME,JMS:JME)
3373	real :: hbms(IDS:IDE-1,JDS:JDE-1)
3374	real :: hse(IDS:IDE-1,JDS:JDE-1)
3375	real :: hne(IDS:IDE-1,JDS:JDE-1)
3376	integer :: IPS,IPE,JPS,JPE,KPS,KPE
3377	integer :: ihl, ihh, m2l, ibas,jmelin
3378	integer :: I,J,KS,IOFFSET,JSTART,JEND
3379	character (len=255) :: message
3380
3381	ips=its
3382	ipe=ite
3383	jps=jts
3384	jpe=jte
3385	kps=kts
3386	kpe=kte
3387
3388	do j= JTS,min(JTE,JDE-1)
3389	ihw(J)=-mod(J,2)
3390	ihe(j)=ihw(J)+1
3391	end do
3392
3393	do J=JTS,min(JTE,JDE-1)
3394	do I=ITS,min(ITE,IDE-1)
3395	hbms(I,J)=landmask(I,J)
3396	enddo
3397	enddo
3398
3399	jmelin=(JDE-1)-nrow+1
3400	ibas=nrow/2
3401	m2l=mod(nrow,2)
3402
3403	do j=jts,min(jte,jde-1)
3404	ihl=ibas+mod(j,2)+m2l*mod(J+1,2)
3405	ihh=(IDE-1)-ibas-m2l*mod(J+1,2)
3406	do i=its,min(ite,ide-1)
3407	if (I .ge. ihl .and. I .le. ihh .and. J .ge. nrow .and. J .le. jmelin) then
3408	hbms(I,J)=0.
3409	endif
3410	end do
3411	end do
3412
3413	634 format(30(f2.0,1x))
3414
3415	do KS=1,nsmth
3416
3417	grid%ht_gc=h
3418	#ifdef DM_PARALLEL
3419	# include "HALO_NMM_MG.inc"
3420	#endif
3421	h=grid%ht_gc
3422	h_old=grid%ht_gc
3423
3424	do J=JTS,min(JTE,JDE-1)
3425	do I=ITS, min(ITE,IDE-1)
3426	if (I .ge. (IDS+mod(J,2)) .and. J .gt. JDS .and. J .lt. JDE-1 .and. I .lt. IDE-1) then
3427	h(i,j)= ( h_old(i+ihe(j),j+1) + h_old(i+ihw(j),j-1) + h_old(i+ihe(j),j-1) + h_old(i+ihw(j),j+1) - &
3428	4. h_old(i,j) )hbms(i,j)*0.125+h_old(i,j)
3429	endif
3430
3431	enddo
3432	enddo
3433
3434	! special treatment for four corners
3435
3436	if (hbms(1,1) .eq. 1 .and. ITS .le. 1 .and. JTS .le. 1) then
3437	h(1,1)=0.75h(1,1)+0.125h(1+ihe(1),2)+ &
3438	0.0625*(h(2,1)+h(1,3))
3439	endif
3440
3441	if (hbms(IDE-1,1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTS .le. 1) then
3442	h(IDE-1,1)=0.75h(IDE-1,1)+0.125h(IDE-1+ihw(1),2)+ &
3443	0.0625*(h(IDE-1-1,1)+h(IDE-1,3))
3444	endif
3445
3446	if (hbms(1,JDE-1) .eq. 1 .and. ITS .le. 1 .and. JTE .ge. JDE-2) then
3447	h(1,JDE-1)=0.75h(1,JDE-1)+0.125h(1+ihe(JDE-1),JDE-1-1)+ &
3448	0.0625*(h(2,JDE-1)+h(1,JDE-1-2))
3449	endif
3450
3451	if (hbms(IDE-1,JDE-1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTE .ge. JDE-2) then
3452	h(IDE-1,JDE-1)=0.75h(IDE-1,JDE-1)+0.125h(IDE-1+ihw(JDE-1),JDE-1-1)+ &
3453	0.0625*(h(IDE-1-1,JDE-1)+h(IDE-1,JDE-1-2))
3454	endif
3455
3456	do J=JMS,JME
3457	do I=IMS,IME
3458	grid%ht_gc(I,J)=h(I,J)
3459	enddo
3460	enddo
3461	#ifdef DM_PARALLEL
3462	# include "HALO_NMM_MG.inc"
3463	#endif
3464	do J=JMS,JME
3465	do I=IMS,IME
3466	h(I,J)=grid%ht_gc(I,J)
3467	enddo
3468	enddo
3469
3470
3471	! S bound
3472	if (JTS .eq. JDS) then
3473	J=JTS
3474
3475	do I=ITS,ITE
3476	if (I .ge. IDS+1 .and. I .le. IDE-2) then
3477	if (hbms(I,J) .eq. 1) then
3478	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
3479	endif
3480	endif
3481	enddo
3482
3483	endif
3484
3485	! N bound
3486	if (JTE .eq. JDE) then
3487	J=JDE-1
3488	write(message,*) 'DOING N BOUND SMOOTHING for J= ', J
3489	CALL wrf_message(message)
3490	do I=ITS,min(ITE,IDE-1)
3491	if (hbms(I,J) .eq. 1 .and. I .ge. IDS+1 .and. I .le. IDE-2) then
3492	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
3493	endif
3494	enddo
3495	endif
3496
3497	! W bound
3498	if (ITS .eq. IDS) then
3499	I=ITS
3500	do J=JTS,min(JTE,JDE-1)
3501	if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
3502	h(I,J)=0.75h(I,J)+0.125(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
3503	endif
3504	enddo
3505	endif
3506
3507	! E bound
3508	if (ITE .eq. IDE) then
3509	write(message,*) 'DOING E BOUND SMOOTHING for I= ', min(ITE,IDE-1)
3510	CALL wrf_message(message)
3511	I=min(ITE,IDE-1)
3512	do J=JTS,min(JTE,JDE-1)
3513	if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
3514	h(I,J)=0.75h(I,J)+0.125(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
3515	endif
3516	enddo
3517	endif
3518
3519	enddo ! end ks loop
3520
3521	do J=JMS,JME
3522	do I=IMS,IME
3523	grid%ht_gc(I,J)=h(I,J)
3524	enddo
3525	enddo
3526	#ifdef DM_PARALLEL
3527	#include "HALO_NMM_MG.inc"
3528	#endif
3529	do J=JMS,JME
3530	do I=IMS,IME
3531	h(I,J)=grid%ht_gc(I,J)
3532	enddo
3533	enddo
3534
3535	! extra smoothing along inner boundary
3536
3537	if (JTS .eq. JDS) then
3538	if (ITE .eq. IDE) then
3539	IOFFSET=1
3540	else
3541	IOFFSET=0
3542	endif
3543	! Southern Boundary
3544	do i=its,min(ITE,IDE-1)-IOFFSET
3545	h(i,2)=0.25*(h(i,1)+h(i+1,1)+ &
3546	h(i,3)+h(i+1,3))
3547	enddo
3548	endif
3549
3550
3551	if (JTE .eq. JDE) then
3552	if (ITE .eq. IDE) then
3553	IOFFSET=1
3554	else
3555	IOFFSET=0
3556	endif
3557	do i=its,min(ITE,IDE-1)-IOFFSET
3558	h(i,(JDE-1)-1)=0.25*(h(i,(JDE-1)-2)+h(i+1,(JDE-1)-2)+ &
3559	h(i,JDE-1)+h(i+1,JDE-1))
3560	enddo
3561	endif
3562
3563	if (JTS .eq. 1) then
3564	JSTART=4
3565	else
3566	JSTART=JTS+mod(JTS,2) ! needs to be even
3567	endif
3568
3569	if (JTE .eq. JDE) then
3570	JEND=(JDE-1)-3
3571	else
3572	JEND=JTE
3573	endif
3574
3575	if (ITS .eq. IDS) then
3576
3577	! Western Boundary
3578	do j=JSTART,JEND,2
3579	h(1,j)=0.25*(h(1,j-1)+h(2,j-1)+ &
3580	h(1,j+1)+h(2,j+1))
3581
3582	enddo
3583	endif
3584
3585
3586	if (ITE .eq. IDE) then
3587	! Eastern Boundary
3588	do j=JSTART,JEND,2
3589	h((IDE-1)-1,j)=0.25*(h((IDE-1)-1,j-1)+h((IDE-1),j-1)+ &
3590	h((IDE-1)-1,j+1)+h((IDE-1),j+1))
3591	enddo
3592	endif
3593
3594
3595	END SUBROUTINE boundary_smooth
3596
3597	!--------------------------------------------------------------------
3598
3599	SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
3600	ids , ide , jds , jde , kds , kde , &
3601	ims , ime , jms , jme , kms , kme , &
3602	its , ite , jts , jte , kts , kte )
3603
3604	! Linrarly in time interpolate data to a current valid time. The data is
3605	! assumed to come in "monthly", valid at the 15th of every month.
3606
3607	IMPLICIT NONE
3608
3609	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3610	ims , ime , jms , jme , kms , kme , &
3611	its , ite , jts , jte , kts , kte
3612
3613	CHARACTER (LEN=24) , INTENT(IN) :: date_str
3614	REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
3615	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
3616
3617	! Local vars
3618
3619	INTEGER :: i , j , l
3620	INTEGER , DIMENSION(0:13) :: middle
3621	INTEGER :: target_julyr , target_julday , target_date
3622	INTEGER :: julyr , julday , int_month, next_month
3623	REAL :: gmt
3624	CHARACTER (LEN=4) :: yr
3625	CHARACTER (LEN=2) :: mon , day15
3626
3627
3628	WRITE(day15,FMT='(I2.2)') 15
3629	DO l = 1 , 12
3630	WRITE(mon,FMT='(I2.2)') l
3631	CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
3632	middle(l) = julyr*1000 + julday
3633	END DO
3634
3635	l = 0
3636	middle(l) = middle( 1) - 31
3637
3638	l = 13
3639	middle(l) = middle(12) + 31
3640
3641	CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
3642	target_date = target_julyr * 1000 + target_julday
3643	find_month : DO l = 0 , 12
3644	IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
3645	DO j = jts , MIN ( jde-1 , jte )
3646	DO i = its , MIN (ide-1 , ite )
3647	int_month = MOD ( l , 12 )
3648	IF ( int_month .EQ. 0 ) int_month = 12
3649
3650	IF (int_month == 12) THEN
3651	next_month=1
3652	ELSE
3653	next_month=int_month+1
3654	ENDIF
3655
3656	field_out(i,j) = ( field_in(i,j,next_month) * ( target_date - middle(l) ) + &
3657	field_in(i,j,int_month ) * ( middle(l+1) - target_date ) ) / &
3658	( middle(l+1) - middle(l) )
3659	END DO
3660	END DO
3661	EXIT find_month
3662	END IF
3663	END DO find_month
3664	END SUBROUTINE monthly_interp_to_date
3665
3666	!---------------------------------------------------------------------
3667	SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
3668	ids , ide , jds , jde , kds , kde , &
3669	ims , ime , jms , jme , kms , kme , &
3670	its , ite , jts , jte , kts , kte )
3671
3672	! Plow through each month, find the max, min values for each i,j.
3673
3674	IMPLICIT NONE
3675
3676	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3677	ims , ime , jms , jme , kms , kme , &
3678	its , ite , jts , jte , kts , kte
3679
3680	REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
3681	REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
3682
3683	! Local vars
3684
3685	INTEGER :: i , j , l
3686	REAL :: minner , maxxer
3687
3688	DO j = jts , MIN(jde-1,jte)
3689	DO i = its , MIN(ide-1,ite)
3690	minner = field_in(i,j,1)
3691	maxxer = field_in(i,j,1)
3692	DO l = 2 , 12
3693	IF ( field_in(i,j,l) .LT. minner ) THEN
3694	minner = field_in(i,j,l)
3695	END IF
3696	IF ( field_in(i,j,l) .GT. maxxer ) THEN
3697	maxxer = field_in(i,j,l)
3698	END IF
3699	END DO
3700	field_min(i,j) = minner
3701	field_max(i,j) = maxxer
3702	END DO
3703	END DO
3704
3705	END SUBROUTINE monthly_min_max
3706
3707	!-----------------------------------------------------------------------
3708
3709	SUBROUTINE reverse_vert_coord ( field, start_z, end_z &
3710	&, IDS,IDE,JDS,JDE,KDS,KDE &
3711	&, IMS,IME,JMS,JME,KMS,KME &
3712	&, ITS,ITE,JTS,JTE,KTS,KTE )
3713
3714	IMPLICIT NONE
3715
3716	INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3717	ims , ime , jms , jme , kms , kme , &
3718	its , ite , jts , jte , kts , kte, &
3719	start_z, end_z
3720
3721	REAL, INTENT(INOUT) :: field(IMS:IME,JMS:JME,end_z)
3722	! local
3723
3724	INTEGER :: I,J,L
3725	REAL, ALLOCATABLE :: dum3d(:,:,:)
3726
3727	allocate(dum3d(IMS:IME,JMS:JME,end_z))
3728
3729	DO L=start_z,end_z
3730	DO J=jts,min(jte,jde-1)
3731	DO I=its,min(ite,ide-1)
3732	dum3d(I,J,L)=field(I,J,end_z-L+start_z)
3733	END DO
3734	END DO
3735	END DO
3736
3737	DO L=start_z,end_z
3738	DO J=jts,min(jte,jde-1)
3739	DO I=its,min(ite,ide-1)
3740	field(I,J,L)=dum3d(I,J,L)
3741	END DO
3742	END DO
3743	END DO
3744
3745	DEALLOCATE(dum3d)
3746
3747	END SUBROUTINE reverse_vert_coord
3748
3749
3750	!--------------------------------------------------------------------
3751
3752	SUBROUTINE compute_nmm_levels(ninterface, ptop, eta_levels)
3753
3754	USE module_model_constants
3755
3756	IMPLICIT NONE
3757
3758	character(len=132):: message
3759	integer :: ninterface,Lthick,L
3760	real, parameter:: gamma=.0065
3761	real, parameter:: t_stand=288.
3762	real, parameter:: p_stand=101325.
3763
3764	real :: maxdz_compute, ptop
3765	real :: plower,pupper,tlay, sum
3766
3767	real :: eta_levels(ninterface)
3768	real, allocatable:: Z(:)
3769	real, allocatable:: deta_levels_spline(:)
3770
3771	logical:: print_pbl_warn
3772
3773	!----------------------------------------------------
3774
3775	allocate(Z(ninterface))
3776	allocate(deta_levels_spline(ninterface-1))
3777
3778	CALL compute_eta_spline(ninterface-1,deta_levels_spline,ptop)
3779
3780	sum=0.
3781	DO L=1,ninterface-1
3782	sum=sum+deta_levels_spline(L)
3783	ENDDO
3784
3785	eta_levels(1)=1.00
3786
3787	DO L=2,ninterface
3788	eta_levels(L)=eta_levels(L-1)-deta_levels_spline(L-1)
3789	ENDDO
3790
3791	eta_levels(ninterface)=0.00
3792
3793	DO L=2,ninterface-1
3794	eta_levels(L)=0.5(eta_levels(L))+0.25(eta_levels(L-1)+eta_levels(L+1))
3795	ENDDO
3796
3797	Z(1)=0.
3798	maxdz_compute=0.
3799	print_pbl_warn=.false.
3800
3801	DO L=2,ninterface
3802	tlay=max( t_stand-gamma*Z(L-1), 216.5)
3803	plower=ptop+(p_stand-ptop)*eta_levels(L-1)
3804	pupper=ptop+(p_stand-ptop)*eta_levels(L)
3805	Z(L)=Z(L-1)+(tlayr_d/g)(log(plower)-log(pupper))
3806
3807	if (plower .gt. 85000. .and. pupper .lt. 85000. .and. L .lt. 10) then
3808	print_pbl_warn=.true.
3809	endif
3810
3811	write(message,*) 'L, eta(l), pupper, Z(L): ', L, eta_levels(L),pupper,Z(L)
3812	CALL wrf_debug(100,message)
3813
3814	if (Z(L)-Z(L-1) .gt. maxdz_compute) then
3815	Lthick=L
3816	endif
3817
3818	maxdz_compute=max(maxdz_compute,Z(L)-Z(L-1))
3819	ENDDO
3820
3821	if (print_pbl_warn) then
3822	write(message,*) 'WARNING - PBL MAY BE POORLY RESOLVED WITH NUMBER OF VERTICAL LEVELS'
3823	CALL wrf_message(message)
3824	write(message,*) ' - CONSIDER INCREASING THE VERTICAL RESOLUTION'
3825	CALL wrf_message(message)
3826	endif
3827
3828	write(message,*) 'thickest layer was: ', maxdz_compute , 'meters thick at level: ', Lthick
3829	CALL wrf_message(message)
3830
3831	END SUBROUTINE compute_nmm_levels
3832
3833	!---------------------------
3834
3835	SUBROUTINE compute_eta_spline(LM, dsg, ptop)
3836
3837	IMPLICIT NONE
3838
3839	real:: dsg(LM), ptop, sum, rsum
3840	real, allocatable:: xold(:),dold(:)
3841	real, allocatable:: xnew(:),sgm(:)
3842	real, allocatable:: pps(:),qqs(:),y2s(:)
3843	integer nlev,LM,L,KOLD
3844
3845	IF (LM .ge. 46) THEN
3846	KOLD=9
3847	allocate(xold(KOLD))
3848	allocate(dold(KOLD))
3849
3850	xold(1)=.00
3851	dold(1)=.006
3852	xold(2)=.13
3853	dold(2)=.009
3854	xold(3)=.19
3855	dold(3)=.012
3856	xold(4)=.30
3857	dold(4)=.036
3858	xold(5)=.42
3859	dold(5)=.041
3860	xold(6)=.56
3861	dold(6)=.040
3862	xold(7)=.69
3863	dold(7)=.018
3864
3865	if (ptop .ge. 2000.) then
3866	xold(8)=.90
3867	dold(8)=.012
3868	xold(9)=1.0
3869	dold(9)=.006
3870	else
3871	xold(8)=.90
3872	dold(8)=.008
3873	xold(9)=1.0
3874	dold(9)=.003
3875	endif
3876
3877	ELSE
3878
3879	KOLD=8
3880	allocate(xold(KOLD))
3881	allocate(dold(KOLD))
3882
3883	xold(1)=.00
3884	dold(1)=.006
3885	xold(2)=.18
3886	dold(2)=.015
3887	xold(3)=.32
3888	dold(3)=.035
3889	xold(4)=.50
3890	dold(4)=.040
3891	xold(5)=.68
3892	dold(5)=.030
3893	xold(6)=.75
3894	dold(6)=.017
3895	xold(7)=.85
3896	dold(7)=.012
3897
3898	if (ptop .ge. 2000.) then
3899	xold(8)=1.0
3900	dold(8)=.013
3901	else
3902	xold(8)=1.0
3903	dold(8)=.008
3904	endif
3905
3906	ENDIF
3907
3908	allocate(xnew(lm))
3909	allocate(sgm(lm))
3910	allocate(pps(lm))
3911	allocate(qqs(lm))
3912	allocate(y2s(lm))
3913
3914	DO L=1,LM
3915	xnew(l)=float(l-1)/float(lm-1)
3916	ENDDO
3917
3918	y2s=0.
3919
3920	CALL spline(kold,kold,xold,dold,y2s,lm,xnew,dsg,pps,qqs)
3921
3922	sum=0.
3923	DO l=1,lm
3924	sum=sum+dsg(l)
3925	ENDDO
3926
3927	rsum=1./sum
3928	sgm(1)=0.
3929
3930	DO L=1,lm-1
3931	dsg(l)=dsg(l)*rsum
3932	sgm(l+1)=sgm(l)+dsg(l)
3933	ENDDO
3934	sgm(lm+1)=1.
3935	dsg(lm)=sgm(lm+1)-sgm(lm)
3936
3937	END SUBROUTINE compute_eta_spline
3938
3939	! -------------------------------------------------------------------
3940
3941	subroutine spline(JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
3942	! ********************************************************************
3943	! * *
3944	! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
3945	! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
3946	! * *
3947	! * PROGRAMER Z. JANJIC *
3948	! * *
3949	! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
3950	! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
3951	! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
3952	! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
3953	! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
3954	! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
3955	! * SPECIFIED. *
3956	! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
3957	! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
3958	! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
3959	! * AND LE XOLD(NOLD). *
3960	! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
3961	! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
3962	! * *
3963	! ********************************************************************
3964	!
3965	! LOG:
3966	!
3967	! PYLE - June 2007 - eliminated use of GO TO statements.
3968	!
3969	!-----------------------------------------------------------------------
3970	IMPLICIT NONE
3971	!-----------------------------------------------------------------------
3972	INTEGER,INTENT(IN) :: JTBX,NNEW,NOLD
3973	REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
3974	REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
3975	REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
3976	!
3977	INTEGER :: K,K1,K2,KOLD,NOLDM1
3978	REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
3979	& ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
3980	!-----------------------------------------------------------------------
3981	NOLDM1=NOLD-1
3982
3983	DXL=XOLD(2)-XOLD(1)
3984	DXR=XOLD(3)-XOLD(2)
3985	DYDXL=(YOLD(2)-YOLD(1))/DXL
3986	DYDXR=(YOLD(3)-YOLD(2))/DXR
3987	RTDXC=0.5/(DXL+DXR)
3988
3989	P(1)= RTDXC(6.(DYDXR-DYDXL)-DXL*Y2(1))
3990	Q(1)=-RTDXC*DXR
3991
3992	first_loop: DO K=3,NOLD-1
3993	IF(NOLD==3) exit first_loop ! should be impossible to execute
3994	DXL=DXR
3995	DYDXL=DYDXR
3996	DXR=XOLD(K+1)-XOLD(K)
3997	DYDXR=(YOLD(K+1)-YOLD(K))/DXR
3998	DXC=DXL+DXR
3999	DEN=1./(DXL*Q(K-2)+DXC+DXC)
4000	P(K-1)= DEN(6.(DYDXR-DYDXL)-DXL*P(K-2))
4001	Q(K-1)=-DEN*DXR
4002	END DO first_loop
4003
4004	DO K=NOLDM1,2,-1
4005	Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
4006	END DO
4007
4008	!-----------------------------------------------------------------------
4009	second_loop: DO K1=1,NNEW+1
4010	XK=XNEW(K1)
4011	third_loop: DO K2=2,NOLD
4012	IF(XOLD(K2)>XK)THEN
4013	KOLD=K2-1
4014	exit third_loop
4015	ENDIF
4016	END DO third_loop
4017
4018	IF (XOLD(K2) .le. XK) THEN
4019	YNEW(K1)=YOLD(NOLD)
4020	CYCLE second_loop
4021	ENDIF
4022
4023	IF (K1 .eq. 1 .or. K .ne. KOLD) THEN
4024	K=KOLD
4025	Y2K=Y2(K)
4026	Y2KP1=Y2(K+1)
4027	DX=XOLD(K+1)-XOLD(K)
4028	RDX=1./DX
4029	AK=.1666667RDX(Y2KP1-Y2K)
4030	BK=0.5*Y2K
4031	CK=RDX(YOLD(K+1)-YOLD(K))-.1666667DX*(Y2KP1+Y2K+Y2K)
4032	ENDIF
4033
4034	X=XK-XOLD(K)
4035	XSQ=X*X
4036	YNEW(K1)=AKXSQX+BKXSQ+CKX+YOLD(K)
4037	END DO second_loop
4038
4039	END SUBROUTINE SPLINE
4040
4041	!--------------------------------------------------------------------
4042	SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,&
4043	NSOIL,ISLTPK, &
4044	SM,SICE,STC,SMC,SH2O)
4045
4046	!! INTEGER, PARAMETER:: NSOTYP=9
4047	! INTEGER, PARAMETER:: NSOTYP=16
4048	INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE???
4049
4050	REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP)
4051	REAL :: STC(IMS:IME,NSOIL,JMS:JME), &
4052	SMC(IMS:IME,NSOIL,JMS:JME)
4053	REAL :: SH2O(IMS:IME,NSOIL,JMS:JME),SICE(IMS:IME,JMS:JME),&
4054	SM(IMS:IME,JMS:JME)
4055	REAL :: HLICE,GRAV,T0,BLIM
4056	INTEGER :: ISLTPK(IMS:IME,JMS:JME)
4057	CHARACTER(LEN=255) :: message
4058
4059	! Constants used in cold start SH2O initialization
4060	DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/
4061	DATA BLIM/5.5/
4062	! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
4063	! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
4064	! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, &
4065	! 0.465,0.404,0.439,0.421/
4066
4067
4068	!!! NOT SURE...PSIS=SATPSI, BETA=BB??
4069
4070	DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, &
4071	0.135, 0.617, 0.263, 0.098, 0.324, 0.468, &
4072	0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 /
4073
4074	DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, &
4075	6.66, 8.72, 8.17, 10.73, 10.39, 11.55, &
4076	5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 /
4077
4078	DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
4079	0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
4080	0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/
4081
4082	DO K=1,NSOIL
4083	DO J=JSTART,JM
4084	DO I=ISTART,IM
4085
4086	!tst
4087	IF (SMC(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
4088	if (K .eq. 1) then
4089	write(message,*) 'I,J,reducing SMC from ' ,I,J,SMC(I,K,J), 'to ', SMCMAX(ISLTPK(I,J))
4090	CALL wrf_debug(100,message)
4091	endif
4092	SMC(I,K,J)=SMCMAX(ISLTPK(I,J))
4093	ENDIF
4094	!tst
4095
4096	IF ( (SM(I,J) .lt. 0.5) .and. (SICE(I,J) .lt. 0.5) ) THEN
4097
4098	IF (ISLTPK(I,J) .gt. 19) THEN
4099	WRITE(message,*) 'FORCING ISLTPK at : ', I,J
4100	CALL wrf_message(message)
4101	ISLTPK(I,J)=9
4102	ELSEIF (ISLTPK(I,J) .le. 0) then
4103	WRITE(message,*) 'FORCING ISLTPK at : ', I,J
4104	CALL wrf_message(message)
4105	ISLTPK(I,J)=1
4106	ENDIF
4107
4108
4109	! cold start: determine liquid soil water content (SH2O)
4110	! SH2O <= SMC for T < 273.149K (-0.001C)
4111
4112	IF (STC(I,K,J) .LT. 273.149) THEN
4113
4114	! first guess following explicit solution for Flerchinger Eqn from Koren
4115	! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
4116
4117	BX = BETA(ISLTPK(I,J))
4118	IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM
4119
4120	if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then
4121	write(message,*) 'TROUBLE'
4122	CALL wrf_message(message)
4123	write(message,*) 'I,J: ', i,J
4124	CALL wrf_message(message)
4125	write(message,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),&
4126	psis(isltpk(I,J))
4127	CALL wrf_message(message)
4128	endif
4129
4130	if (BX .eq. 0 .or. STC(I,K,J) .eq. 0) then
4131	write(message,*) 'TROUBLE -- I,J,BX, STC: ', I,J,BX,STC(I,K,J)
4132	CALL wrf_message(message)
4133	endif
4134	FK = (((HLICE/(GRAV(-PSIS(ISLTPK(I,J))))) &
4135	((STC(I,K,J)-T0)/STC(I,K,J)))** &
4136	(-1/BX))*SMCMAX(ISLTPK(I,J))
4137	IF (FK .LT. 0.02) FK = 0.02
4138	SH2O(I,K,J) = MIN ( FK, SMC(I,K,J) )
4139	! ----------------------------------------------------------------------
4140	! now use iterative solution for liquid soil water content using
4141	! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
4142	! initial guess for SH2O from above explicit first guess.
4143
4144	SH2O(I,K,J)=FRH2O_init(STC(I,K,J),SMC(I,K,J),SH2O(I,K,J), &
4145	SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), &
4146	PSIS(ISLTPK(I,J)))
4147
4148	ELSE ! above freezing
4149	SH2O(I,K,J)=SMC(I,K,J)
4150	ENDIF
4151
4152
4153	ELSE ! water point
4154	SH2O(I,K,J)=SMC(I,K,J)
4155
4156	ENDIF ! test on land/ice/sea
4157	if (SH2O(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
4158	write(message,*) 'SH2O > THAN SMCMAX ', I,J,SH2O(I,K,J),SMCMAX(ISLTPK(I,J)),SMC(I,K,J)
4159	CALL wrf_message(message)
4160	endif
4161
4162	ENDDO
4163	ENDDO
4164	ENDDO
4165
4166	END SUBROUTINE NMM_SH2O
4167
4168	!-------------------------------------------------------------------
4169
4170	FUNCTION FRH2O_init(TKELV,SMC,SH2O,SMCMAX,B,PSIS)
4171
4172	IMPLICIT NONE
4173
4174	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4175	! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
4176	! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION
4177	! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
4178	! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
4179	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4180	!
4181	! New version (JUNE 2001): much faster and more accurate newton iteration
4182	! achieved by first taking log of eqn cited above -- less than 4
4183	! (typically 1 or 2) iterations achieves convergence. Also, explicit
4184	! 1-step solution option for special case of parameter Ck=0, which reduces
4185	! the original implicit equation to a simpler explicit form, known as the
4186	! ""Flerchinger Eqn". Improved handling of solution in the limit of
4187	! freezing point temperature T0.
4188	!
4189	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4190	!
4191	! INPUT:
4192	!
4193	! TKELV.........Temperature (Kelvin)
4194	! SMC...........Total soil moisture content (volumetric)
4195	! SH2O..........Liquid soil moisture content (volumetric)
4196	! SMCMAX........Saturation soil moisture content (from REDPRM)
4197	! B.............Soil type "B" parameter (from REDPRM)
4198	! PSIS..........Saturated soil matric potential (from REDPRM)
4199	!
4200	! OUTPUT:
4201	! FRH2O.........supercooled liquid water content.
4202	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4203
4204	REAL B
4205	REAL BLIM
4206	REAL BX
4207	REAL CK
4208	REAL DENOM
4209	REAL DF
4210	REAL DH2O
4211	REAL DICE
4212	REAL DSWL
4213	REAL ERROR
4214	REAL FK
4215	REAL FRH2O_init
4216	REAL GS
4217	REAL HLICE
4218	REAL PSIS
4219	REAL SH2O
4220	REAL SMC
4221	REAL SMCMAX
4222	REAL SWL
4223	REAL SWLK
4224	REAL TKELV
4225	REAL T0
4226
4227	INTEGER NLOG
4228	INTEGER KCOUNT
4229	PARAMETER (CK=8.0)
4230	! PARAMETER (CK=0.0)
4231	PARAMETER (BLIM=5.5)
4232	! PARAMETER (BLIM=7.0)
4233	PARAMETER (ERROR=0.005)
4234
4235	PARAMETER (HLICE=3.335E5)
4236	PARAMETER (GS = 9.81)
4237	PARAMETER (DICE=920.0)
4238	PARAMETER (DH2O=1000.0)
4239	PARAMETER (T0=273.15)
4240
4241	! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) ####
4242	! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT ####
4243	! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES ####
4244	! ################################################################
4245	!
4246	BX = B
4247	IF ( B .GT. BLIM ) BX = BLIM
4248	! ------------------------------------------------------------------
4249
4250	! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
4251	NLOG=0
4252	KCOUNT=0
4253
4254	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4255	! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
4256	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4257
4258
4259	IF (TKELV .GT. (T0 - 1.E-3)) THEN
4260
4261	FRH2O_init=SMC
4262
4263	ELSE
4264
4265	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4266	IF (CK .NE. 0.0) THEN
4267
4268	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4269	! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
4270	! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
4271
4272	! INITIAL GUESS FOR SWL (frozen content)
4273	SWL = SMC-SH2O
4274	! KEEP WITHIN BOUNDS.
4275	IF (SWL .GT. (SMC-0.02)) SWL=SMC-0.02
4276	IF(SWL .LT. 0.) SWL=0.
4277	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4278	! C START OF ITERATIONS
4279	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4280	DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
4281	NLOG = NLOG+1
4282	DF = ALOG(( PSISGS/HLICE ) ( ( 1.+CKSWL )2. ) &
4283	( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
4284	DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
4285	SWLK = SWL - DF/DENOM
4286	! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
4287	IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
4288	IF(SWLK .LT. 0.) SWLK = 0.
4289	! MATHEMATICAL SOLUTION BOUNDS APPLIED.
4290	DSWL=ABS(SWLK-SWL)
4291	SWL=SWLK
4292	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4293	! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
4294	! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
4295	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4296	IF ( DSWL .LE. ERROR ) THEN
4297	KCOUNT=KCOUNT+1
4298	END IF
4299	END DO
4300	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4301	! C END OF ITERATIONS
4302	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4303	! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
4304	FRH2O_init = SMC - SWL
4305
4306	! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC
4307
4308	ENDIF
4309
4310	IF (KCOUNT .EQ. 0) THEN
4311	! Print*,'Flerchinger used in NEW version. Iterations=',NLOG
4312
4313	! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4314	! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
4315	! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC
4316
4317	FK=(((HLICE/(GS(-PSIS)))((TKELV-T0)/TKELV))*(-1/BX))SMCMAX
4318	! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
4319	IF (FK .LT. 0.02) FK = 0.02
4320	FRH2O_init = MIN ( FK, SMC )
4321
4322	! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC
4323
4324	ENDIF
4325
4326	ENDIF
4327
4328	RETURN
4329
4330	END FUNCTION FRH2O_init
4331
4332
4333	!--------------------------------------------------------------------
4334
4335	SUBROUTINE init_module_initialize
4336	END SUBROUTINE init_module_initialize
4337
4338	!---------------------------------------------------------------------
4339
4340	END MODULE module_initialize_real

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

Download in other formats:

Original Format