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