wrfv2_fire/dyn_nmm/module_initialize_real.F

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