174,177d173
<     REAL, ALLOCATABLE           :: rainwat_in(:,:,:)            
<     REAL, ALLOCATABLE           :: snowmr_in(:,:,:)            
<     REAL, ALLOCATABLE           :: cloudice_in(:,:,:)            
<     REAL, ALLOCATABLE           :: cloudwat_in(:,:,:)            
182d177
<     REAL, ALLOCATABLE           :: q1d_out(:)
194c189
<     REAL, ALLOCATABLE,DIMENSION(:)  :: uin,vin,qin,dfl
---
>     REAL, ALLOCATABLE,DIMENSION(:)  :: uin,vin,dfl
210c205
<     REAL           :: qest,tnew,pdo_0d
---
>     REAL           :: qest,tnew
223c218
<     ALLOCATE (z3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim ))
---
>     ALLOCATE(z3d_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim ))
226,229d220
<     ALLOCATE (rainwat_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim))
<     ALLOCATE (snowmr_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim))
<     ALLOCATE (cloudice_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim))
<     ALLOCATE (cloudwat_in(dom_meta%xdim,dom_meta%ydim,dom_meta%zdim))
258,259d248
< 
< 
274c263,265
< 
---
> 	write(6,*) 'z3d(1,1,1), z3d(1,1,2): ', z3d_in(1,1,1),&
> 		z3d_in(1,1,2)
>       
291,358d281
< 
< !!!! try not getting these variables
< 
< !	goto 979
< 
<     IF(.not.Dryrun)PRINT '(A,A)', & 
< 		'INTERP_STATE_NMMH: Getting 3D rain water field...'
<     CALL get_variable(setup_info%input_prefix, rainwat_id, dom_meta%id, &
<       setup_info%current_date,status)
<     IF (status.EQ.0) THEN
<       rainwat_in = real_array(:,:,:,1,1)
<       num_processed = num_processed + 1
<       processed_var_list(num_processed) = rainwat_id
<     ELSE
<       PRINT *,'Problem getting ', rainwat_id
<       STOP 'No_RAINWAT3D_data'
<     ENDIF
<     if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(rainwat_in), &
<       MAXVAL(rainwat_in)
< 
< 
<     IF(.not.Dryrun)PRINT '(A,A)', & 
< 		'INTERP_STATE_NMMH: Getting 3D snow mixing ratio field...'
<     CALL get_variable(setup_info%input_prefix, snowmr_id, dom_meta%id, &
<       setup_info%current_date,status)
<     IF (status.EQ.0) THEN
<       snowmr_in = real_array(:,:,:,1,1)
<       num_processed = num_processed + 1
<       processed_var_list(num_processed) = snowmr_id
<     ELSE
<       PRINT *,'Problem getting ', snowmr_id
<       STOP 'No_SNOWMR_data'
<     ENDIF
<     if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(snowmr_in), &
<       MAXVAL(snowmr_in)
< 
<     IF(.not.Dryrun)PRINT '(A,A)', & 
< 		'INTERP_STATE_NMMH: Getting 3D cloud ice field...'
<     CALL get_variable(setup_info%input_prefix, cloudice_id, dom_meta%id, &
<       setup_info%current_date,status)
<     IF (status.EQ.0) THEN
<       cloudice_in = real_array(:,:,:,1,1)
<       num_processed = num_processed + 1
<       processed_var_list(num_processed) = cloudice_id
<     ELSE
<       PRINT *,'Problem getting ', cloudice_id
<       STOP 'No_CLOUDICE_data'
<     ENDIF
<     if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(cloudice_in), &
<       MAXVAL(cloudice_in)
< 
<     IF(.not.Dryrun)PRINT '(A,A)', & 
< 		'INTERP_STATE_NMMH: Getting 3D cloud water field...'
<     CALL get_variable(setup_info%input_prefix, cloudwat_id, dom_meta%id, &
<       setup_info%current_date,status)
<     IF (status.EQ.0) THEN
<       cloudwat_in = real_array(:,:,:,1,1)
<       num_processed = num_processed + 1
<       processed_var_list(num_processed) = cloudwat_id
<     ELSE
<       PRINT *,'Problem getting ', cloudwat_id
<       STOP 'No_CLOUDWAT_data'
<     ENDIF
<     if(.not.DryRun)print '(A,2F9.1)','Min/Max values = ',MINVAL(cloudwat_in), &
<       MAXVAL(cloudwat_in)
< 
<   979	continue
< 
397,399d319
< 
< 	write(6,*) 'press3d_in(1,1,1): ', press3d_in(1,1,1)
< 
409,413d328
< 	write(6,*) 'at 1,1:'
< 	do L=1,dom_meta%zdim
< 	write(6,*) 'L,p3d,z3d,t3d: ', L, &
<                press3d_in(1,1,L),z3d_in(1,1,L),t3d_in(1,1,L)
< 	enddo
415,432d329
< 	write(6,*) 'at 50,50:'
< 	do L=1,dom_meta%zdim
< 	write(6,*) 'L,p3d,z3d,t3d: ', L, &
<                press3d_in(50,50,L),z3d_in(50,50,L),t3d_in(50,50,L)
< 	enddo
< 
< 	write(6,*) 'at 100,100:'
< 	do L=1,dom_meta%zdim
< 	write(6,*) 'L,p3d,z3d,t3d: ', L, &
<                press3d_in(100,100,L),z3d_in(100,100,L),t3d_in(100,100,L)
< 	enddo
< 
< 	write(6,*) 'at 200,250:'
< 	do L=1,dom_meta%zdim
< 	write(6,*) 'L,p3d,z3d,t3d: ', L, &
<                press3d_in(200,250,L),z3d_in(200,250,L),t3d_in(200,250,L)
< 	enddo
< 
440a338
> 
450,453d347
<              rainwat_in(I,J,L)=rainwat_in(I,J,L+1)
<              snowmr_in(I,J,L)=snowmr_in(I,J,L+1)
<              cloudice_in(I,J,L)=cloudice_in(I,J,L+1)
<              cloudwat_in(I,J,L)=cloudwat_in(I,J,L+1)
458,479d351
< 	elseif(press3d_in(1,1,1) .eq. 100000.) then
< 	elim_lowest=.true.
< 	LMAX=dom_meta%zdim-1
< 	write(6,*) 'in unexpected branch of elimination'
< 	write(6,*) 'reduced LMAX to: ', LMAX
< 
< 	do L=1,LMAX
< 	 do J=1,dom_meta%ydim
< 	  do I=1,dom_meta%xdim
<              z3d_in(I,J,L)=z3d_in(I,J,L+1)
<              t3d_in(I,J,L)=t3d_in(I,J,L+1)
<              rh3d_in(I,J,L)=rh3d_in(I,J,L+1)
<              rainwat_in(I,J,L)=rainwat_in(I,J,L+1)
<              snowmr_in(I,J,L)=snowmr_in(I,J,L+1)
<              cloudice_in(I,J,L)=cloudice_in(I,J,L+1)
<              cloudwat_in(I,J,L)=cloudwat_in(I,J,L+1)
< 	  enddo
< 	 enddo
< 	enddo
< 
< 
< 
482,487d353
< 	write(6,*) 'at 50,50, post elim:'
< 	do L=1,LMAX
< 	write(6,*) 'L,p3d,z3d,t3d: ', L, &
<                press3d_in(50,50,L),z3d_in(50,50,L),t3d_in(50,50,L)
< 	enddo
< 
501d366
< 
504a370,432
> !    IF (MAXVAL(press3d_in(:,:,1)).LT.200000.) THEN
> !      use_sfc = .false.
> !      skip_lowest = .false.
> !    ELSE
> 
> !        CALL shift_sfc_slab(z3d_in)
>         ! Here is the sanity check...
> !        adjustz_cnt = 0
> !        DO j = 1,dom_meta%ydim
> !          DO i = 1,dom_meta%xdim 
> !            k = sfc_k_ind(i,j)
> !            IF ( k .EQ. 1)  THEN
> !              IF ( z3d_in(i,j,k).GT.z3d_in(i,j,k+1) ) THEN
> !                ! PRINT *, 'Height problem at i/j/k/k+1',i,j,k,k+1, 'Values =',&
> !                !       z3d_in(i,j,k), z3d_in(i,j,k+1)
> !                ! If we are within 10 mb, do a simple correction, otherwise
> !                ! stop with an error
> !                deltap = press3d_in(i,j,k) - press3d_in(i,j,k+1)
> !                IF ((deltap .LT. 1000.).AND.(deltap.GE.0)) THEN
> !                  ! PRINT *, 'Applying simple correction.'
> !                  adjustz_cnt = adjustz_cnt + 1
> !                  z3d_in(i,j,k) = z3d_in(i,j,k+1) - deltap/10.
> !                ELSE
> !                  PRINT *, 'Cannot correct.  Need to debug...'
> !                  STOP 'bad_surface_height_vs_press'
> !                ENDIF
> !              ENDIF
> !            ELSE
> !              IF ( (z3d_in(i,j,k).LT.z3d_in(i,j,k-1)) .OR. &
> !                   (z3d_in(i,j,k).GT.z3d_in(i,j,k+1)) ) THEN
> !                adjustz_cnt = adjustz_cnt+1
> !                ! PRINT *, 'Bad height at i/j/k = ', i,j,k
> !                ! PRINT *, 'Z0/Z/Z1 = ', z3d_in(i,j,k-1:k+1) 
> ! PRINT *, 'Replacing with interpolated value'
> !                wgt0 = 1.0 - ( ALOG( press3d_in(i,j,k-1)/press3d_in(i,j,k))/&
> !                               ALOG( press3d_in(i,j,k-1)/press3d_in(i,j,k+1)))  
> !                z3d_in(i,j,k) = wgt0*z3d_in(i,j,k-1) + &
> !                                (1.-wgt0)*z3d_in(i,j,k+1)     
> !                ! PRINT *, 'New Z0/Z/Z1 = ',z3d_in(i,j,k-1:k+1)
> !              ENDIF
> !            ENDIF
> !          ENDDO
> !        ENDDO
> ! IF(.not.Dryrun)PRINT *, 'Number of columns where Zsfc had to be adjusted: ', &
> !                  adjustz_cnt
>     ! ELSE
>         ! Sort by height and derive surface pressure 
>         ! COMPLETE THIS WORK LATER.  FOR NOW, JUST MENTION
>         ! NOT SUPPORTED AND SET USE_SFC = .FALSE. 
>     !   PRINT *, 'Sorting by ZSFC and deriving PSFC'
>     !   use_sfc = .true.
>     !   skip_lowest = .false.
> 
> !      ENDIF
> 
> !    ENDIF 
> 
>     ! Generate the 3d array of Eta values on the background gridpoints.
>     ! Subroutine also returns the 2D mu array, which is computed using
>     ! the WRF topography to ensure a pressure field consistent with the
>     ! forecast model.
> 
> 
697c625
< 	write(6,*) 'H points'
---
> !	write(6,*) 'H points'
699c627
< 	write(6,633) (terrain_hgt_t(I,J),I=1,dom_meta%xdim,dom_meta%xdim/10)
---
> !	write(6,633) (terrain_hgt_t(I,J),I=1,dom_meta%xdim,dom_meta%xdim/10)
702,711d629
< 	write(6,*) 'z3d_in(1) points'
< 	do J=dom_meta%ydim,1,-dom_meta%ydim/20
< 	write(6,633) (z3d_in(I,J,1),I=1,dom_meta%xdim,dom_meta%xdim/10)
< 	enddo
< 
< 	write(6,*) 'press3d_in(1) points'
< 	do J=dom_meta%ydim,1,-dom_meta%ydim/20
< 	write(6,635) (press3d_in(I,J,1),I=1,dom_meta%xdim,dom_meta%xdim/10)
< 	enddo
< 
713d630
<   635	format(20(f7.0,1x))
735a653,661
> !
> !	NOTE:  The above ln P formulation and below P formulation
> !		agree to within less than 1 mb.  Typically within
> !		.01-.1 mb
> 
> !	dum2ds(I,J)= (press3d_in(i,j,2)-press3d_in(i,j,1))/ &
> !		(z3d_in(i,j,2)-z3d_in(i,j,1))* &
> !		(terrain_hgt_t(i,j)-z3d_in(i,j,1))+press3d_in(i,j,1)
> 
757c683
< 	write(6,*) 'must have flukey situation ', I,J
---
> 	write(6,*) 'must have flukey situation'
787d712
< 	allocate(qin(LMAX))
1498d1422
<     ALLOCATE (q1d_out(dom_out%zdim))
1512c1436,1438
< 	PDO_0D=PDO(I,J)
---
>        DO LD=1,LMAX-1 !LMIN-1
>         ALPR(LD)=ALOG(press3d_in(I,J,LD)/press3d_in(I,J,LD+1))
>        ENDDO
1514,1517c1440,1442
< 	do L=1,LMAX
< 	pin(L)=press3d_in(I,J,L)
< 	qin(L)=qv3d_in(I,J,L)
< 	enddo
---
>        DO LD=1,LMAX !LMIN
>         ALP(LD)=ALOG(press3d_in(I,J,LD))
>        ENDDO
1519,1520c1444
< 	CALL INTERP_Q_LNP(pin(1:LMAX),qin,LMAX,SGML1,SGML2,&
<                           q1d_out,LKO,PDO_0D,PDTOP,PT)
---
> !	integrating from TOA downward...
1522c1446,1448
< 	DO L=1,LKO
---
>         DO L=LKO,1,-1
> 	  ALPOUT=log(pt+sgml1(l)*pdtop+sgml2(l)*pdo(i,j))
> 	  PMIDO=pt+sgml1(l)*pdtop+sgml2(l)*pdo(i,j)
1524c1450
< 	q3d_out(I,J,L)=q1d_out(L)
---
> 	Q3D_OUT(I,J,L)=-9999.
1525a1452,1531
> 	IF (PMIDO .gt. press3d_in(I,J,1)) THEN  ! extrap downward
> 
> 	DLNQDLNP= (log(qv3d_in(i,j,1))-log(qv3d_in(i,j,2))) / &
> 		(log(press3d_in(I,J,1))-log(press3d_in(I,J,2)))
> 
> 	q3d_out(I,J,L)= exp( &
> 			log(qv3d_in(i,j,1)) + DLNQDLNP * &
> 			(log(pmido)-log(press3d_in(i,j,1))) &
> 			)
> 
> 	   if (I .eq. dom_meta%xdim/2 .and. J .eq. dom_meta%ydim/2) then
>            write(6,*) 'extrap down Q3DOUT, I,J,L,q3d_out: ',  &
>                                 I,J,L,q3d_out(I,J,L)
>            write(6,*) 'qv3d_in,DLNQDLNP,pmido,press3d_in: ', &
>                            qv3d_in(i,j,1),DLNQDLNP,pmido,press3d_in(i,j,1)
> 	   endif
> 
> 	   if (q3d_out(I,J,L) .gt. .05 ) then ! 50 g/kg would be extreme, yes?
> 		write(6,*) 'CRAZY Q3DOUT, I,J,L,q3d_out: ',  &
> 				I,J,L,q3d_out(I,J,L)	
> 		write(6,*) 'qv3d_in,DLNQDLNP,pmido,press3d_in: ', &
>                            qv3d_in(i,j,1),DLNQDLNP,pmido,press3d_in(i,j,1)
> 	   endif
> 
> 	ELSEIF (pmido .lt. press3d_in(I,J,LMAX)) THEN ! extrap up
> 
> 	   DLNQDLNP= (log(qv3d_in(i,j,LMAX-1))- &
> 			log(qv3d_in(i,j,LMAX))) / &
> 		(log(press3d_in(i,j,LMAX-1))- &
> 			log(press3d_in(i,j,LMAX)))
> 
> 	   q3d_out(I,J,L)= exp( &
> 			log(qv3d_in(i,j,LMAX)) + DLNQDLNP * &
> 			(log(pmido)-log(press3d_in(i,j,LMAX))) &
> 			)
> 
> 	   if (q3d_out(I,J,L) .gt. .05 ) then ! 50 g/kg would be extreme, yes?
> 		write(6,*) 'CRAZY Q3DOUT, I,J,L,q3d_out: ',  &
> 				I,J,L,q3d_out(I,J,L)	
> 	   endif
> 
> 
> 	ELSE ! normal case
> 
> 	DO LD=LMAX-1,1,-1
> 
> 	if (pmido .le. press3d_in(I,J,LD) .and.  &
>             pmido .ge. press3d_in(I,J,LD+1)) then
> 
> !	  QLD=qv3d_in(I,J,LD)
>         QLD=AMAX1( qv3d_in(I,J,LD+1),1.e-12 )  ! to avoid computing log(0.)
> 	qv3d_in(I,J,LD)=amax1(qv3d_in(I,J,LD), 1.e-12)
> 
> !old	       CF=(ALP(LD)-ALPOUT)/ALPR(LD+1)
> 	       CF=(ALPOUT-ALP(LD+1))/ALPR(LD)
> 
> !old        q3d_out(I,J,L)=EXP(ALOG(QLD)+(ALOG(qv3d_in(I,J,LD+1))-ALOG(QLD))*CF)
>         q3d_out(I,J,L)=EXP(ALOG(QLD)+(ALOG(qv3d_in(I,J,LD))-ALOG(QLD))*CF)
> 
> 
> !	if (I .eq. dom_meta%xdim/2 .and. J .eq. dom_meta%ydim/2) then
> !	 write(6,*) 'LD, QLD, qv3d_in(LD+1): ', L, QLD,qv3d_in(I,J,LD+1)
> !	write(6,*) 'CF, q3d_out(I,J,L): ', L, CF, q3d_out(I,J,L)
> !	endif
> 
> 	   if (q3d_out(I,J,L) .gt. .05 ) then ! 50 g/kg would be extreme, yes?
> 		write(6,*) 'CRAZY Q3DOUT, I,J,L,q3d_out: ',  &
> 				I,J,L,q3d_out(I,J,L)	
> 		write(6,*) 'qld,qv3d(LD+1), cf: ', &
> 			QLD,qv3d_in(I,J,LD+1),CF
> 		write(6,*) 'ALP,ALPOUT,ALPR(LD+1): ', ALP(LD),ALPOUT, &
> 			ALPR(LD+1)
> 	   endif
> 	
> 	endif  ! check of pressure within range
> 
> 	ENDDO !LD
> 
> 	ENDIF
> 
1567c1573
<     DEALLOCATE(press3d_in,qin)
---
>     DEALLOCATE(press3d_in)
1576,1577c1582
< 	deallocate(rainwat_in, snowmr_in, cloudice_in, cloudwat_in)
< 	deallocate(t3d_out,q3d_out,q1d_out,u3d_out,v3d_out,pdo)
---
> 	deallocate(t3d_out,q3d_out,u3d_out,v3d_out,pdo)
2320,2387d2324
< 
< !-------------------------------------------------------------------------------
< 
< 	SUBROUTINE INTERP_Q_LNP(press1d_in,q1d_in,LMAX,SGML1,SGML2,Q1D_OUT,LM,&
<                                 PDO,PDTOP,PT )
< 
<         INTEGER:: LMAX,LD,LM,L
<         REAL:: press1d_in(LMAX),q1d_in(LMAX)
<         REAL:: ALPR(LMAX),ALP(LMAX),QLD
<         REAL:: SGML1(LM),SGML2(LM),PDO,PDTOP,PT
<         REAL:: Q1D_OUT(LM),DLNQDLNP,PMIDO,ALPOUT,CF
< 
<        DO LD=1,LMAX-1 
<         ALPR(LD)=ALOG(press1d_in(LD)/press1d_in(LD+1))
<        ENDDO
< 
<        DO LD=1,LMAX 
<         ALP(LD)=ALOG(press1d_in(LD))
<        ENDDO
< 
< !	integrating from TOA downward...
< 
<         DO L=LM,1,-1
< 
< 	ALPOUT=log(pt+sgml1(l)*pdtop+sgml2(l)*pdo)
< 	PMIDO=pt+sgml1(l)*pdtop+sgml2(l)*pdo
< 	Q1D_OUT(L)=-9999.
< 
< 	IF (PMIDO .gt. press1d_in(1)) THEN  ! extrap downward
< 
< 	DLNQDLNP= (log(q1d_in(1))-log(q1d_in(2))) / &
< 	  	  (log(press1d_in(1))-log(press1d_in(2)))
< 
< 	q1d_out(L)= exp(log(q1d_in(1)) + DLNQDLNP * &
< 			(log(pmido)-log(press1d_in(1))) )
< 
< 	ELSEIF (pmido .lt. press1d_in(LMAX)) THEN ! extrap up
< 
< 	   DLNQDLNP= (log(q1d_in(LMAX-1))- &
< 			log(q1d_in(LMAX))) / &
< 		(log(press1d_in(LMAX-1))- &
< 			log(press1d_in(LMAX)))
< 
< 	   q1d_out(L)= exp( &
< 			log(q1d_in(LMAX)) + DLNQDLNP * &
< 			(log(pmido)-log(press1d_in(LMAX))) &
< 			)
< 
< 	ELSE ! normal case
< 
< 	DO LD=LMAX-1,1,-1
< 
< 	IF (pmido .le. press1d_in(LD) .and.  &
<             pmido .ge. press1d_in(LD+1)) THEN
< 
<         QLD=AMAX1( q1d_in(LD+1),1.e-12 )  ! to avoid computing log(0.)
< 	q1d_in(LD)=amax1(q1d_in(LD), 1.e-12)
<         CF=(ALPOUT-ALP(LD+1))/ALPR(LD)
<         q1d_out(L)=EXP(ALOG(QLD)+(ALOG(q1d_in(LD))-ALOG(QLD))*CF)
< 	
< 	ENDIF  ! check of pressure within range
< 
< 	ENDDO !LD
<         ENDIF
<         ENDDO !LM
<  
< 	END SUBROUTINE INTERP_Q_LNP
<