source: lmdz_wrf/trunk/WRFV3/dyn_nmm/module_initialize_real.F @ 1465

Last change on this file since 1465 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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