wrfv2_fire/share/module_soil_pre.F

   1 #if ( ! NMM_CORE == 1 )
   2
   3 MODULE module_soil_pre
   4
   5    USE module_date_time
   6    USE module_state_description
   7
   8 CONTAINS
   9
  10    SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
  11                                       xland , landusef , isltyp , soilcat , soilctop , &
  12                                       soilcbot , tmn , &
  13                                       seaice_threshold , &
  14                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  15                                       iswater , isice , &
  16                                       scheme , &
  17                                       ids , ide , jds , jde , kds , kde , &
  18                                       ims , ime , jms , jme , kms , kme , &
  19                                       its , ite , jts , jte , kts , kte )
  20
  21       IMPLICIT NONE
  22
  23       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  24                               ims , ime , jms , jme , kms , kme , &
  25                               its , ite , jts , jte , kts , kte , &
  26                               iswater , isice
  27       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
  28
  29       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
  30       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
  31       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
  32       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  33       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
  34                                                            vegcat, xland , soilcat , tmn
  35       REAL , INTENT(IN) :: seaice_threshold
  36
  37       INTEGER :: i , j , num_seaice_changes , loop
  38       CHARACTER (LEN=132) :: message
  39
  40       num_seaice_changes = 0
  41       fix_seaice : SELECT CASE ( scheme )
  42
  43          CASE ( SLABSCHEME )
  44             DO j = jts , MIN(jde-1,jte)
  45                DO i = its , MIN(ide-1,ite)
  46                   IF ( xice(i,j) .GT. 200.0 ) THEN
  47                      xice(i,j) = 0.
  48                      num_seaice_changes = num_seaice_changes + 1
  49                   END IF
  50                END DO
  51             END DO
  52             IF ( num_seaice_changes .GT. 0 ) THEN
  53                WRITE ( message , FMT='(A,I6)' ) &
  54                'Total pre number of sea ice locations removed (due to FLAG values) = ', &
  55                num_seaice_changes
  56                CALL wrf_debug ( 0 , message )
  57             END IF
  58             num_seaice_changes = 0
  59             DO j = jts , MIN(jde-1,jte)
  60                DO i = its , MIN(ide-1,ite)
  61                   IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
  62                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  63                      xice(i,j) = 1.
  64                      num_seaice_changes = num_seaice_changes + 1
  65                      if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
  66                      vegcat(i,j)=isice
  67                      ivgtyp(i,j)=isice
  68                      lu_index(i,j)=isice
  69                      landmask(i,j)=1.
  70                      xland(i,j)=1.
  71                      DO loop=1,num_veg_cat
  72                         landusef(i,loop,j)=0.
  73                      END DO
  74                      landusef(i,ivgtyp(i,j),j)=1.
  75
  76                      isltyp(i,j) = 16
  77                      soilcat(i,j)=isltyp(i,j)
  78                      DO loop=1,num_soil_top_cat
  79                         soilctop(i,loop,j)=0
  80                      END DO
  81                      DO loop=1,num_soil_bot_cat
  82                         soilcbot(i,loop,j)=0
  83                      END DO
  84                      soilctop(i,isltyp(i,j),j)=1.
  85                      soilcbot(i,isltyp(i,j),j)=1.
  86                   END IF
  87                END DO
  88             END DO
  89             IF ( num_seaice_changes .GT. 0 ) THEN
  90                WRITE ( message , FMT='(A,I6)' ) &
  91                'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
  92                CALL wrf_debug ( 0 , message )
  93             END IF
  94
  95          CASE ( LSMSCHEME , RUCLSMSCHEME )
  96             num_seaice_changes = 0
  97             DO j = jts , MIN(jde-1,jte)
  98                DO i = its , MIN(ide-1,ite)
  99                   IF ( landmask(i,j) .GT. 0.5 ) THEN
 100                      if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
 101                      xice(i,j) = 0.
 102                   END IF
 103                END DO
 104             END DO
 105             IF ( num_seaice_changes .GT. 0 ) THEN
 106                WRITE ( message , FMT='(A,I6)' ) &
 107                'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
 108                CALL wrf_debug ( 0 , message )
 109             END IF
 110
 111       END SELECT fix_seaice
 112
 113    END SUBROUTINE adjust_for_seaice_pre
 114
 115    SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
 116                                       xland , landusef , isltyp , soilcat , soilctop , &
 117                                       soilcbot , tmn , vegfra , &
 118                                       tslb , smois , sh2o , &
 119                                       seaice_threshold , &
 120                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
 121                                       num_soil_layers , &
 122                                       iswater , isice , &
 123                                       scheme , &
 124                                       ids , ide , jds , jde , kds , kde , &
 125                                       ims , ime , jms , jme , kms , kme , &
 126                                       its , ite , jts , jte , kts , kte )
 127
 128       IMPLICIT NONE
 129
 130       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
 131                               ims , ime , jms , jme , kms , kme , &
 132                               its , ite , jts , jte , kts , kte , &
 133                               iswater , isice
 134       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
 135       INTEGER , INTENT(IN) :: num_soil_layers
 136
 137       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
 138       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
 139       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
 140       REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
 141       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
 142       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
 143                                                            vegcat, xland , soilcat , tmn , &
 144                                                            tsk_old , vegfra
 145       REAL , INTENT(IN) :: seaice_threshold
 146       REAL :: total_depth , mid_point_depth
 147
 148       INTEGER :: i , j , num_seaice_changes , loop
 149       CHARACTER (LEN=132) :: message
 150
 151       num_seaice_changes = 0
 152       fix_seaice : SELECT CASE ( scheme )
 153
 154          CASE ( SLABSCHEME )
 155
 156          CASE ( LSMSCHEME , RUCLSMSCHEME )
 157             DO j = jts , MIN(jde-1,jte)
 158                DO i = its , MIN(ide-1,ite)
 159                   IF ( xice(i,j) .GT. 200.0 ) THEN
 160                      xice(i,j) = 0.
 161                      num_seaice_changes = num_seaice_changes + 1
 162                   END IF
 163                END DO
 164             END DO
 165             IF ( num_seaice_changes .GT. 0 ) THEN
 166                WRITE ( message , FMT='(A,I6)' ) &
 167                'Total post number of sea ice locations removed (due to FLAG values) = ', &
 168                num_seaice_changes
 169                CALL wrf_debug ( 0 , message )
 170             END IF
 171             num_seaice_changes = 0
 172             DO j = jts , MIN(jde-1,jte)
 173                DO i = its , MIN(ide-1,ite)
 174                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
 175                        ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
 176                      tsk(i,j) = tsk_old(i,j)
 177                   END IF
 178                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
 179                        ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
 180                      print *,'TSK woes in seaice post, i,j=',i,j,'  tsk = ',tsk(i,j), tsk_old(i,j)
 181                      CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
 182                   ELSE IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
 183                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
 184                      xice(i,j) = 1.
 185                      num_seaice_changes = num_seaice_changes + 1
 186                      if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
 187                      vegcat(i,j)=isice
 188                      ivgtyp(i,j)=isice
 189                      lu_index(i,j)=isice
 190                      landmask(i,j)=1.
 191                      xland(i,j)=1.
 192                      vegfra(i,j)=0.
 193                      DO loop=1,num_veg_cat
 194                         landusef(i,loop,j)=0.
 195                      END DO
 196                      landusef(i,ivgtyp(i,j),j)=1.
 197
 198                      tsk_old(i,j) = tsk(i,j)
 199
 200                      isltyp(i,j) = 16
 201                      soilcat(i,j)=isltyp(i,j)
 202                      DO loop=1,num_soil_top_cat
 203                         soilctop(i,loop,j)=0
 204                      END DO
 205                      DO loop=1,num_soil_bot_cat
 206                         soilcbot(i,loop,j)=0
 207                      END DO
 208                      soilctop(i,isltyp(i,j),j)=1.
 209                      soilcbot(i,isltyp(i,j),j)=1.
 210
 211                      total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
 212                      DO loop = 1,num_soil_layers
 213                         mid_point_depth=(total_depth/num_soil_layers)/2. + &
 214                                         (loop-1)*(total_depth/num_soil_layers)
 215                         tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
 216                                             mid_point_depth*tmn(i,j) ) / total_depth
 217                      END DO
 218
 219                      DO loop=1,num_soil_layers
 220                         smois(i,loop,j) = 1.0
 221                         sh2o(i,loop,j)  = 0.0
 222                      END DO
 223                   ELSE IF ( xice(i,j) .LT. 0.5 ) THEN
 224                      xice(i,j) = 0.
 225                   END IF
 226                END DO
 227             END DO
 228             IF ( num_seaice_changes .GT. 0 ) THEN
 229                WRITE ( message , FMT='(A,I6)' ) &
 230                'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
 231                CALL wrf_debug ( 0 , message )
 232             END IF
 233
 234       END SELECT fix_seaice
 235
 236    END SUBROUTINE adjust_for_seaice_post
 237
 238    SUBROUTINE process_percent_cat_new ( landmask ,  &
 239                                 landuse_frac , soil_top_cat , soil_bot_cat , &
 240                                 isltyp , ivgtyp , &
 241                                 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
 242                                 ids , ide , jds , jde , kds , kde , &
 243                                 ims , ime , jms , jme , kms , kme , &
 244                                 its , ite , jts , jte , kts , kte , &
 245                                 iswater )
 246
 247       IMPLICIT NONE
 248
 249       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
 250                               ims , ime , jms , jme , kms , kme , &
 251                               its , ite , jts , jte , kts , kte , &
 252                               iswater
 253       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
 254       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
 255       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
 256       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
 257       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
 258       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
 259
 260       INTEGER :: i , j , l , ll, dominant_index
 261       REAL :: dominant_value
 262
 263 #ifdef WRF_CHEM
 264 !      REAL :: lwthresh = .99
 265       REAL :: lwthresh = .50
 266 #else
 267       REAL :: lwthresh = .50
 268 #endif
 269
 270       INTEGER , PARAMETER :: iswater_soil = 14
 271       INTEGER :: iforce
 272       CHARACTER (LEN=132) :: message
 273 integer :: change_water , change_land
 274 change_water = 0
 275 change_land = 0
 276
 277       !  Sanity check on the 50/50 points
 278
 279       DO j = jts , MIN(jde-1,jte)
 280          DO i = its , MIN(ide-1,ite)
 281             dominant_value = landuse_frac(i,iswater,j)
 282             IF ( dominant_value .EQ. lwthresh ) THEN
 283                DO l = 1 , num_veg_cat
 284                   IF ( l .EQ. iswater ) CYCLE
 285                   IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
 286                      PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
 287                      landuse_frac(i,l,j) = lwthresh - .01
 288                      landuse_frac(i,iswater,j) = lwthresh + 0.01
 289                   ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
 290                      PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
 291                      landuse_frac(i,l,j) = lwthresh + .01
 292                      landuse_frac(i,iswater,j) = lwthresh - 0.01
 293                   END IF
 294                END DO
 295             END IF
 296          END DO
 297       END DO
 298
 299       !  Compute the dominant VEGETATION INDEX.
 300
 301       DO j = jts , MIN(jde-1,jte)
 302          DO i = its , MIN(ide-1,ite)
 303             dominant_value = landuse_frac(i,1,j)
 304             dominant_index = 1
 305             DO l = 2 , num_veg_cat
 306                IF        ( l .EQ. iswater ) THEN
 307                   ! wait a bit
 308                ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
 309                   dominant_value = landuse_frac(i,l,j)
 310                   dominant_index = l
 311                END IF
 312             END DO
 313             IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
 314                dominant_value = landuse_frac(i,iswater,j)
 315                dominant_index = iswater
 316             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
 317                       ( landmask(i,j) .LT. 0.5) .AND. &
 318                       ( dominant_value .EQ. lwthresh) ) THEN
 319                dominant_value = landuse_frac(i,iswater,j)
 320                dominant_index = iswater
 321             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
 322                       ( landmask(i,j) .GT. 0.5) .AND. &
 323                       ( dominant_value .EQ. lwthresh) ) THEN
 324                !no op
 325             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
 326                       ( dominant_value .LT. lwthresh ) ) THEN
 327                dominant_value = landuse_frac(i,iswater,j)
 328                dominant_index = iswater
 329             END IF
 330             IF      ( dominant_index .EQ. iswater ) THEN
 331 if(landmask(i,j).gt.lwthresh) then
 332 !print *,'changing to water at point ',i,j
 333 !print '(24(i3,1x))',1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16, 17,18,19,20,21, 22, 23,24
 334 !print '(24(i3,1x))',nint(landuse_frac(i,:,j)*100)
 335 change_water=change_water+1
 336 endif
 337                landmask(i,j) = 0
 338             ELSE IF ( dominant_index .NE. iswater ) THEN
 339 if(landmask(i,j).lt.lwthresh) then
 340 !print *,'changing to land at point ',i,j
 341 !print '(24(i3,1x))',1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16, 17,18,19,20,21, 22, 23,24
 342 !print '(24(i3,1x))',nint(landuse_frac(i,:,j)*100)
 343 change_land=change_land+1
 344 endif
 345                landmask(i,j) = 1
 346             END IF
 347             ivgtyp(i,j) = dominant_index
 348          END DO
 349       END DO
 350
 351       !  Compute the dominant SOIL TEXTURE INDEX, TOP.
 352
 353       iforce = 0
 354       DO i = its , MIN(ide-1,ite)
 355          DO j = jts , MIN(jde-1,jte)
 356             dominant_value = soil_top_cat(i,1,j)
 357             dominant_index = 1
 358             IF ( landmask(i,j) .GT. lwthresh ) THEN
 359                DO l = 2 , num_soil_top_cat
 360                   IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
 361                      dominant_value = soil_top_cat(i,l,j)
 362                      dominant_index = l
 363                   END IF
 364                END DO
 365                IF ( dominant_value .LT. 0.01 ) THEN
 366                   iforce = iforce + 1
 367                   WRITE ( message , FMT = '(A,I4,I4)' ) &
 368                   'based on landuse, changing soil to land at point ',i,j
 369                   CALL wrf_debug(1,message)
 370                   WRITE ( message , FMT = '(16(i3,1x))' ) &
 371                   1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
 372                   CALL wrf_debug(1,message)
 373                   WRITE ( message , FMT = '(16(i3,1x))' ) &
 374                   nint(soil_top_cat(i,:,j)*100)
 375                   CALL wrf_debug(1,message)
 376                   dominant_index = 8
 377                END IF
 378             ELSE
 379                dominant_index = iswater_soil
 380             END IF
 381             isltyp(i,j) = dominant_index
 382          END DO
 383       END DO
 384
 385 if(iforce.ne.0)then
 386 WRITE(message,FMT='(A,I4,A,I6)' ) &
 387 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
 388 (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
 389 CALL wrf_debug(0,message)
 390 endif
 391 print *,'LAND  CHANGE = ',change_land
 392 print *,'WATER CHANGE = ',change_water
 393
 394    END SUBROUTINE process_percent_cat_new
 395
 396    SUBROUTINE process_soil_real ( tsk , tmn , &
 397                                 landmask , sst , &
 398                                 st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
 399                                 zs , dzs , tslb , smois , sh2o , &
 400                                 flag_sst , flag_soilt000, flag_soilm000, &
 401                                 ids , ide , jds , jde , kds , kde , &
 402                                 ims , ime , jms , jme , kms , kme , &
 403                                 its , ite , jts , jte , kts , kte , &
 404                                 sf_surface_physics , num_soil_layers , real_data_init_type , &
 405                                 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 406                                 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
 407
 408       IMPLICIT NONE
 409
 410       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
 411                               ims , ime , jms , jme , kms , kme , &
 412                               its , ite , jts , jte , kts , kte , &
 413                               sf_surface_physics , num_soil_layers , real_data_init_type , &
 414                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 415                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
 416
 417       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
 418
 419       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
 420
 421       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
 422       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
 423       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
 424       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
 425       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
 426       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
 427
 428       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
 429       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
 430       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
 431       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
 432
 433       INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
 434       REAL :: dominant_value
 435
 436       !  Initialize the soil depth, and the soil temperature and moisture.
 437
 438       IF      ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 439          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
 440          CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
 441                                  landmask , sst , flag_sst , &
 442                                  ids , ide , jds , jde , kds , kde , &
 443                                  ims , ime , jms , jme , kms , kme , &
 444                                  its , ite , jts , jte , kts , kte )
 445
 446       ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 447          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
 448          CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
 449                                  st_input , sm_input , sw_input , landmask , sst , &
 450                                  zs , dzs , &
 451                                  st_levels_input , sm_levels_input , sw_levels_input , &
 452                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 453                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
 454                                  flag_sst , flag_soilt000 , flag_soilm000 , &
 455                                  ids , ide , jds , jde , kds , kde , &
 456                                  ims , ime , jms , jme , kms , kme , &
 457                                  its , ite , jts , jte , kts , kte )
 458 !        CALL init_soil_old ( tsk , tmn , &
 459 !                             smois , tslb , zs , dzs , num_soil_layers , &
 460 !                             st000010_input , st010040_input , st040100_input , st100200_input , &
 461 !                             st010200_input , &
 462 !                             sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
 463 !                             sm010200_input , &
 464 !                             landmask_input , sst_input , &
 465 !                             ids , ide , jds , jde , kds , kde , &
 466 !                             ims , ime , jms , jme , kms , kme , &
 467 !                             its , ite , jts , jte , kts , kte )
 468       ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 469          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
 470          CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
 471                                  st_input , sm_input , landmask , sst , &
 472                                  zs , dzs , &
 473                                  st_levels_input , sm_levels_input , &
 474                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , &
 475                                  num_st_levels_alloc , num_sm_levels_alloc , &
 476                                  flag_sst , flag_soilt000 , flag_soilm000 , &
 477                                  ids , ide , jds , jde , kds , kde , &
 478                                  ims , ime , jms , jme , kms , kme , &
 479                                  its , ite , jts , jte , kts , kte )
 480       ELSE IF ( ( sf_surface_physics .EQ. 7 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 481          CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
 482          CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
 483                                  st_input , sm_input , sw_input, landmask , sst , &
 484                                  zs , dzs , &
 485                                  st_levels_input , sm_levels_input , sw_levels_input, &
 486                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 487                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
 488                                  flag_sst , flag_soilt000 , flag_soilm000 , &
 489                                  ids , ide , jds , jde , kds , kde , &
 490                                  ims , ime , jms , jme , kms , kme , &
 491                                  its , ite , jts , jte , kts , kte )
 492       END IF
 493
 494    END SUBROUTINE process_soil_real
 495
 496    SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat,  &
 497                                    ivgtyp,isltyp,tslb,smois, &
 498                                    tsk,tmn,zs,dzs,           &
 499                                    num_soil_layers,          &
 500                                    sf_surface_physics ,      &
 501                                    ids,ide, jds,jde, kds,kde,&
 502                                    ims,ime, jms,jme, kms,kme,&
 503                                    its,ite, jts,jte, kts,kte )
 504
 505       IMPLICIT NONE
 506
 507       INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde,  &
 508                             ims,ime, jms,jme, kms,kme,  &
 509                             its,ite, jts,jte, kts,kte
 510
 511       INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
 512
 513       REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
 514
 515       REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
 516
 517       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
 518       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
 519       INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
 520
 521       !  Local variables.
 522
 523       INTEGER :: itf,jtf
 524
 525       itf=MIN(ite,ide-1)
 526       jtf=MIN(jte,jde-1)
 527
 528       IF      ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 529          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
 530          CALL init_soil_1_ideal(tsk,tmn,tslb,xland,                      &
 531                                 ivgtyp,zs,dzs,num_soil_layers,           &
 532                                 ids,ide, jds,jde, kds,kde,               &
 533                                 ims,ime, jms,jme, kms,kme,               &
 534                                 its,ite, jts,jte, kts,kte                )
 535       ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 536          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
 537          CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,         &
 538                                   ivgtyp,isltyp,tslb,smois,tmn,          &
 539                                   num_soil_layers,                       &
 540                                   ids,ide, jds,jde, kds,kde,             &
 541                                   ims,ime, jms,jme, kms,kme,             &
 542                                   its,ite, jts,jte, kts,kte              )
 543       ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
 544          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
 545
 546       END IF
 547
 548    END SUBROUTINE process_soil_ideal
 549
 550    SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
 551                                  tsk , ter , toposoil , landmask , flag_toposoil , flag_tavgsfc , &
 552                                       st000010 ,      st010040 ,      st040100 ,      st100200 ,      st010200 , &
 553                                  flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
 554                                       st000007 ,      st007028 ,      st028100 ,      st100255 , &
 555                                  flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255 , &
 556                                       soilt000 ,      soilt005 ,      soilt020 ,      soilt040 ,      soilt160 ,      soilt300 , &
 557                                  flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
 558                                  ids , ide , jds , jde , kds , kde , &
 559                                  ims , ime , jms , jme , kms , kme , &
 560                                  its , ite , jts , jte , kts , kte )
 561
 562       IMPLICIT NONE
 563
 564       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
 565                               ims , ime , jms , jme , kms , kme , &
 566                               its , ite , jts , jte , kts , kte
 567
 568       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: ter , toposoil , landmask
 569       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
 570       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
 571       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000007 , st007028 , st028100 , st100255
 572       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
 573
 574       INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
 575       INTEGER , INTENT(IN) :: flag_st000007 , flag_st007028 , flag_st028100 , flag_st100255
 576       INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
 577       INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
 578
 579       INTEGER :: i , j
 580
 581       REAL :: soil_elev_min_val ,  soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
 582
 583       !  Adjust the annual mean temperature as if it is based on from a sea-level elevation
 584       !  if the value used is from the actual annula mean data set.  If the input field to
 585       !  be used for tmn is one of the first-guess input temp fields, need to do an adjustment
 586       !  only on the diff in topo from the model terrain and the first-guess terrain.
 587
 588        SELECT CASE ( sf_surface_physics )
 589
 590          CASE (LSMSCHEME)
 591             DO j = jts , MIN(jde-1,jte)
 592                DO i = its , MIN(ide-1,ite)
 593                   IF (landmask(i,j) .GT. 0.5 ) THEN
 594                      tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
 595                   END IF
 596                END DO
 597             END DO
 598
 599          CASE (RUCLSMSCHEME)
 600             DO j = jts , MIN(jde-1,jte)
 601                DO i = its , MIN(ide-1,ite)
 602                   IF (landmask(i,j) .GT. 0.5 ) THEN
 603                      tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
 604                   END IF
 605                END DO
 606             END DO
 607
 608       END SELECT
 609
 610
 611       !  Do we have a soil field with which to modify soil temperatures?
 612
 613       IF ( flag_toposoil .EQ. 1 ) THEN
 614
 615          DO j = jts , MIN(jde-1,jte)
 616             DO i = its , MIN(ide-1,ite)
 617
 618                !  Is the toposoil field OK, or is it a subversive soil elevation field.  We can tell
 619                !  usually by looking at values.  Anything less than -1000 m (lower than the Dead Sea) is
 620                !  bad.  Anything larger than 10 km (taller than Everest) is toast.  Also, anything where
 621                !  the difference between the soil elevation and the terrain is greater than 3 km means
 622                !  that the soil data is either all zeros or that the data are inconsistent.  Any of these
 623                !  three conditions is grievous enough to induce a WRF fatality.  However, if they are at
 624                !  a water point, then we can safely ignore them.
 625
 626                soil_elev_min_val = toposoil(i,j)
 627                soil_elev_max_val = toposoil(i,j)
 628                soil_elev_min_dif = ter(i,j) - toposoil(i,j)
 629                soil_elev_max_dif = ter(i,j) - toposoil(i,j)
 630
 631                IF      ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
 632                   CYCLE
 633                ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
 634 !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
 635 cycle
 636 !                 CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
 637                ENDIF
 638
 639                IF      ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
 640                   CYCLE
 641                ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
 642 print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
 643 cycle
 644                   CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
 645                ENDIF
 646
 647                IF      ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
 648                            ( landmask(i,j) .LT. 0.5 ) ) THEN
 649                   CYCLE
 650                ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
 651                            ( landmask(i,j) .GT. 0.5 ) ) THEN
 652 print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
 653 cycle
 654                   CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
 655                ENDIF
 656
 657                !  For each of the fields that we would like to modify, check to see if it came in from the SI.
 658                !  If so, then use a -6.5 K/km lapse rate (based on the elevation diffs).  We only adjust when we
 659                !  are not at a water point.
 660
 661                IF (landmask(i,j) .GT. 0.5 ) THEN
 662
 663                   IF ( sf_surface_physics .EQ. SLABSCHEME ) THEN
 664                      IF ( ( flag_tavgsfc  .EQ. 1 ) .OR. &
 665                           ( flag_st010040 .EQ. 1 ) .OR. &
 666                           ( flag_st000010 .EQ. 1 ) .OR. &
 667                           ( flag_soilt020 .EQ. 1 ) .OR. &
 668                           ( flag_st007028 .EQ. 1 ) ) THEN
 669                         tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 670                      ELSE
 671                         tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
 672                      END IF
 673                   END IF
 674
 675                   tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 676
 677                   IF ( flag_st000010 .EQ. 1 ) THEN
 678                      st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 679                   END IF
 680                   IF ( flag_st010040 .EQ. 1 ) THEN
 681                      st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 682                   END IF
 683                   IF ( flag_st040100 .EQ. 1 ) THEN
 684                      st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 685                   END IF
 686                   IF ( flag_st100200 .EQ. 1 ) THEN
 687                      st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 688                   END IF
 689                   IF ( flag_st010200 .EQ. 1 ) THEN
 690                      st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 691                   END IF
 692
 693                   IF ( flag_st000007 .EQ. 1 ) THEN
 694                      st000007(i,j) = st000007(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 695                   END IF
 696                   IF ( flag_st007028 .EQ. 1 ) THEN
 697                      st007028(i,j) = st007028(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 698                   END IF
 699                   IF ( flag_st028100 .EQ. 1 ) THEN
 700                      st028100(i,j) = st028100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 701                   END IF
 702                   IF ( flag_st100255 .EQ. 1 ) THEN
 703                      st100255(i,j) = st100255(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 704                   END IF
 705
 706                   IF ( flag_soilt000 .EQ. 1 ) THEN
 707                      soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 708                   END IF
 709                   IF ( flag_soilt005 .EQ. 1 ) THEN
 710                      soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 711                   END IF
 712                   IF ( flag_soilt020 .EQ. 1 ) THEN
 713                      soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 714                   END IF
 715                   IF ( flag_soilt040 .EQ. 1 ) THEN
 716                      soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 717                   END IF
 718                   IF ( flag_soilt160 .EQ. 1 ) THEN
 719                      soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 720                   END IF
 721                   IF ( flag_soilt300 .EQ. 1 ) THEN
 722                      soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
 723                   END IF
 724
 725                END IF
 726             END DO
 727          END DO
 728
 729       END IF
 730
 731    END SUBROUTINE adjust_soil_temp_new
 732
 733
 734    SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
 735
 736       IMPLICIT NONE
 737
 738       INTEGER, INTENT(IN) :: num_soil_layers
 739
 740       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
 741
 742       INTEGER                   ::      l
 743
 744       !  Define layers (top layer = 0.01 m).  Double the thicknesses at each step (dzs values).
 745       !  The distance from the ground level to the midpoint of the layer is given by zs.
 746
 747       !    -------   Ground Level   ----------      ||      ||   ||  ||
 748       !                                             ||      ||   ||  || zs(1) = 0.005 m
 749       !    --  --  --  --  --  --  --  --  --       ||      ||   ||  \/
 750       !                                             ||      ||   ||
 751       !    -----------------------------------      ||  ||  ||   \/   dzs(1) = 0.01 m
 752       !                                             ||  ||  ||
 753       !                                             ||  ||  || zs(2) = 0.02
 754       !    --  --  --  --  --  --  --  --  --       ||  ||  \/
 755       !                                             ||  ||
 756       !                                             ||  ||
 757       !    -----------------------------------  ||  ||  \/   dzs(2) = 0.02 m
 758       !                                         ||  ||
 759       !                                         ||  ||
 760       !                                         ||  ||
 761       !                                         ||  || zs(3) = 0.05
 762       !    --  --  --  --  --  --  --  --  --   ||  \/
 763       !                                         ||
 764       !                                         ||
 765       !                                         ||
 766       !                                         ||
 767       !    -----------------------------------  \/   dzs(3) = 0.04 m
 768
 769       IF ( num_soil_layers .NE. 5 ) THEN
 770          PRINT '(A)','Usually, the 5-layer diffusion uses 5 layers.  Change this in the namelist.'
 771          CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
 772       END IF
 773
 774       dzs(1)=.01
 775       zs(1)=.5*dzs(1)
 776
 777       DO l=2,num_soil_layers
 778          dzs(l)=2*dzs(l-1)
 779          zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
 780       ENDDO
 781
 782    END SUBROUTINE init_soil_depth_1
 783
 784    SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
 785
 786       IMPLICIT NONE
 787
 788       INTEGER, INTENT(IN) :: num_soil_layers
 789
 790       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
 791
 792       INTEGER                   ::      l
 793
 794       dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
 795
 796       IF ( num_soil_layers .NE. 4 ) THEN
 797          PRINT '(A)','Usually, the LSM uses 4 layers.  Change this in the namelist.'
 798          CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
 799       END IF
 800
 801       zs(1)=.5*dzs(1)
 802
 803       DO l=2,num_soil_layers
 804          zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
 805       ENDDO
 806
 807    END SUBROUTINE init_soil_depth_2
 808
 809    SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
 810
 811       IMPLICIT NONE
 812
 813       INTEGER, INTENT(IN) :: num_soil_layers
 814
 815       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
 816
 817       INTEGER                   ::      l
 818
 819       CHARACTER (LEN=132) :: message
 820
 821 ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
 822 ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
 823 ! Other options with number of levels are possible, but
 824 ! WRF users should change consistently the namelist entry with the
 825 !    ZS array in this subroutine.
 826
 827      IF ( num_soil_layers .EQ. 6) THEN
 828       zs  = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
 829 !      dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
 830      ELSEIF ( num_soil_layers .EQ. 9) THEN
 831       zs  = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
 832 !      dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
 833      ENDIF
 834
 835       IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
 836          write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels.  Change this in the namelist.'
 837          CALL wrf_error_fatal ( message )
 838       END IF
 839
 840    END SUBROUTINE init_soil_depth_3
 841
 842    SUBROUTINE init_soil_depth_7 ( zs , dzs , num_soil_layers )
 843
 844       IMPLICIT NONE
 845
 846       INTEGER, INTENT(IN) :: num_soil_layers
 847
 848       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
 849
 850       INTEGER                   ::      l
 851
 852       dzs = (/ 0.01 , 0.99 /)
 853
 854       IF ( num_soil_layers .NE. 2 ) THEN
 855          PRINT '(A)','Usually, the PX LSM uses 2 layers.  Change this in the namelist.'
 856          CALL wrf_error_fatal ( 'PXLSM_uses_2_layers' )
 857       END IF
 858
 859       zs(1) = 0.5 * dzs(1)
 860       zs(2) = dzs(1) + 0.5 * dzs(2)
 861
 862    END SUBROUTINE init_soil_depth_7
 863
 864    SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
 865                                  num_soil_layers , real_data_init_type , &
 866                                  landmask , sst , flag_sst , &
 867                                  ids , ide , jds , jde , kds , kde , &
 868                                  ims , ime , jms , jme , kms , kme , &
 869                                  its , ite , jts , jte , kts , kte )
 870
 871       IMPLICIT NONE
 872
 873       INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
 874                               ids , ide , jds , jde , kds , kde , &
 875                               ims , ime , jms , jme , kms , kme , &
 876                               its , ite , jts , jte , kts , kte
 877
 878       INTEGER , INTENT(IN) :: flag_sst
 879
 880       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
 881       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
 882
 883       REAL , DIMENSION(num_soil_layers) :: zs , dzs
 884
 885       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
 886
 887       INTEGER :: i , j , l
 888
 889       !  Soil temperature is linearly interpolated between the skin temperature (taken to be at a
 890       !  depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
 891       !  The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
 892       !  annual mean temperature.
 893
 894       DO j = jts , MIN(jde-1,jte)
 895          DO i = its , MIN(ide-1,ite)
 896             IF ( landmask(i,j) .GT. 0.5 ) THEN
 897                DO l = 1 , num_soil_layers
 898                   tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) )   + &
 899                                  tmn(i,j) * ( zs(              l) - zs(1) ) ) / &
 900                                             ( zs(num_soil_layers) - zs(1) )
 901                END DO
 902             ELSE
 903                IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
 904                   DO l = 1 , num_soil_layers
 905                      tslb(i,l,j)= sst(i,j)
 906                   END DO
 907                ELSE
 908                   DO l = 1 , num_soil_layers
 909                      tslb(i,l,j)= tsk(i,j)
 910                   END DO
 911                END IF
 912             END IF
 913          END DO
 914       END DO
 915
 916    END SUBROUTINE init_soil_1_real
 917
 918    SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland,             &
 919                        ivgtyp,ZS,DZS,num_soil_layers,           &
 920                        ids,ide, jds,jde, kds,kde,               &
 921                        ims,ime, jms,jme, kms,kme,               &
 922                        its,ite, jts,jte, kts,kte                )
 923
 924       IMPLICIT NONE
 925
 926       INTEGER, INTENT(IN   )    ::      ids,ide, jds,jde, kds,kde, &
 927                                         ims,ime, jms,jme, kms,kme, &
 928                                         its,ite, jts,jte, kts,kte
 929
 930       INTEGER, INTENT(IN   )    ::      num_soil_layers
 931
 932       REAL, DIMENSION( ims:ime , num_soil_layers , jms:jme ), INTENT(OUT) :: tslb
 933       REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
 934       INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
 935
 936       REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
 937
 938       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
 939
 940       !  Lcal variables.
 941
 942       INTEGER :: l,j,i,itf,jtf
 943
 944       itf=MIN(ite,ide-1)
 945       jtf=MIN(jte,jde-1)
 946
 947       IF (num_soil_layers.NE.1)THEN
 948          DO j=jts,jtf
 949             DO l=1,num_soil_layers
 950                DO i=its,itf
 951                  tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
 952                              ( zs(num_soil_layers)-zs(1) )
 953                ENDDO
 954             ENDDO
 955          ENDDO
 956       ENDIF
 957
 958 !     DO j=jts,jtf
 959 !        DO i=its,itf
 960 !          xland(i,j)  = 2
 961 !          ivgtyp(i,j) = 7
 962 !        ENDDO
 963 !     ENDDO
 964
 965    END SUBROUTINE init_soil_1_ideal
 966
 967    SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
 968                                  st_input , sm_input , sw_input , landmask , sst , &
 969                                  zs , dzs , &
 970                                  st_levels_input , sm_levels_input , sw_levels_input , &
 971                                  num_soil_layers , num_st_levels_input , num_sm_levels_input ,  num_sw_levels_input , &
 972                                  num_st_levels_alloc , num_sm_levels_alloc ,  num_sw_levels_alloc , &
 973                                  flag_sst , flag_soilt000 , flag_soilm000 , &
 974                                  ids , ide , jds , jde , kds , kde , &
 975                                  ims , ime , jms , jme , kms , kme , &
 976                                  its , ite , jts , jte , kts , kte )
 977
 978       IMPLICIT NONE
 979
 980       INTEGER , INTENT(IN) :: num_soil_layers , &
 981                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
 982                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
 983                               ids , ide , jds , jde , kds , kde , &
 984                               ims , ime , jms , jme , kms , kme , &
 985                               its , ite , jts , jte , kts , kte
 986
 987       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
 988
 989       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
 990       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
 991       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
 992
 993       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
 994       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
 995       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
 996       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
 997
 998       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
 999       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1000       REAL , DIMENSION(num_soil_layers) :: zs , dzs
1001
1002       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1003
1004       REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1005
1006       INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1007       REAL :: temp
1008       LOGICAL :: found_levels
1009
1010       CHARACTER (LEN=132) :: message
1011
1012       !  Are there any soil temp and moisture levels - ya know, they are mandatory.
1013
1014       num = num_st_levels_input * num_sm_levels_input
1015
1016       IF ( num .GE. 1 ) THEN
1017
1018          !  Ordered levels that we have data for.
1019
1020 !tgs add option to initialize from RUCLSM
1021          IF ( flag_soilt000 .eq. 1 ) THEN
1022            write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1023            CALL wrf_message ( message )
1024            ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input)  ) )
1025          ELSE
1026            write(message, FMT='(A)') ' Assume Noah LSM input'
1027            CALL wrf_message ( message )
1028          ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1029          END IF
1030
1031
1032          !  Sort the levels for temperature.
1033
1034          outert : DO lout = 1 , num_st_levels_input-1
1035             innert : DO lin = lout+1 , num_st_levels_input
1036                IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1037                   temp = st_levels_input(lout)
1038                   st_levels_input(lout) = st_levels_input(lin)
1039                   st_levels_input(lin) = NINT(temp)
1040                   DO j = jts , MIN(jde-1,jte)
1041                      DO i = its , MIN(ide-1,ite)
1042                         temp = st_input(i,lout+1,j)
1043                         st_input(i,lout+1,j) = st_input(i,lin+1,j)
1044                         st_input(i,lin+1,j) = temp
1045                      END DO
1046                   END DO
1047                END IF
1048             END DO innert
1049          END DO outert
1050 !tgs add IF
1051       IF ( flag_soilt000 .NE. 1 ) THEN
1052          DO j = jts , MIN(jde-1,jte)
1053             DO i = its , MIN(ide-1,ite)
1054                st_input(i,1,j) = tsk(i,j)
1055                st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1056             END DO
1057          END DO
1058       ENDIF
1059
1060          !  Sort the levels for moisture.
1061
1062          outerm: DO lout = 1 , num_sm_levels_input-1
1063             innerm : DO lin = lout+1 , num_sm_levels_input
1064                IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1065                   temp = sm_levels_input(lout)
1066                   sm_levels_input(lout) = sm_levels_input(lin)
1067                   sm_levels_input(lin) = NINT(temp)
1068                   DO j = jts , MIN(jde-1,jte)
1069                      DO i = its , MIN(ide-1,ite)
1070                         temp = sm_input(i,lout+1,j)
1071                         sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1072                         sm_input(i,lin+1,j) = temp
1073                      END DO
1074                   END DO
1075                END IF
1076             END DO innerm
1077          END DO outerm
1078 !tgs add IF
1079       IF ( flag_soilm000 .NE. 1 ) THEN
1080          DO j = jts , MIN(jde-1,jte)
1081             DO i = its , MIN(ide-1,ite)
1082                sm_input(i,1,j) = sm_input(i,2,j)
1083                sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1084             END DO
1085          END DO
1086       ENDIF
1087
1088          !  Sort the levels for liquid moisture.
1089
1090          outerw: DO lout = 1 , num_sw_levels_input-1
1091             innerw : DO lin = lout+1 , num_sw_levels_input
1092                IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1093                   temp = sw_levels_input(lout)
1094                   sw_levels_input(lout) = sw_levels_input(lin)
1095                   sw_levels_input(lin) = NINT(temp)
1096                   DO j = jts , MIN(jde-1,jte)
1097                      DO i = its , MIN(ide-1,ite)
1098                         temp = sw_input(i,lout+1,j)
1099                         sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1100                         sw_input(i,lin+1,j) = temp
1101                      END DO
1102                   END DO
1103                END IF
1104             END DO innerw
1105          END DO outerw
1106          IF ( num_sw_levels_input .GT. 1 ) THEN
1107             DO j = jts , MIN(jde-1,jte)
1108                DO i = its , MIN(ide-1,ite)
1109                   sw_input(i,1,j) = sw_input(i,2,j)
1110                   sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1111                END DO
1112             END DO
1113          END IF
1114
1115          found_levels = .TRUE.
1116
1117       ELSE IF ( ( num .LE. 0 ) .AND. (  start_date .NE. current_date ) ) THEN
1118
1119          found_levels = .FALSE.
1120
1121       ELSE
1122          CALL wrf_error_fatal ( &
1123          'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
1124       END IF
1125
1126       !  Is it OK to continue?
1127
1128       IF ( found_levels ) THEN
1129
1130          !tgs:  Here are the levels that we have from the input for temperature.
1131
1132          IF ( flag_soilt000 .EQ. 1 ) THEN
1133             DO l = 1 , num_st_levels_input
1134                zhave(l) = st_levels_input(l) / 100.
1135             END DO
1136
1137          !  Interpolate between the layers we have (zhave) and those that we want (zs).
1138
1139          z_wantt : DO lwant = 1 , num_soil_layers
1140             z_havet : DO lhave = 1 , num_st_levels_input -1
1141                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1142                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1143                   DO j = jts , MIN(jde-1,jte)
1144                      DO i = its , MIN(ide-1,ite)
1145                         tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1146                                            st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1147                                                                    ( zhave(lhave+1) - zhave(lhave) )
1148                      END DO
1149                   END DO
1150                   EXIT z_havet
1151                END IF
1152             END DO z_havet
1153          END DO z_wantt
1154
1155          ELSE
1156
1157          !  Here are the levels that we have from the input for temperature.  The input levels plus
1158          !  two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1159
1160          zhave(1) = 0.
1161          DO l = 1 , num_st_levels_input
1162             zhave(l+1) = st_levels_input(l) / 100.
1163          END DO
1164          zhave(num_st_levels_input+2) = 300. / 100.
1165
1166          !  Interpolate between the layers we have (zhave) and those that we want (zs).
1167
1168          z_wantt_2: DO lwant = 1 , num_soil_layers
1169             z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
1170                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1171                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1172                   DO j = jts , MIN(jde-1,jte)
1173                      DO i = its , MIN(ide-1,ite)
1174                         tslb(i,lwant,j)= ( st_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1175                                            st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1176                                                                    ( zhave(lhave+1) - zhave(lhave) )
1177                      END DO
1178                   END DO
1179                   EXIT z_havet_2
1180                END IF
1181             END DO z_havet_2
1182          END DO z_wantt_2
1183
1184       END IF
1185
1186 !tgs:
1187       IF ( flag_soilm000 .EQ. 1 ) THEN
1188          DO l = 1 , num_sm_levels_input
1189             zhave(l) = sm_levels_input(l) / 100.
1190          END DO
1191
1192       !  Interpolate between the layers we have (zhave) and those that we want (zs).
1193
1194       z_wantm : DO lwant = 1 , num_soil_layers
1195          z_havem : DO lhave = 1 , num_sm_levels_input -1
1196             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1197                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1198                DO j = jts , MIN(jde-1,jte)
1199                   DO i = its , MIN(ide-1,ite)
1200                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1201                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1202                                                                  ( zhave(lhave+1) - zhave(lhave) )
1203                   END DO
1204                END DO
1205                EXIT z_havem
1206             END IF
1207          END DO z_havem
1208       END DO z_wantm
1209
1210       ELSE
1211          !  Here are the levels that we have from the input for moisture.  The input levels plus
1212          !  two more: a value at 0 cm and one at 300 cm.  The 0 cm value is taken to be identical
1213          !  to the most shallow layer's value.  Similarly, the 300 cm value is taken to be the same
1214          !  as the most deep layer's value.
1215
1216          zhave(1) = 0.
1217          DO l = 1 , num_sm_levels_input
1218             zhave(l+1) = sm_levels_input(l) / 100.
1219          END DO
1220          zhave(num_sm_levels_input+2) = 300. / 100.
1221
1222          !  Interpolate between the layers we have (zhave) and those that we want (zs).
1223
1224          z_wantm_2 : DO lwant = 1 , num_soil_layers
1225             z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
1226                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1227                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1228                   DO j = jts , MIN(jde-1,jte)
1229                      DO i = its , MIN(ide-1,ite)
1230                         smois(i,lwant,j)= ( sm_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1231                                             sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1232                                                                     ( zhave(lhave+1) - zhave(lhave) )
1233                      END DO
1234                   END DO
1235                   EXIT z_havem_2
1236                END IF
1237             END DO z_havem_2
1238          END DO z_wantm_2
1239        ENDIF
1240
1241          !  Any liquid soil moisture to worry about?
1242
1243          IF ( num_sw_levels_input .GT. 1 ) THEN
1244
1245             zhave(1) = 0.
1246             DO l = 1 , num_sw_levels_input
1247                zhave(l+1) = sw_levels_input(l) / 100.
1248             END DO
1249             zhave(num_sw_levels_input+2) = 300. / 100.
1250
1251             !  Interpolate between the layers we have (zhave) and those that we want (zs).
1252
1253             z_wantw : DO lwant = 1 , num_soil_layers
1254                z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1255                   IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1256                        ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1257                      DO j = jts , MIN(jde-1,jte)
1258                         DO i = its , MIN(ide-1,ite)
1259                            sh2o(i,lwant,j)= ( sw_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1260                                                sw_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1261                                                                        ( zhave(lhave+1) - zhave(lhave) )
1262                         END DO
1263                      END DO
1264                      EXIT z_havew
1265                   END IF
1266                END DO z_havew
1267             END DO z_wantw
1268
1269          END IF
1270
1271
1272          !  Over water, put in reasonable values for soil temperature and moisture.  These won't be
1273          !  used, but they will make a more continuous plot.
1274
1275          IF ( flag_sst .EQ. 1 ) THEN
1276             DO j = jts , MIN(jde-1,jte)
1277                DO i = its , MIN(ide-1,ite)
1278                   IF ( landmask(i,j) .LT. 0.5 ) THEN
1279                      DO l = 1 , num_soil_layers
1280                         tslb(i,l,j)= sst(i,j)
1281                         smois(i,l,j)= 1.0
1282                         sh2o (i,l,j)= 1.0
1283                      END DO
1284                   END IF
1285                END DO
1286             END DO
1287          ELSE
1288             DO j = jts , MIN(jde-1,jte)
1289                DO i = its , MIN(ide-1,ite)
1290                   IF ( landmask(i,j) .LT. 0.5 ) THEN
1291                      DO l = 1 , num_soil_layers
1292                         tslb(i,l,j)= tsk(i,j)
1293                         smois(i,l,j)= 1.0
1294                         sh2o (i,l,j)= 1.0
1295                      END DO
1296                   END IF
1297                END DO
1298             END DO
1299          END IF
1300
1301          DEALLOCATE (zhave)
1302
1303       END IF
1304
1305    END SUBROUTINE init_soil_2_real
1306
1307    SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,     &
1308                      ivgtyp,isltyp,tslb,smois,tmn,                  &
1309                      num_soil_layers,                               &
1310                      ids,ide, jds,jde, kds,kde,                     &
1311                      ims,ime, jms,jme, kms,kme,                     &
1312                      its,ite, jts,jte, kts,kte                      )
1313
1314       IMPLICIT NONE
1315
1316       INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde,  &
1317                             ims,ime, jms,jme, kms,kme,  &
1318                             its,ite, jts,jte, kts,kte
1319
1320       INTEGER, INTENT(IN) ::num_soil_layers
1321
1322       REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
1323
1324       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT)  :: xland, snow, canwat, xice, vegfra, tmn
1325
1326       INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
1327
1328       INTEGER :: icm,jcm,itf,jtf
1329       INTEGER ::  i,j,l
1330
1331       itf=min0(ite,ide-1)
1332       jtf=min0(jte,jde-1)
1333
1334       icm = ide/2
1335       jcm = jde/2
1336
1337       DO j=jts,jtf
1338          DO l=1,num_soil_layers
1339             DO i=its,itf
1340                if (xland(i,j) .lt. 1.5) then
1341                smois(i,1,j)=0.30
1342                smois(i,2,j)=0.30
1343                smois(i,3,j)=0.30
1344                smois(i,4,j)=0.30
1345
1346                tslb(i,1,j)=290.
1347                tslb(i,2,j)=290.
1348                tslb(i,3,j)=290.
1349                tslb(i,4,j)=290.
1350                endif
1351             ENDDO
1352          ENDDO
1353       ENDDO
1354
1355    END SUBROUTINE init_soil_2_ideal
1356
1357    SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
1358                                  st_input , sm_input , landmask, sst, &
1359                                  zs , dzs , &
1360                                  st_levels_input , sm_levels_input , &
1361                                  num_soil_layers , num_st_levels_input , num_sm_levels_input ,  &
1362                                  num_st_levels_alloc , num_sm_levels_alloc , &
1363                                  flag_sst , flag_soilt000 , flag_soilm000 , &
1364                                  ids , ide , jds , jde , kds , kde , &
1365                                  ims , ime , jms , jme , kms , kme , &
1366                                  its , ite , jts , jte , kts , kte )
1367
1368       IMPLICIT NONE
1369
1370       INTEGER , INTENT(IN) :: num_soil_layers , &
1371                               num_st_levels_input , num_sm_levels_input , &
1372                               num_st_levels_alloc , num_sm_levels_alloc , &
1373                               ids , ide , jds , jde , kds , kde , &
1374                               ims , ime , jms , jme , kms , kme , &
1375                               its , ite , jts , jte , kts , kte
1376
1377       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
1378
1379       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1380       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1381
1382       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1383       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1384       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1385
1386       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1387       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1388       REAL , DIMENSION(num_soil_layers) :: zs , dzs
1389
1390       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
1391
1392       REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1393
1394       INTEGER :: i , j , l , lout , lin , lwant , lhave, k
1395       REAL :: temp
1396
1397       CHARACTER (LEN=132) :: message
1398
1399       !  Allocate the soil layer array used for interpolating.
1400
1401       IF ( ( num_st_levels_input .LE. 0 ) .OR. &
1402            ( num_sm_levels_input .LE. 0 ) ) THEN
1403          write (message, FMT='(A)')&
1404 'No input soil level data (either temperature or moisture, or both are missing).  Required for RUC LSM.'
1405          CALL wrf_error_fatal ( message )
1406       ELSE
1407          IF ( flag_soilt000 .eq. 1 ) THEN
1408            write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
1409            CALL wrf_message ( message )
1410            ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input)  ) )
1411          ELSE
1412            write(message, FMT='(A)') ' Assume non-RUC LSM input'
1413            CALL wrf_message ( message )
1414            ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers)  ) )
1415          END IF
1416       END IF
1417
1418       !  Sort the levels for temperature.
1419
1420       outert : DO lout = 1 , num_st_levels_input-1
1421          innert : DO lin = lout+1 , num_st_levels_input
1422             IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1423                temp = st_levels_input(lout)
1424                st_levels_input(lout) = st_levels_input(lin)
1425                st_levels_input(lin) = NINT(temp)
1426                DO j = jts , MIN(jde-1,jte)
1427                   DO i = its , MIN(ide-1,ite)
1428                      temp = st_input(i,lout,j)
1429                      st_input(i,lout,j) = st_input(i,lin,j)
1430                      st_input(i,lin,j) = temp
1431                   END DO
1432                END DO
1433             END IF
1434          END DO innert
1435       END DO outert
1436
1437       IF ( flag_soilt000 .NE. 1 ) THEN
1438       DO j = jts , MIN(jde-1,jte)
1439          DO i = its , MIN(ide-1,ite)
1440             st_input(i,1,j) = tsk(i,j)
1441             st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1442          END DO
1443       END DO
1444       END IF
1445
1446       !  Sort the levels for moisture.
1447
1448       outerm: DO lout = 1 , num_sm_levels_input-1
1449          innerm : DO lin = lout+1 , num_sm_levels_input
1450             IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1451                temp = sm_levels_input(lout)
1452                sm_levels_input(lout) = sm_levels_input(lin)
1453                sm_levels_input(lin) = NINT(temp)
1454                DO j = jts , MIN(jde-1,jte)
1455                   DO i = its , MIN(ide-1,ite)
1456                      temp = sm_input(i,lout,j)
1457                      sm_input(i,lout,j) = sm_input(i,lin,j)
1458                      sm_input(i,lin,j) = temp
1459                   END DO
1460                END DO
1461             END IF
1462          END DO innerm
1463       END DO outerm
1464
1465       IF ( flag_soilm000 .NE. 1 ) THEN
1466       DO j = jts , MIN(jde-1,jte)
1467          DO i = its , MIN(ide-1,ite)
1468             sm_input(i,1,j) = (sm_input(i,2,j)-sm_input(i,3,j))/   &
1469                               (st_levels_input(2)-st_levels_input(1))*st_levels_input(1)+  &
1470                               sm_input(i,2,j)
1471 !           sm_input(i,1,j) = sm_input(i,2,j)
1472             sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1473          END DO
1474       END DO
1475       END IF
1476
1477       !  Here are the levels that we have from the input for temperature.
1478
1479       IF ( flag_soilt000 .EQ. 1 ) THEN
1480          DO l = 1 , num_st_levels_input
1481             zhave(l) = st_levels_input(l) / 100.
1482          END DO
1483
1484       !  Interpolate between the layers we have (zhave) and those that we want (zs).
1485
1486
1487       z_wantt : DO lwant = 1 , num_soil_layers
1488          z_havet : DO lhave = 1 , num_st_levels_input -1
1489             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1490                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1491                DO j = jts , MIN(jde-1,jte)
1492                   DO i = its , MIN(ide-1,ite)
1493                      tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1494                                         st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1495                                                                 ( zhave(lhave+1) - zhave(lhave) )
1496                   END DO
1497                END DO
1498                EXIT z_havet
1499             END IF
1500          END DO z_havet
1501       END DO z_wantt
1502
1503       ELSE
1504
1505          zhave(1) = 0.
1506          DO l = 1 , num_st_levels_input
1507             zhave(l+1) = st_levels_input(l) / 100.
1508          END DO
1509          zhave(num_st_levels_input+2) = 300. / 100.
1510
1511       !  Interpolate between the layers we have (zhave) and those that we want (zs).
1512
1513       z_wantt_2 : DO lwant = 1 , num_soil_layers
1514          z_havet_2 : DO lhave = 1 , num_st_levels_input +2
1515             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1516                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1517                DO j = jts , MIN(jde-1,jte)
1518                   DO i = its , MIN(ide-1,ite)
1519                      tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1520                                         st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1521                                                                 ( zhave(lhave+1) - zhave(lhave) )
1522                   END DO
1523                END DO
1524                EXIT z_havet_2
1525             END IF
1526          END DO z_havet_2
1527       END DO z_wantt_2
1528
1529       END IF
1530
1531       !  Here are the levels that we have from the input for moisture.
1532
1533       IF ( flag_soilm000 .EQ. 1 ) THEN
1534          DO l = 1 , num_sm_levels_input
1535             zhave(l) = sm_levels_input(l) / 100.
1536          END DO
1537
1538       !  Interpolate between the layers we have (zhave) and those that we want (zs).
1539
1540       z_wantm : DO lwant = 1 , num_soil_layers
1541          z_havem : DO lhave = 1 , num_sm_levels_input -1
1542             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1543                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1544                DO j = jts , MIN(jde-1,jte)
1545                   DO i = its , MIN(ide-1,ite)
1546                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1547                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1548                                                                  ( zhave(lhave+1) - zhave(lhave) )
1549                   END DO
1550                END DO
1551                EXIT z_havem
1552             END IF
1553          END DO z_havem
1554       END DO z_wantm
1555
1556       ELSE
1557
1558          zhave(1) = 0.
1559          DO l = 1 , num_sm_levels_input
1560             zhave(l+1) = sm_levels_input(l) / 100.
1561          END DO
1562          zhave(num_sm_levels_input+2) = 300. / 100.
1563
1564       z_wantm_2 : DO lwant = 1 , num_soil_layers
1565          z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
1566             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1567                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1568                DO j = jts , MIN(jde-1,jte)
1569                   DO i = its , MIN(ide-1,ite)
1570                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
1571                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1572                                                                  ( zhave(lhave+1) - zhave(lhave) )
1573                   END DO
1574                END DO
1575                EXIT z_havem_2
1576             END IF
1577          END DO z_havem_2
1578       END DO z_wantm_2
1579
1580       END IF
1581       !  Over water, put in reasonable values for soil temperature and moisture.  These won't be
1582       !  used, but they will make a more continuous plot.
1583
1584       IF ( flag_sst .EQ. 1 ) THEN
1585          DO j = jts , MIN(jde-1,jte)
1586             DO i = its , MIN(ide-1,ite)
1587                IF ( landmask(i,j) .LT. 0.5 ) THEN
1588                   DO l = 1 , num_soil_layers
1589                      tslb(i,l,j) = sst(i,j)
1590                      tsk(i,j)    = sst(i,j)
1591                      smois(i,l,j)= 1.0
1592                   END DO
1593                END IF
1594             END DO
1595          END DO
1596       ELSE
1597          DO j = jts , MIN(jde-1,jte)
1598             DO i = its , MIN(ide-1,ite)
1599                IF ( landmask(i,j) .LT. 0.5 ) THEN
1600                   DO l = 1 , num_soil_layers
1601                      tslb(i,l,j)= tsk(i,j)
1602                      smois(i,l,j)= 1.0
1603                   END DO
1604                END IF
1605             END DO
1606          END DO
1607       END IF
1608
1609       DEALLOCATE (zhave)
1610
1611    END SUBROUTINE init_soil_3_real
1612    SUBROUTINE init_soil_7_real ( tsk , tmn , smois , sh2o , tslb , &
1613                                  st_input , sm_input , sw_input , landmask , sst ,     &
1614                                  zs , dzs ,                                            &
1615                                  st_levels_input , sm_levels_input , sw_levels_input , &
1616                                  num_soil_layers , num_st_levels_input ,               &
1617                                  num_sm_levels_input ,  num_sw_levels_input ,          &
1618                                  num_st_levels_alloc , num_sm_levels_alloc ,           &
1619                                  num_sw_levels_alloc ,                                 &
1620                                  flag_sst , flag_soilt000 , flag_soilmt000 ,           &
1621                                  ids , ide , jds , jde , kds , kde ,                   &
1622                                  ims , ime , jms , jme , kms , kme ,                   &
1623                                  its , ite , jts , jte , kts , kte )
1624
1625       ! for soil temperature and moisture initialization for the PX LSM
1626
1627       IMPLICIT NONE
1628
1629       INTEGER , INTENT(IN) :: num_soil_layers , &
1630                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
1631                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
1632                               ids , ide , jds , jde , kds , kde , &
1633                               ims , ime , jms , jme , kms , kme , &
1634                               its , ite , jts , jte , kts , kte
1635
1636       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilmt000
1637
1638       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
1639       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
1640       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
1641
1642       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
1643       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
1644       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
1645       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
1646
1647       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
1648       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
1649       REAL , DIMENSION(num_soil_layers) :: zs , dzs
1650
1651       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
1652
1653       REAL , ALLOCATABLE , DIMENSION(:) :: zhave
1654
1655       INTEGER :: i , j , l , lout , lin , lwant , lhave , num
1656       REAL    :: temp
1657       LOGICAL :: found_levels
1658
1659       !  Are there any soil temp and moisture levels - ya know, they are mandatory.
1660
1661       num = num_st_levels_input * num_sm_levels_input
1662
1663       IF ( num .GE. 1 ) THEN
1664
1665          !  Ordered levels that we have data for.
1666
1667          ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
1668
1669          !  Sort the levels for temperature.
1670          outert : DO lout = 1 , num_st_levels_input-1
1671             innert : DO lin = lout+1 , num_st_levels_input
1672                IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
1673                   temp = st_levels_input(lout)
1674                   st_levels_input(lout) = st_levels_input(lin)
1675                   st_levels_input(lin) = NINT(temp)
1676                   DO j = jts , MIN(jde-1,jte)
1677                      DO i = its , MIN(ide-1,ite)
1678                         temp = st_input(i,lout+1,j)
1679                         st_input(i,lout+1,j) = st_input(i,lin+1,j)
1680                         st_input(i,lin+1,j) = temp
1681                      END DO
1682                   END DO
1683                END IF
1684             END DO innert
1685          END DO outert
1686          DO j = jts , MIN(jde-1,jte)
1687             DO i = its , MIN(ide-1,ite)
1688                st_input(i,1,j) = tsk(i,j)
1689                st_input(i,num_st_levels_input+2,j) = tmn(i,j)
1690             END DO
1691          END DO
1692
1693          !  Sort the levels for moisture.
1694
1695          outerm: DO lout = 1 , num_sm_levels_input-1
1696             innerm : DO lin = lout+1 , num_sm_levels_input
1697               IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
1698                   temp = sm_levels_input(lout)
1699                   sm_levels_input(lout) = sm_levels_input(lin)
1700                   sm_levels_input(lin) = NINT(temp)
1701                   DO j = jts , MIN(jde-1,jte)
1702                      DO i = its , MIN(ide-1,ite)
1703                         temp = sm_input(i,lout+1,j)
1704                         sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
1705                         sm_input(i,lin+1,j) = temp
1706                      END DO
1707                   END DO
1708                END IF
1709             END DO innerm
1710          END DO outerm
1711          DO j = jts , MIN(jde-1,jte)
1712             DO i = its , MIN(ide-1,ite)
1713                sm_input(i,1,j) = sm_input(i,2,j)
1714                sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
1715             END DO
1716          END DO
1717
1718          !  Sort the levels for liquid moisture.
1719
1720          outerw: DO lout = 1 , num_sw_levels_input-1
1721             innerw : DO lin = lout+1 , num_sw_levels_input
1722                IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
1723                   temp = sw_levels_input(lout)
1724                   sw_levels_input(lout) = sw_levels_input(lin)
1725                   sw_levels_input(lin) = NINT(temp)
1726                   DO j = jts , MIN(jde-1,jte)
1727                      DO i = its , MIN(ide-1,ite)
1728                         temp = sw_input(i,lout+1,j)
1729                         sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
1730                         sw_input(i,lin+1,j) = temp
1731                      END DO
1732                   END DO
1733                END IF
1734             END DO innerw
1735          END DO outerw
1736         IF ( num_sw_levels_input .GT. 1 ) THEN
1737             DO j = jts , MIN(jde-1,jte)
1738                DO i = its , MIN(ide-1,ite)
1739                   sw_input(i,1,j) = sw_input(i,2,j)
1740                   sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
1741                END DO
1742             END DO
1743          END IF
1744
1745          found_levels = .TRUE.
1746
1747       ELSE IF ( ( num .LE. 0 ) .AND. (  start_date .NE. current_date ) ) THEN
1748
1749          found_levels = .FALSE.
1750
1751       ELSE
1752          CALL wrf_error_fatal ( &
1753          'No input soil level data (temperature, moisture or liquid, or all are missing). Required for PX LSM.' )
1754       END IF
1755
1756       !  Is it OK to continue?
1757
1758       IF ( found_levels ) THEN
1759
1760          !  Here are the levels that we have from the input for temperature.  The input levels plus
1761          !  two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
1762
1763          zhave(1) = 0.
1764          DO l = 1 , num_st_levels_input
1765             zhave(l+1) = st_levels_input(l) / 100.
1766         END DO
1767          zhave(num_st_levels_input+2) = 300. / 100.
1768
1769          !  Interpolate between the layers we have (zhave) and those that we want (zs).
1770
1771          z_wantt : DO lwant = 1 , num_soil_layers
1772             z_havet : DO lhave = 1 , num_st_levels_input +2 -1
1773                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1774                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1775                   DO j = jts , MIN(jde-1,jte)
1776                      DO i = its , MIN(ide-1,ite)
1777                         tslb(i,lwant,j)= ( st_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1778                                            st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1779                                                                    ( zhave(lhave+1) - zhave(lhave) )
1780                      END DO
1781                   END DO
1782                   EXIT z_havet
1783                END IF
1784             END DO z_havet
1785          END DO z_wantt
1786
1787          !  Here are the levels that we have from the input for moisture.  The input levels plus
1788          !  two more: a value at 0 cm and one at 300 cm.  The 0 cm value is taken to be identical
1789          !  to the most shallow layer's value.  Similarly, the 300 cm value is taken to be the same
1790          !  as the most deep layer's value.
1791
1792          zhave(1) = 0.
1793          DO l = 1 , num_sm_levels_input
1794             zhave(l+1) = sm_levels_input(l) / 100.
1795          END DO
1796         zhave(num_sm_levels_input+2) = 300. / 100.
1797
1798          !  Interpolate between the layers we have (zhave) and those that we want (zs).
1799
1800          z_wantm : DO lwant = 1 , num_soil_layers
1801             z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
1802                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1803                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1804                   DO j = jts , MIN(jde-1,jte)
1805                      DO i = its , MIN(ide-1,ite)
1806                         smois(i,lwant,j)= ( sm_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1807                                             sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1808                                                                     ( zhave(lhave+1) - zhave(lhave) )
1809                      END DO
1810                   END DO
1811                   EXIT z_havem
1812                END IF
1813             END DO z_havem
1814          END DO z_wantm
1815
1816          !  Any liquid soil moisture to worry about?
1817
1818          IF ( num_sw_levels_input .GT. 1 ) THEN
1819
1820             zhave(1) = 0.
1821             DO l = 1 , num_sw_levels_input
1822                zhave(l+1) = sw_levels_input(l) / 100.
1823             END DO
1824             zhave(num_sw_levels_input+2) = 300. / 100.
1825
1826           !  Interpolate between the layers we have (zhave) and those that we want (zs).
1827
1828             z_wantw : DO lwant = 1 , num_soil_layers
1829                z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
1830                   IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
1831                        ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
1832                      DO j = jts , MIN(jde-1,jte)
1833                         DO i = its , MIN(ide-1,ite)
1834                            sh2o(i,lwant,j)= ( sw_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
1835                                                sw_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
1836                                                                        ( zhave(lhave+1) - zhave(lhave) )
1837                         END DO
1838                      END DO
1839                      EXIT z_havew
1840                   END IF
1841                END DO z_havew
1842             END DO z_wantw
1843
1844          END IF
1845
1846
1847         !  Over water, put in reasonable values for soil temperature and moisture.  These won't be
1848         !  used, but they will make a more continuous plot.
1849
1850      DO j = jts , MIN(jde-1,jte)
1851         DO i = its , MIN(ide-1,ite)
1852              tslb(i,1,j)= tsk(i,j)
1853              tslb(i,2,j)= tmn(i,j)
1854         END DO
1855      END DO
1856
1857          IF ( flag_sst .EQ. 1 ) THEN
1858             DO j = jts , MIN(jde-1,jte)
1859                DO i = its , MIN(ide-1,ite)
1860                   IF ( landmask(i,j) .LT. 0.5 ) THEN
1861                      DO l = 1 , num_soil_layers
1862                         tslb(i,l,j)= sst(i,j)
1863                         smois(i,l,j)= 1.0
1864                         sh2o (i,l,j)= 1.0
1865                      END DO
1866                   END IF
1867                END DO
1868             END DO
1869          ELSE
1870             DO j = jts , MIN(jde-1,jte)
1871                DO i = its , MIN(ide-1,ite)
1872                   IF ( landmask(i,j) .LT. 0.5 ) THEN
1873                      DO l = 1 , num_soil_layers
1874                         tslb(i,l,j)= tsk(i,j)
1875                         smois(i,l,j)= 1.0
1876                         sh2o (i,l,j)= 1.0
1877                      END DO
1878                   END IF
1879                END DO
1880             END DO
1881          END IF
1882
1883          DEALLOCATE (zhave)
1884
1885       END IF
1886
1887    END SUBROUTINE init_soil_7_real
1888
1889
1890 END MODULE module_soil_pre
1891
1892 #else
1893
1894 MODULE module_soil_pre
1895
1896    USE module_date_time
1897    USE module_state_description
1898
1899 CONTAINS
1900
1901    SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1902                                       xland , landusef , isltyp , soilcat , soilctop , &
1903                                       soilcbot , tmn , &
1904                                       seaice_threshold , &
1905                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1906                                       iswater , isice , &
1907                                       scheme , &
1908                                       ids , ide , jds , jde , kds , kde , &
1909                                       ims , ime , jms , jme , kms , kme , &
1910                                       its , ite , jts , jte , kts , kte )
1911
1912       IMPLICIT NONE
1913
1914       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1915                               ims , ime , jms , jme , kms , kme , &
1916                               its , ite , jts , jte , kts , kte , &
1917                               iswater , isice
1918       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
1919
1920       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
1921       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
1922       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
1923       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1924       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
1925                                                            vegcat, xland , soilcat , tmn
1926       REAL , INTENT(IN) :: seaice_threshold
1927
1928       INTEGER :: i , j , num_seaice_changes , loop
1929       CHARACTER (LEN=132) :: message
1930
1931       fix_seaice : SELECT CASE ( scheme )
1932
1933          CASE ( SLABSCHEME )
1934             DO j = jts , MIN(jde-1,jte)
1935                DO i = its , MIN(ide-1,ite)
1936                   IF ( xice(i,j) .GT. 200.0 ) THEN
1937                      xice(i,j) = 0.
1938                      num_seaice_changes = num_seaice_changes + 1
1939                   END IF
1940                END DO
1941             END DO
1942             IF ( num_seaice_changes .GT. 0 ) THEN
1943                WRITE ( message , FMT='(A,I6)' ) &
1944                'Total pre number of sea ice locations removed (due to FLAG values) = ', &
1945                num_seaice_changes
1946                CALL wrf_debug ( 0 , message )
1947             END IF
1948             num_seaice_changes = 0
1949             DO j = jts , MIN(jde-1,jte)
1950                DO i = its , MIN(ide-1,ite)
1951                   IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
1952                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
1953                      xice(i,j) = 1.
1954                      num_seaice_changes = num_seaice_changes + 1
1955                      tmn(i,j) = 271.4
1956                      vegcat(i,j)=isice
1957                      lu_index(i,j)=ivgtyp(i,j)
1958                      landmask(i,j)=1.
1959                      xland(i,j)=1.
1960                      DO loop=1,num_veg_cat
1961                         landusef(i,loop,j)=0.
1962                      END DO
1963                      landusef(i,ivgtyp(i,j),j)=1.
1964
1965                      isltyp(i,j) = 16
1966                      soilcat(i,j)=isltyp(i,j)
1967                      DO loop=1,num_soil_top_cat
1968                         soilctop(i,loop,j)=0
1969                      END DO
1970                      DO loop=1,num_soil_bot_cat
1971                         soilcbot(i,loop,j)=0
1972                      END DO
1973                      soilctop(i,isltyp(i,j),j)=1.
1974                      soilcbot(i,isltyp(i,j),j)=1.
1975                   END IF
1976                END DO
1977             END DO
1978             IF ( num_seaice_changes .GT. 0 ) THEN
1979                WRITE ( message , FMT='(A,I6)' ) &
1980                'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
1981                CALL wrf_debug ( 0 , message )
1982             END IF
1983
1984          CASE ( LSMSCHEME )
1985          CASE ( RUCLSMSCHEME )
1986
1987       END SELECT fix_seaice
1988
1989    END SUBROUTINE adjust_for_seaice_pre
1990
1991    SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
1992                                       xland , landusef , isltyp , soilcat , soilctop , &
1993                                       soilcbot , tmn , &
1994                                       tslb , smois , sh2o , &
1995                                       seaice_threshold , &
1996                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1997                                       num_soil_layers , &
1998                                       iswater , isice , &
1999                                       scheme , &
2000                                       ids , ide , jds , jde , kds , kde , &
2001                                       ims , ime , jms , jme , kms , kme , &
2002                                       its , ite , jts , jte , kts , kte )
2003
2004       IMPLICIT NONE
2005
2006       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2007                               ims , ime , jms , jme , kms , kme , &
2008                               its , ite , jts , jte , kts , kte , &
2009                               iswater , isice
2010       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
2011       INTEGER , INTENT(IN) :: num_soil_layers
2012
2013       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
2014       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
2015       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
2016       REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
2017       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
2018       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
2019                                                            vegcat, xland , soilcat , tmn
2020       REAL , INTENT(IN) :: seaice_threshold
2021       REAL :: total_depth , mid_point_depth
2022
2023       INTEGER :: i , j , num_seaice_changes , loop
2024       CHARACTER (LEN=132) :: message
2025
2026       fix_seaice : SELECT CASE ( scheme )
2027
2028          CASE ( SLABSCHEME )
2029
2030          CASE ( LSMSCHEME )
2031             DO j = jts , MIN(jde-1,jte)
2032                DO i = its , MIN(ide-1,ite)
2033                   IF ( xice(i,j) .GT. 200.0 ) THEN
2034                      xice(i,j) = 0.
2035                      num_seaice_changes = num_seaice_changes + 1
2036                   END IF
2037                END DO
2038             END DO
2039             IF ( num_seaice_changes .GT. 0 ) THEN
2040                WRITE ( message , FMT='(A,I6)' ) &
2041                'Total post number of sea ice locations removed (due to FLAG values) = ', &
2042                num_seaice_changes
2043                CALL wrf_debug ( 0 , message )
2044             END IF
2045             num_seaice_changes = 0
2046             DO j = jts , MIN(jde-1,jte)
2047                DO i = its , MIN(ide-1,ite)
2048                   IF ( ( xice(i,j) .GE. 0.5 ) .OR. &
2049                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
2050                      xice(i,j) = 1.
2051                      num_seaice_changes = num_seaice_changes + 1
2052                      tmn(i,j) = 271.16
2053                      vegcat(i,j)=isice
2054                      lu_index(i,j)=ivgtyp(i,j)
2055                      landmask(i,j)=1.
2056                      xland(i,j)=1.
2057                      DO loop=1,num_veg_cat
2058                         landusef(i,loop,j)=0.
2059                      END DO
2060                      landusef(i,ivgtyp(i,j),j)=1.
2061
2062                      isltyp(i,j) = 16
2063                      soilcat(i,j)=isltyp(i,j)
2064                      DO loop=1,num_soil_top_cat
2065                         soilctop(i,loop,j)=0
2066                      END DO
2067                      DO loop=1,num_soil_bot_cat
2068                         soilcbot(i,loop,j)=0
2069                      END DO
2070                      soilctop(i,isltyp(i,j),j)=1.
2071                      soilcbot(i,isltyp(i,j),j)=1.
2072
2073                      total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
2074                      DO loop = 1,num_soil_layers
2075                         mid_point_depth=(total_depth/num_soil_layers)/2. + &
2076                                         (loop-1)*(total_depth/num_soil_layers)
2077                         tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
2078                                             mid_point_depth*tmn(i,j) ) / total_depth
2079                      END DO
2080
2081                      DO loop=1,num_soil_layers
2082                         smois(i,loop,j) = 1.0
2083                         sh2o(i,loop,j)  = 0.0
2084                      END DO
2085                   END IF
2086                END DO
2087             END DO
2088             IF ( num_seaice_changes .GT. 0 ) THEN
2089                WRITE ( message , FMT='(A,I6)' ) &
2090                'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
2091                CALL wrf_debug ( 0 , message )
2092             END IF
2093
2094          CASE ( RUCLSMSCHEME )
2095
2096       END SELECT fix_seaice
2097
2098    END SUBROUTINE adjust_for_seaice_post
2099
2100    SUBROUTINE process_percent_cat_new ( landmask ,  &
2101                                 landuse_frac , soil_top_cat , soil_bot_cat , &
2102                                 isltyp , ivgtyp , &
2103                                 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
2104                                 ids , ide , jds , jde , kds , kde , &
2105                                 ims , ime , jms , jme , kms , kme , &
2106                                 its , ite , jts , jte , kts , kte , &
2107                                 iswater )
2108
2109       IMPLICIT NONE
2110
2111       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2112                               ims , ime , jms , jme , kms , kme , &
2113                               its , ite , jts , jte , kts , kte , &
2114                               iswater
2115       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
2116       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
2117       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
2118       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
2119       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
2120       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
2121
2122       INTEGER :: i , j , l , ll, dominant_index
2123       REAL :: dominant_value
2124
2125 #ifdef WRF_CHEM
2126 !      REAL :: lwthresh = .99
2127       REAL :: lwthresh = .50
2128 #else
2129       REAL :: lwthresh = .50
2130 #endif
2131
2132       INTEGER , PARAMETER :: iswater_soil = 14
2133       INTEGER :: iforce
2134       CHARACTER (LEN=1024) :: message
2135
2136       !  Sanity check on the 50/50 points
2137
2138       DO j = jts , MIN(jde-1,jte)
2139          DO i = its , MIN(ide-1,ite)
2140             dominant_value = landuse_frac(i,iswater,j)
2141             IF ( dominant_value .EQ. lwthresh ) THEN
2142                DO l = 1 , num_veg_cat
2143                   IF ( l .EQ. iswater ) CYCLE
2144                   IF ( landuse_frac(i,l,j) .EQ. lwthresh ) THEN
2145                 write(message, FMT='(I4,I4,A,I4,A,F12.4)') i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
2146                 call wrf_message ( message )
2147                      landuse_frac(i,l,j) = lwthresh - .01
2148                      landuse_frac(i,iswater,j) = 1. - (lwthresh + 0.01)
2149                   END IF
2150                END DO
2151             END IF
2152          END DO
2153       END DO
2154
2155       !  Compute the dominant VEGETATION INDEX.
2156
2157       DO j = jts , MIN(jde-1,jte)
2158          DO i = its , MIN(ide-1,ite)
2159             dominant_value = landuse_frac(i,1,j)
2160             dominant_index = 1
2161             DO l = 2 , num_veg_cat
2162                IF        ( l .EQ. iswater ) THEN
2163                   ! wait a bit
2164                ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
2165                   dominant_value = landuse_frac(i,l,j)
2166                   dominant_index = l
2167                END IF
2168             END DO
2169             IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
2170                dominant_value = landuse_frac(i,iswater,j)
2171                dominant_index = iswater
2172 #if 0
2173 si needs to beef up consistency checks before we can use this part
2174             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
2175                       ( dominant_value .EQ. lwthresh) ) THEN
2176                ! no op
2177 #else
2178 ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.LT.iswater))THEN
2179 write(message,*)'temporary SI landmask fix'
2180 call wrf_message(trim(message))
2181 ! no op
2182 ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.GT.iswater))THEN
2183 write(message,*)'temporary SI landmask fix'
2184 call wrf_message(trim(message))
2185 dominant_value = landuse_frac(i,iswater,j)
2186 dominant_index = iswater
2187 #endif
2188             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
2189                       ( dominant_value .LT. lwthresh ) ) THEN
2190                dominant_value = landuse_frac(i,iswater,j)
2191                dominant_index = iswater
2192             END IF
2193             IF      ( dominant_index .EQ. iswater ) THEN
2194 if(landmask(i,j).gt.lwthresh) then
2195 write(message,*) 'changing to water at point ',i,j
2196 call wrf_message(trim(message))
2197 write(message,*)  nint(landuse_frac(i,:,j)*100)
2198 call wrf_message(trim(message))
2199 endif
2200                landmask(i,j) = 0
2201             ELSE IF ( dominant_index .NE. iswater ) THEN
2202 if(landmask(i,j).lt.lwthresh) then
2203 write(message,*) 'changing to land at point ',i,j
2204 call wrf_message(trim(message))
2205 write(message,*) nint(landuse_frac(i,:,j)*100)
2206 call wrf_message(trim(message))
2207 endif
2208                landmask(i,j) = 1
2209             END IF
2210             ivgtyp(i,j) = dominant_index
2211          END DO
2212       END DO
2213
2214       !  Compute the dominant SOIL TEXTURE INDEX, TOP.
2215
2216       iforce = 0
2217       DO i = its , MIN(ide-1,ite)
2218          DO j = jts , MIN(jde-1,jte)
2219             dominant_value = soil_top_cat(i,1,j)
2220             dominant_index = 1
2221             IF ( landmask(i,j) .GT. lwthresh ) THEN
2222                DO l = 2 , num_soil_top_cat
2223                   IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
2224                      dominant_value = soil_top_cat(i,l,j)
2225                      dominant_index = l
2226                   END IF
2227                END DO
2228                IF ( dominant_value .LT. 0.01 ) THEN
2229                   iforce = iforce + 1
2230                   WRITE ( message , FMT = '(A,I4,I4)' ) &
2231                   'based on landuse, changing soil to land at point ',i,j
2232                   CALL wrf_debug(1,message)
2233                   WRITE ( message , FMT = '(16(i3,1x))' ) &
2234                   1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
2235                   CALL wrf_debug(1,message)
2236                   WRITE ( message , FMT = '(16(i3,1x))' ) &
2237                   nint(soil_top_cat(i,:,j)*100)
2238                   CALL wrf_debug(1,message)
2239                   dominant_index = 8
2240                END IF
2241             ELSE
2242                dominant_index = iswater_soil
2243             END IF
2244             isltyp(i,j) = dominant_index
2245          END DO
2246       END DO
2247
2248 if(iforce.ne.0)then
2249 WRITE(message,FMT='(A,I4,A,I6)' ) &
2250 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
2251 (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
2252 CALL wrf_debug(0,message)
2253 endif
2254
2255    END SUBROUTINE process_percent_cat_new
2256
2257    SUBROUTINE process_soil_real ( tsk , tmn , &
2258                                 landmask , sst , &
2259                                 st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
2260                                 zs , dzs , tslb , smois , sh2o , &
2261                                 flag_sst , flag_soilt000, flag_soilm000, &
2262                                 ids , ide , jds , jde , kds , kde , &
2263                                 ims , ime , jms , jme , kms , kme , &
2264                                 its , ite , jts , jte , kts , kte , &
2265                                 sf_surface_physics , num_soil_layers , real_data_init_type , &
2266                                 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2267                                 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
2268
2269       IMPLICIT NONE
2270
2271       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2272                               ims , ime , jms , jme , kms , kme , &
2273                               its , ite , jts , jte , kts , kte , &
2274                               sf_surface_physics , num_soil_layers , real_data_init_type , &
2275                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2276                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
2277
2278       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2279
2280       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2281
2282       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2283       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2284       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2285       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2286       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2287       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2288
2289       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
2290       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2291       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2292       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2293
2294       INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
2295       REAL :: dominant_value
2296
2297       !  Initialize the soil depth, and the soil temperature and moisture.
2298
2299       IF      ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2300          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2301          CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
2302                                  landmask , sst , flag_sst , &
2303                                  ids , ide , jds , jde , kds , kde , &
2304                                  ims , ime , jms , jme , kms , kme , &
2305                                  its , ite , jts , jte , kts , kte )
2306
2307       ELSE IF ( ( sf_surface_physics .EQ. 2 .or. sf_surface_physics .EQ. 99 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2308          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2309          CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2310                                  st_input , sm_input , sw_input , landmask , sst , &
2311                                  zs , dzs , &
2312                                  st_levels_input , sm_levels_input , sw_levels_input , &
2313                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2314                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2315                                  flag_sst , flag_soilt000 , flag_soilm000 , &
2316                                  ids , ide , jds , jde , kds , kde , &
2317                                  ims , ime , jms , jme , kms , kme , &
2318                                  its , ite , jts , jte , kts , kte )
2319 !        CALL init_soil_old ( tsk , tmn , &
2320 !                             smois , tslb , zs , dzs , num_soil_layers , &
2321 !                             st000010_input , st010040_input , st040100_input , st100200_input , &
2322 !                             st010200_input , &
2323 !                             sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
2324 !                             sm010200_input , &
2325 !                             landmask_input , sst_input , &
2326 !                             ids , ide , jds , jde , kds , kde , &
2327 !                             ims , ime , jms , jme , kms , kme , &
2328 !                             its , ite , jts , jte , kts , kte )
2329       ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2330          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2331          CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
2332                                  st_input , sm_input , landmask , sst , &
2333                                  zs , dzs , &
2334                                  st_levels_input , sm_levels_input , &
2335                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , &
2336                                  num_st_levels_alloc , num_sm_levels_alloc , &
2337                                  flag_sst , flag_soilt000 , flag_soilm000 , &
2338                                  ids , ide , jds , jde , kds , kde , &
2339                                  ims , ime , jms , jme , kms , kme , &
2340                                  its , ite , jts , jte , kts , kte )
2341       END IF
2342
2343    END SUBROUTINE process_soil_real
2344
2345    SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat,  &
2346                                    ivgtyp,isltyp,tslb,smois, &
2347                                    tsk,tmn,zs,dzs,           &
2348                                    num_soil_layers,          &
2349                                    sf_surface_physics ,      &
2350                                    ids,ide, jds,jde, kds,kde,&
2351                                    ims,ime, jms,jme, kms,kme,&
2352                                    its,ite, jts,jte, kts,kte )
2353
2354       IMPLICIT NONE
2355
2356       INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde,  &
2357                             ims,ime, jms,jme, kms,kme,  &
2358                             its,ite, jts,jte, kts,kte
2359
2360       INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
2361
2362       REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
2363
2364       REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
2365
2366       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
2367       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
2368       INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
2369
2370       !  Local variables.
2371
2372       INTEGER :: itf,jtf
2373
2374       itf=MIN(ite,ide-1)
2375       jtf=MIN(jte,jde-1)
2376
2377       IF      ( ( sf_surface_physics .EQ. 1 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2378          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
2379          CALL init_soil_1_ideal(tsk,tmn,tslb,xland,                      &
2380                                 ivgtyp,zs,dzs,num_soil_layers,           &
2381                                 ids,ide, jds,jde, kds,kde,               &
2382                                 ims,ime, jms,jme, kms,kme,               &
2383                                 its,ite, jts,jte, kts,kte                )
2384       ELSE IF ( ( sf_surface_physics .EQ. 2 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2385          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
2386          CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,         &
2387                                   ivgtyp,isltyp,tslb,smois,tmn,          &
2388                                   num_soil_layers,                       &
2389                                   ids,ide, jds,jde, kds,kde,             &
2390                                   ims,ime, jms,jme, kms,kme,             &
2391                                   its,ite, jts,jte, kts,kte              )
2392       ELSE IF ( ( sf_surface_physics .EQ. 3 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
2393          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
2394
2395       END IF
2396
2397    END SUBROUTINE process_soil_ideal
2398
2399    SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
2400                                  tsk , ter , toposoil , landmask , flag_toposoil , &
2401                                       st000010 ,      st010040 ,      st040100 ,      st100200 ,      st010200 , &
2402                                  flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
2403                                       soilt000 ,      soilt005 ,      soilt020 ,      soilt040 ,      soilt160 ,      soilt300 , &
2404                                  flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
2405                                  ids , ide , jds , jde , kds , kde , &
2406                                  ims , ime , jms , jme , kms , kme , &
2407                                  its , ite , jts , jte , kts , kte )
2408
2409       IMPLICIT NONE
2410
2411       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
2412                               ims , ime , jms , jme , kms , kme , &
2413                               its , ite , jts , jte , kts , kte
2414
2415       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: ter , toposoil , landmask
2416       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
2417       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
2418       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
2419
2420       INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
2421       INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
2422       INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
2423
2424       INTEGER :: i , j
2425
2426       REAL :: soil_elev_min_val ,  soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
2427
2428       !  Do we have a soil field with which to modify soil temperatures?
2429
2430       IF ( flag_toposoil .EQ. 1 ) THEN
2431
2432          DO j = jts , MIN(jde-1,jte)
2433             DO i = its , MIN(ide-1,ite)
2434
2435                !  Is the toposoil field OK, or is it a subversive soil elevation field.  We can tell
2436                !  usually by looking at values.  Anything less than -1000 m (lower than the Dead Sea) is
2437                !  bad.  Anything larger than 10 km (taller than Everest) is toast.  Also, anything where
2438                !  the difference between the soil elevation and the terrain is greater than 3 km means
2439                !  that the soil data is either all zeros or that the data are inconsistent.  Any of these
2440                !  three conditions is grievous enough to induce a WRF fatality.  However, if they are at
2441                !  a water point, then we can safely ignore them.
2442
2443                soil_elev_min_val = toposoil(i,j)
2444                soil_elev_max_val = toposoil(i,j)
2445                soil_elev_min_dif = ter(i,j) - toposoil(i,j)
2446                soil_elev_max_dif = ter(i,j) - toposoil(i,j)
2447
2448                IF      ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2449                   CYCLE
2450                ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2451 !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
2452 cycle
2453 !                 CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
2454                ENDIF
2455
2456                IF      ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
2457                   CYCLE
2458                ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
2459 print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
2460 cycle
2461                   CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
2462                ENDIF
2463
2464                IF      ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2465                            ( landmask(i,j) .LT. 0.5 ) ) THEN
2466                   CYCLE
2467                ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
2468                            ( landmask(i,j) .GT. 0.5 ) ) THEN
2469 print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
2470 cycle
2471                   CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
2472                ENDIF
2473
2474                !  For each of the fields that we would like to modify, check to see if it came in from the SI.
2475                !  If so, then use a -6.5 K/km lapse rate (based on the elevation diffs).  We only adjust when we
2476                !  are not at a water point.
2477
2478                IF (landmask(i,j) .GT. 0.5 ) THEN
2479
2480                   IF ( sf_surface_physics .EQ. 1 .OR. sf_surface_physics .EQ. 7) THEN
2481                      tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2482                   END IF
2483
2484                   tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2485
2486                   IF ( flag_st000010 .EQ. 1 ) THEN
2487                      st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2488                   END IF
2489                   IF ( flag_st010040 .EQ. 1 ) THEN
2490                      st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2491                   END IF
2492                   IF ( flag_st040100 .EQ. 1 ) THEN
2493                      st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2494                   END IF
2495                   IF ( flag_st100200 .EQ. 1 ) THEN
2496                      st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2497                   END IF
2498                   IF ( flag_st010200 .EQ. 1 ) THEN
2499                      st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2500                   END IF
2501
2502                   IF ( flag_soilt000 .EQ. 1 ) THEN
2503                      soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2504                   END IF
2505                   IF ( flag_soilt005 .EQ. 1 ) THEN
2506                      soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2507                   END IF
2508                   IF ( flag_soilt020 .EQ. 1 ) THEN
2509                      soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2510                   END IF
2511                   IF ( flag_soilt040 .EQ. 1 ) THEN
2512                      soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2513                   END IF
2514                   IF ( flag_soilt160 .EQ. 1 ) THEN
2515                      soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2516                   END IF
2517                   IF ( flag_soilt300 .EQ. 1 ) THEN
2518                      soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
2519                   END IF
2520
2521                END IF
2522             END DO
2523          END DO
2524
2525       END IF
2526
2527    END SUBROUTINE adjust_soil_temp_new
2528
2529
2530    SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
2531
2532       IMPLICIT NONE
2533
2534       INTEGER, INTENT(IN) :: num_soil_layers
2535
2536       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
2537
2538       INTEGER                   ::      l
2539
2540       CHARACTER (LEN=132) :: message
2541
2542       !  Define layers (top layer = 0.01 m).  Double the thicknesses at each step (dzs values).
2543       !  The distance from the ground level to the midpoint of the layer is given by zs.
2544
2545       !    -------   Ground Level   ----------      ||      ||   ||  ||
2546       !                                             ||      ||   ||  || zs(1) = 0.005 m
2547       !    --  --  --  --  --  --  --  --  --       ||      ||   ||  \/
2548       !                                             ||      ||   ||
2549       !    -----------------------------------      ||  ||  ||   \/   dzs(1) = 0.01 m
2550       !                                             ||  ||  ||
2551       !                                             ||  ||  || zs(2) = 0.02
2552       !    --  --  --  --  --  --  --  --  --       ||  ||  \/
2553       !                                             ||  ||
2554       !                                             ||  ||
2555       !    -----------------------------------  ||  ||  \/   dzs(2) = 0.02 m
2556       !                                         ||  ||
2557       !                                         ||  ||
2558       !                                         ||  ||
2559       !                                         ||  || zs(3) = 0.05
2560       !    --  --  --  --  --  --  --  --  --   ||  \/
2561       !                                         ||
2562       !                                         ||
2563       !                                         ||
2564       !                                         ||
2565       !    -----------------------------------  \/   dzs(3) = 0.04 m
2566
2567       IF ( num_soil_layers .NE. 5 ) THEN
2568          write(message,FMT= '(A)') 'Usually, the 5-layer diffusion uses 5 layers.  Change this in the namelist.'
2569          CALL wrf_error_fatal ( message )
2570       END IF
2571
2572       dzs(1)=.01
2573       zs(1)=.5*dzs(1)
2574
2575       DO l=2,num_soil_layers
2576          dzs(l)=2*dzs(l-1)
2577          zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
2578       ENDDO
2579
2580    END SUBROUTINE init_soil_depth_1
2581
2582    SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
2583
2584       IMPLICIT NONE
2585
2586       INTEGER, INTENT(IN) :: num_soil_layers
2587
2588       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
2589
2590       INTEGER                   ::      l
2591
2592       CHARACTER (LEN=132) :: message
2593
2594       dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
2595
2596       IF ( num_soil_layers .NE. 4 ) THEN
2597          write(message,FMT='(A)') 'Usually, the LSM uses 4 layers.  Change this in the namelist.'
2598          CALL wrf_error_fatal ( message )
2599       END IF
2600
2601       zs(1)=.5*dzs(1)
2602
2603       DO l=2,num_soil_layers
2604          zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
2605       ENDDO
2606
2607    END SUBROUTINE init_soil_depth_2
2608
2609    SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
2610
2611       IMPLICIT NONE
2612
2613       INTEGER, INTENT(IN) :: num_soil_layers
2614
2615       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
2616
2617       INTEGER                   ::      l
2618
2619       CHARACTER (LEN=132) :: message
2620
2621 ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
2622 ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
2623 ! Other options with number of levels are possible, but
2624 ! WRF users should change consistently the namelist entry with the
2625 !    ZS array in this subroutine.
2626
2627      IF ( num_soil_layers .EQ. 6) THEN
2628       zs  = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
2629 !      dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2630      ELSEIF ( num_soil_layers .EQ. 9) THEN
2631       zs  = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
2632 !      dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
2633      ENDIF
2634
2635       IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
2636          WRITE(message,FMT= '(A)')'Usually, the RUC LSM uses 6, 9 or more levels.  Change this in the namelist.'
2637          CALL wrf_error_fatal ( message )
2638       END IF
2639
2640    END SUBROUTINE init_soil_depth_3
2641
2642    SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
2643                                  num_soil_layers , real_data_init_type , &
2644                                  landmask , sst , flag_sst , &
2645                                  ids , ide , jds , jde , kds , kde , &
2646                                  ims , ime , jms , jme , kms , kme , &
2647                                  its , ite , jts , jte , kts , kte )
2648
2649       IMPLICIT NONE
2650
2651       INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
2652                               ids , ide , jds , jde , kds , kde , &
2653                               ims , ime , jms , jme , kms , kme , &
2654                               its , ite , jts , jte , kts , kte
2655
2656       INTEGER , INTENT(IN) :: flag_sst
2657
2658       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2659       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
2660
2661       REAL , DIMENSION(num_soil_layers) :: zs , dzs
2662
2663       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
2664
2665       INTEGER :: i , j , l
2666
2667       !  Soil temperature is linearly interpolated between the skin temperature (taken to be at a
2668       !  depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
2669       !  The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
2670       !  annual mean temperature.
2671
2672       DO j = jts , MIN(jde-1,jte)
2673          DO i = its , MIN(ide-1,ite)
2674             IF ( landmask(i,j) .GT. 0.5 ) THEN
2675                DO l = 1 , num_soil_layers
2676                   tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) )   + &
2677                                  tmn(i,j) * ( zs(              l) - zs(1) ) ) / &
2678                                             ( zs(num_soil_layers) - zs(1) )
2679                END DO
2680             ELSE
2681                IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
2682                   DO l = 1 , num_soil_layers
2683                      tslb(i,l,j)= sst(i,j)
2684                   END DO
2685                ELSE
2686                   DO l = 1 , num_soil_layers
2687                      tslb(i,l,j)= tsk(i,j)
2688                   END DO
2689                END IF
2690             END IF
2691          END DO
2692       END DO
2693
2694    END SUBROUTINE init_soil_1_real
2695
2696    SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland,             &
2697                        ivgtyp,ZS,DZS,num_soil_layers,           &
2698                        ids,ide, jds,jde, kds,kde,               &
2699                        ims,ime, jms,jme, kms,kme,               &
2700                        its,ite, jts,jte, kts,kte                )
2701
2702       IMPLICIT NONE
2703
2704       INTEGER, INTENT(IN   )    ::      ids,ide, jds,jde, kds,kde, &
2705                                         ims,ime, jms,jme, kms,kme, &
2706                                         its,ite, jts,jte, kts,kte
2707
2708       INTEGER, INTENT(IN   )    ::      num_soil_layers
2709
2710       REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
2711       REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
2712       INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
2713
2714       REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
2715
2716       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
2717
2718       !  Lcal variables.
2719
2720       INTEGER :: l,j,i,itf,jtf
2721
2722       itf=MIN(ite,ide-1)
2723       jtf=MIN(jte,jde-1)
2724
2725       IF (num_soil_layers.NE.1)THEN
2726          DO j=jts,jtf
2727             DO l=1,num_soil_layers
2728                DO i=its,itf
2729                  tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
2730                              ( zs(num_soil_layers)-zs(1) )
2731                ENDDO
2732             ENDDO
2733          ENDDO
2734       ENDIF
2735       DO j=jts,jtf
2736          DO i=its,itf
2737            xland(i,j)  = 2
2738            ivgtyp(i,j) = 7
2739          ENDDO
2740       ENDDO
2741
2742    END SUBROUTINE init_soil_1_ideal
2743
2744    SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
2745                                  st_input , sm_input , sw_input , landmask , sst , &
2746                                  zs , dzs , &
2747                                  st_levels_input , sm_levels_input , sw_levels_input , &
2748                                  num_soil_layers , num_st_levels_input , num_sm_levels_input ,  num_sw_levels_input , &
2749                                  num_st_levels_alloc , num_sm_levels_alloc ,  num_sw_levels_alloc , &
2750                                  flag_sst , flag_soilt000 , flag_soilm000 , &
2751                                  ids , ide , jds , jde , kds , kde , &
2752                                  ims , ime , jms , jme , kms , kme , &
2753                                  its , ite , jts , jte , kts , kte )
2754
2755       IMPLICIT NONE
2756
2757       INTEGER , INTENT(IN) :: num_soil_layers , &
2758                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
2759                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
2760                               ids , ide , jds , jde , kds , kde , &
2761                               ims , ime , jms , jme , kms , kme , &
2762                               its , ite , jts , jte , kts , kte
2763
2764       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
2765
2766       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
2767       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
2768       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
2769
2770       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
2771       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
2772       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
2773       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
2774
2775       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
2776       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
2777       REAL , DIMENSION(num_soil_layers) :: zs , dzs
2778
2779       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
2780
2781       REAL , ALLOCATABLE , DIMENSION(:) :: zhave
2782
2783       INTEGER :: i , j , l , lout , lin , lwant , lhave , num
2784       REAL :: temp
2785       LOGICAL :: found_levels
2786
2787       CHARACTER (LEN=132) :: message
2788
2789       !  Are there any soil temp and moisture levels - ya know, they are mandatory.
2790
2791       num = num_st_levels_input * num_sm_levels_input
2792
2793       IF ( num .GE. 1 ) THEN
2794
2795          !  Ordered levels that we have data for.
2796          IF ( flag_soilt000 .eq. 1 ) THEN
2797            write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
2798            CALL wrf_message ( message )
2799            ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input)  ) )
2800          ELSE
2801            write(message, FMT='(A)') ' Assume Noah LSM input'
2802            CALL wrf_message ( message )
2803          ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
2804          END IF
2805
2806
2807          !  Sort the levels for temperature.
2808 !print*,'num_st_levels_input, st_levels_input', num_st_levels_input, st_levels_input
2809 !print*,'num_sm_levels_input,num_sw_levels_input',num_sm_levels_input,num_sw_levels_input
2810
2811          outert : DO lout = 1 , num_st_levels_input-1
2812             innert : DO lin = lout+1 , num_st_levels_input
2813                IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
2814                   temp = st_levels_input(lout)
2815                   st_levels_input(lout) = st_levels_input(lin)
2816                   st_levels_input(lin) = NINT(temp)
2817                   DO j = jts , MIN(jde-1,jte)
2818                      DO i = its , MIN(ide-1,ite)
2819                         temp = st_input(i,lout+1,j)
2820                         st_input(i,lout+1,j) = st_input(i,lin+1,j)
2821                         st_input(i,lin+1,j) = temp
2822                      END DO
2823                   END DO
2824                END IF
2825             END DO innert
2826          END DO outert
2827 !tgs add IF
2828       IF ( flag_soilt000 .NE. 1 ) THEN
2829          DO j = jts , MIN(jde-1,jte)
2830             DO i = its , MIN(ide-1,ite)
2831                st_input(i,1,j) = tsk(i,j)
2832                st_input(i,num_st_levels_input+2,j) = tmn(i,j)
2833             END DO
2834          END DO
2835       ENDIF
2836
2837
2838 #if ( NMM_CORE == 1 )
2839 !new
2840 !       write(0,*) 'st_input(1) in init_soil_2_real'
2841 !       DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2842 !       write(0,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2843 !       ENDDO
2844
2845 !       write(0,*) 'st_input(2) in init_soil_2_real'
2846 !       DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
2847 !       write(0,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
2848 !       ENDDO
2849
2850   616   format(20(f4.0,1x))
2851 #endif
2852
2853          !  Sort the levels for moisture.
2854
2855          outerm: DO lout = 1 , num_sm_levels_input-1
2856             innerm : DO lin = lout+1 , num_sm_levels_input
2857                IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
2858                   temp = sm_levels_input(lout)
2859                   sm_levels_input(lout) = sm_levels_input(lin)
2860                   sm_levels_input(lin) = NINT(temp)
2861                   DO j = jts , MIN(jde-1,jte)
2862                      DO i = its , MIN(ide-1,ite)
2863                         temp = sm_input(i,lout+1,j)
2864                         sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
2865                         sm_input(i,lin+1,j) = temp
2866                      END DO
2867                   END DO
2868                END IF
2869             END DO innerm
2870          END DO outerm
2871 !tgs add IF
2872       IF ( flag_soilm000 .NE. 1 ) THEN
2873          DO j = jts , MIN(jde-1,jte)
2874             DO i = its , MIN(ide-1,ite)
2875                sm_input(i,1,j) = sm_input(i,2,j)
2876                sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
2877             END DO
2878          END DO
2879       ENDIF
2880
2881          !  Sort the levels for liquid moisture.
2882
2883          outerw: DO lout = 1 , num_sw_levels_input-1
2884             innerw : DO lin = lout+1 , num_sw_levels_input
2885                IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
2886                   temp = sw_levels_input(lout)
2887                   sw_levels_input(lout) = sw_levels_input(lin)
2888                   sw_levels_input(lin) = NINT(temp)
2889                   DO j = jts , MIN(jde-1,jte)
2890                      DO i = its , MIN(ide-1,ite)
2891                         temp = sw_input(i,lout+1,j)
2892                         sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
2893                         sw_input(i,lin+1,j) = temp
2894                      END DO
2895                   END DO
2896                END IF
2897             END DO innerw
2898          END DO outerw
2899          IF ( num_sw_levels_input .GT. 1 ) THEN
2900             DO j = jts , MIN(jde-1,jte)
2901                DO i = its , MIN(ide-1,ite)
2902                   sw_input(i,1,j) = sw_input(i,2,j)
2903                   sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
2904                END DO
2905             END DO
2906          END IF
2907
2908          found_levels = .TRUE.
2909
2910       ELSE IF ( ( num .LE. 0 ) .AND. (  start_date .NE. current_date ) ) THEN
2911
2912          found_levels = .FALSE.
2913
2914       ELSE
2915          CALL wrf_error_fatal ( &
2916          'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
2917       END IF
2918
2919       !  Is it OK to continue?
2920
2921       IF ( found_levels ) THEN
2922
2923 !tgs add IF
2924       IF ( flag_soilt000 .NE.1) THEN
2925 !tgs initialize from Noah data
2926          !  Here are the levels that we have from the input for temperature.  The input levels plus
2927          !  two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
2928
2929          zhave(1) = 0.
2930          DO l = 1 , num_st_levels_input
2931             zhave(l+1) = st_levels_input(l) / 100.
2932          END DO
2933          zhave(num_st_levels_input+2) = 300. / 100.
2934
2935          !  Interpolate between the layers we have (zhave) and those that we want (zs).
2936
2937          z_wantt : DO lwant = 1 , num_soil_layers
2938             z_havet : DO lhave = 1 , num_st_levels_input +2 -1
2939                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
2940                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2941                   DO j = jts , MIN(jde-1,jte)
2942                      DO i = its , MIN(ide-1,ite)
2943                         tslb(i,lwant,j)= ( st_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
2944                                            st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
2945                                                                    ( zhave(lhave+1) - zhave(lhave) )
2946                      END DO
2947                   END DO
2948                   EXIT z_havet
2949                END IF
2950             END DO z_havet
2951          END DO z_wantt
2952
2953      ELSE
2954
2955 !tgs initialize from RUCLSM data
2956          DO l = 1 , num_st_levels_input
2957             zhave(l) = st_levels_input(l) / 100.
2958          END DO
2959
2960       !  Interpolate between the layers we have (zhave) and those that we want (zs).
2961
2962
2963       z_wantt_2 : DO lwant = 1 , num_soil_layers
2964          z_havet_2 : DO lhave = 1 , num_st_levels_input -1
2965             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
2966                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
2967                DO j = jts , MIN(jde-1,jte)
2968                   DO i = its , MIN(ide-1,ite)
2969                      tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
2970                                         st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
2971                                                                 ( zhave(lhave+1) - zhave(lhave) )
2972                   END DO
2973                END DO
2974                EXIT z_havet_2
2975             END IF
2976          END DO z_havet_2
2977       END DO z_wantt_2
2978
2979       ENDIF
2980
2981
2982       IF ( flag_soilm000 .NE. 1 ) THEN
2983 !tgs initialize from Noah
2984          !  Here are the levels that we have from the input for moisture.  The input levels plus
2985          !  two more: a value at 0 cm and one at 300 cm.  The 0 cm value is taken to be identical
2986          !  to the most shallow layer's value.  Similarly, the 300 cm value is taken to be the same
2987          !  as the most deep layer's value.
2988
2989          zhave(1) = 0.
2990          DO l = 1 , num_sm_levels_input
2991             zhave(l+1) = sm_levels_input(l) / 100.
2992          END DO
2993          zhave(num_sm_levels_input+2) = 300. / 100.
2994
2995          !  Interpolate between the layers we have (zhave) and those that we want (zs).
2996
2997          z_wantm : DO lwant = 1 , num_soil_layers
2998             z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
2999                IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3000                     ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3001                   DO j = jts , MIN(jde-1,jte)
3002                      DO i = its , MIN(ide-1,ite)
3003                         smois(i,lwant,j)= ( sm_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
3004                                             sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3005                                                                     ( zhave(lhave+1) - zhave(lhave) )
3006                      END DO
3007                   END DO
3008                   EXIT z_havem
3009                END IF
3010             END DO z_havem
3011          END DO z_wantm
3012
3013      ELSE
3014
3015 !tgs  initialize from RUCLSM data
3016          DO l = 1 , num_sm_levels_input
3017             zhave(l) = sm_levels_input(l) / 100.
3018          END DO
3019
3020       !  Interpolate between the layers we have (zhave) and those that we want (zs).
3021
3022       z_wantm_2 : DO lwant = 1 , num_soil_layers
3023          z_havem_2 : DO lhave = 1 , num_sm_levels_input -1
3024             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3025                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3026                DO j = jts , MIN(jde-1,jte)
3027                   DO i = its , MIN(ide-1,ite)
3028                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
3029                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3030                                                                  ( zhave(lhave+1) - zhave(lhave) )
3031                   END DO
3032                END DO
3033                EXIT z_havem_2
3034              END IF
3035          END DO z_havem_2
3036       END DO z_wantm_2
3037
3038      ENDIF
3039          !  Any liquid soil moisture to worry about?
3040
3041          IF ( num_sw_levels_input .GT. 1 ) THEN
3042
3043             zhave(1) = 0.
3044             DO l = 1 , num_sw_levels_input
3045                zhave(l+1) = sw_levels_input(l) / 100.
3046             END DO
3047             zhave(num_sw_levels_input+2) = 300. / 100.
3048
3049             !  Interpolate between the layers we have (zhave) and those that we want (zs).
3050
3051             z_wantw : DO lwant = 1 , num_soil_layers
3052                z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
3053                   IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3054                        ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3055                      DO j = jts , MIN(jde-1,jte)
3056                         DO i = its , MIN(ide-1,ite)
3057                            sh2o(i,lwant,j)= ( sw_input(i,lhave  ,j) * ( zhave(lhave+1) - zs   (lwant) ) + &
3058                                                sw_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3059                                                                        ( zhave(lhave+1) - zhave(lhave) )
3060                         END DO
3061                      END DO
3062                      EXIT z_havew
3063                   END IF
3064                END DO z_havew
3065             END DO z_wantw
3066
3067          END IF
3068
3069
3070          !  Over water, put in reasonable values for soil temperature and moisture.  These won't be
3071          !  used, but they will make a more continuous plot.
3072
3073          IF ( flag_sst .EQ. 1 ) THEN
3074             DO j = jts , MIN(jde-1,jte)
3075                DO i = its , MIN(ide-1,ite)
3076                   IF ( landmask(i,j) .LT. 0.5 ) THEN
3077                      DO l = 1 , num_soil_layers
3078                         tslb(i,l,j)= sst(i,j)
3079 !tgs add line for tsk
3080                         tsk(i,j)    = sst(i,j)
3081                         smois(i,l,j)= 1.0
3082                         sh2o (i,l,j)= 1.0
3083                      END DO
3084                   END IF
3085                END DO
3086             END DO
3087          ELSE
3088             DO j = jts , MIN(jde-1,jte)
3089                DO i = its , MIN(ide-1,ite)
3090                   IF ( landmask(i,j) .LT. 0.5 ) THEN
3091                      DO l = 1 , num_soil_layers
3092                         tslb(i,l,j)= tsk(i,j)
3093                         smois(i,l,j)= 1.0
3094                         sh2o (i,l,j)= 1.0
3095                      END DO
3096                   END IF
3097                END DO
3098             END DO
3099          END IF
3100
3101          DEALLOCATE (zhave)
3102
3103       END IF
3104
3105    END SUBROUTINE init_soil_2_real
3106
3107    SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,     &
3108                      ivgtyp,isltyp,tslb,smois,tmn,                  &
3109                      num_soil_layers,                               &
3110                      ids,ide, jds,jde, kds,kde,                     &
3111                      ims,ime, jms,jme, kms,kme,                     &
3112                      its,ite, jts,jte, kts,kte                      )
3113
3114       IMPLICIT NONE
3115
3116       INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde,  &
3117                             ims,ime, jms,jme, kms,kme,  &
3118                             its,ite, jts,jte, kts,kte
3119
3120       INTEGER, INTENT(IN) ::num_soil_layers
3121
3122       REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
3123
3124       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT)  :: xland, snow, canwat, xice, vegfra, tmn
3125
3126       INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
3127
3128       INTEGER :: icm,jcm,itf,jtf
3129       INTEGER ::  i,j,l
3130
3131       itf=min0(ite,ide-1)
3132       jtf=min0(jte,jde-1)
3133
3134       icm = ide/2
3135       jcm = jde/2
3136
3137       DO j=jts,jtf
3138          DO l=1,num_soil_layers
3139             DO i=its,itf
3140
3141                smois(i,1,j)=0.10
3142                smois(i,2,j)=0.10
3143                smois(i,3,j)=0.10
3144                smois(i,4,j)=0.10
3145
3146                tslb(i,1,j)=295.
3147                tslb(i,2,j)=297.
3148                tslb(i,3,j)=293.
3149                tslb(i,4,j)=293.
3150
3151             ENDDO
3152          ENDDO
3153       ENDDO
3154
3155       DO j=jts,jtf
3156          DO i=its,itf
3157             xland(i,j)  =   2
3158             tmn(i,j)    = 294.
3159             xice(i,j)   =   0.
3160             vegfra(i,j) =   0.
3161             snow(i,j)   =   0.
3162             canwat(i,j) =   0.
3163             ivgtyp(i,j) =   7
3164             isltyp(i,j) =   8
3165          ENDDO
3166       ENDDO
3167
3168    END SUBROUTINE init_soil_2_ideal
3169
3170    SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
3171                                  st_input , sm_input , landmask, sst, &
3172                                  zs , dzs , &
3173                                  st_levels_input , sm_levels_input , &
3174                                  num_soil_layers , num_st_levels_input , num_sm_levels_input ,  &
3175                                  num_st_levels_alloc , num_sm_levels_alloc , &
3176                                  flag_sst , flag_soilt000 , flag_soilm000 , &
3177                                  ids , ide , jds , jde , kds , kde , &
3178                                  ims , ime , jms , jme , kms , kme , &
3179                                  its , ite , jts , jte , kts , kte )
3180
3181       IMPLICIT NONE
3182
3183       INTEGER , INTENT(IN) :: num_soil_layers , &
3184                               num_st_levels_input , num_sm_levels_input , &
3185                               num_st_levels_alloc , num_sm_levels_alloc , &
3186                               ids , ide , jds , jde , kds , kde , &
3187                               ims , ime , jms , jme , kms , kme , &
3188                               its , ite , jts , jte , kts , kte
3189
3190       INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
3191
3192       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
3193       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
3194
3195       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
3196       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
3197       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
3198
3199       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
3200       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
3201       REAL , DIMENSION(num_soil_layers) :: zs , dzs
3202
3203       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
3204
3205       REAL , ALLOCATABLE , DIMENSION(:) :: zhave
3206
3207       INTEGER :: i , j , l , lout , lin , lwant , lhave
3208       REAL :: temp
3209
3210       !  Allocate the soil layer array used for interpolating.
3211
3212       IF ( ( num_st_levels_input .LE. 0 ) .OR. &
3213            ( num_sm_levels_input .LE. 0 ) ) THEN
3214          PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing).  Required for RUC LSM.'
3215          CALL wrf_error_fatal ( 'no soil data' )
3216       ELSE
3217          IF ( flag_soilt000 .eq. 1 ) THEN
3218            PRINT '(A)',' Assume RUC LSM 6-level input'
3219            ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input)  ) )
3220          ELSE
3221            PRINT '(A)',' Assume non-RUC LSM input'
3222            ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers)  ) )
3223          END IF
3224       END IF
3225
3226       !  Sort the levels for temperature.
3227
3228       outert : DO lout = 1 , num_st_levels_input-1
3229          innert : DO lin = lout+1 , num_st_levels_input
3230             IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
3231                temp = st_levels_input(lout)
3232                st_levels_input(lout) = st_levels_input(lin)
3233                st_levels_input(lin) = NINT(temp)
3234                DO j = jts , MIN(jde-1,jte)
3235                   DO i = its , MIN(ide-1,ite)
3236                      temp = st_input(i,lout,j)
3237                      st_input(i,lout,j) = st_input(i,lin,j)
3238                      st_input(i,lin,j) = temp
3239                   END DO
3240                END DO
3241             END IF
3242          END DO innert
3243       END DO outert
3244
3245       IF ( flag_soilt000 .NE. 1 ) THEN
3246       DO j = jts , MIN(jde-1,jte)
3247          DO i = its , MIN(ide-1,ite)
3248             st_input(i,1,j) = tsk(i,j)
3249             st_input(i,num_st_levels_input+2,j) = tmn(i,j)
3250          END DO
3251       END DO
3252       END IF
3253
3254       !  Sort the levels for moisture.
3255
3256       outerm: DO lout = 1 , num_sm_levels_input-1
3257          innerm : DO lin = lout+1 , num_sm_levels_input
3258             IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
3259                temp = sm_levels_input(lout)
3260                sm_levels_input(lout) = sm_levels_input(lin)
3261                sm_levels_input(lin) = NINT(temp)
3262                DO j = jts , MIN(jde-1,jte)
3263                   DO i = its , MIN(ide-1,ite)
3264                      temp = sm_input(i,lout,j)
3265                      sm_input(i,lout,j) = sm_input(i,lin,j)
3266                      sm_input(i,lin,j) = temp
3267                   END DO
3268                END DO
3269             END IF
3270          END DO innerm
3271       END DO outerm
3272
3273       IF ( flag_soilm000 .NE. 1 ) THEN
3274       DO j = jts , MIN(jde-1,jte)
3275          DO i = its , MIN(ide-1,ite)
3276             sm_input(i,1,j) = sm_input(i,2,j)
3277             sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
3278          END DO
3279       END DO
3280       END IF
3281
3282       !  Here are the levels that we have from the input for temperature.
3283
3284       IF ( flag_soilt000 .EQ. 1 ) THEN
3285          DO l = 1 , num_st_levels_input
3286             zhave(l) = st_levels_input(l) / 100.
3287          END DO
3288
3289       !  Interpolate between the layers we have (zhave) and those that we want (zs).
3290
3291       z_wantt : DO lwant = 1 , num_soil_layers
3292          z_havet : DO lhave = 1 , num_st_levels_input -1
3293             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3294                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3295                DO j = jts , MIN(jde-1,jte)
3296                   DO i = its , MIN(ide-1,ite)
3297                      tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
3298                                         st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3299                                                                 ( zhave(lhave+1) - zhave(lhave) )
3300                   END DO
3301                END DO
3302                EXIT z_havet
3303             END IF
3304          END DO z_havet
3305       END DO z_wantt
3306
3307       ELSE
3308
3309          zhave(1) = 0.
3310          DO l = 1 , num_st_levels_input
3311             zhave(l+1) = st_levels_input(l) / 100.
3312          END DO
3313          zhave(num_st_levels_input+2) = 300. / 100.
3314
3315       !  Interpolate between the layers we have (zhave) and those that we want (zs).
3316
3317       z_wantt_2 : DO lwant = 1 , num_soil_layers
3318          z_havet_2 : DO lhave = 1 , num_st_levels_input +2
3319             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3320                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3321                DO j = jts , MIN(jde-1,jte)
3322                   DO i = its , MIN(ide-1,ite)
3323                      tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
3324                                         st_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3325                                                                 ( zhave(lhave+1) - zhave(lhave) )
3326                   END DO
3327                END DO
3328                EXIT z_havet_2
3329             END IF
3330          END DO z_havet_2
3331       END DO z_wantt_2
3332
3333       END IF
3334
3335       !  Here are the levels that we have from the input for moisture.
3336
3337       IF ( flag_soilm000 .EQ. 1 ) THEN
3338          DO l = 1 , num_sm_levels_input
3339             zhave(l) = sm_levels_input(l) / 100.
3340          END DO
3341
3342       !  Interpolate between the layers we have (zhave) and those that we want (zs).
3343
3344       z_wantm : DO lwant = 1 , num_soil_layers
3345          z_havem : DO lhave = 1 , num_sm_levels_input -1
3346             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3347                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3348                DO j = jts , MIN(jde-1,jte)
3349                   DO i = its , MIN(ide-1,ite)
3350                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
3351                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3352                                                                  ( zhave(lhave+1) - zhave(lhave) )
3353                   END DO
3354                END DO
3355                EXIT z_havem
3356             END IF
3357          END DO z_havem
3358       END DO z_wantm
3359
3360       ELSE
3361
3362          zhave(1) = 0.
3363          DO l = 1 , num_sm_levels_input
3364             zhave(l+1) = sm_levels_input(l) / 100.
3365          END DO
3366          zhave(num_sm_levels_input+2) = 300. / 100.
3367
3368       z_wantm_2 : DO lwant = 1 , num_soil_layers
3369          z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
3370             IF ( ( zs(lwant) .GE. zhave(lhave  ) ) .AND. &
3371                  ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
3372                DO j = jts , MIN(jde-1,jte)
3373                   DO i = its , MIN(ide-1,ite)
3374                      smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs   (lwant) ) + &
3375                                          sm_input(i,lhave+1,j) * ( zs   (lwant  ) - zhave(lhave) ) ) / &
3376                                                                  ( zhave(lhave+1) - zhave(lhave) )
3377                   END DO
3378                END DO
3379                EXIT z_havem_2
3380             END IF
3381          END DO z_havem_2
3382       END DO z_wantm_2
3383
3384       END IF
3385       !  Over water, put in reasonable values for soil temperature and moisture.  These won't be
3386       !  used, but they will make a more continuous plot.
3387
3388       IF ( flag_sst .EQ. 1 ) THEN
3389          DO j = jts , MIN(jde-1,jte)
3390             DO i = its , MIN(ide-1,ite)
3391                IF ( landmask(i,j) .LT. 0.5 ) THEN
3392                   DO l = 1 , num_soil_layers
3393                      tslb(i,l,j) = sst(i,j)
3394                      tsk(i,j)    = sst(i,j)
3395                      smois(i,l,j)= 1.0
3396                   END DO
3397                END IF
3398             END DO
3399          END DO
3400       ELSE
3401          DO j = jts , MIN(jde-1,jte)
3402             DO i = its , MIN(ide-1,ite)
3403                IF ( landmask(i,j) .LT. 0.5 ) THEN
3404                   DO l = 1 , num_soil_layers
3405                      tslb(i,l,j)= tsk(i,j)
3406                      smois(i,l,j)= 1.0
3407                   END DO
3408                END IF
3409             END DO
3410          END DO
3411       END IF
3412
3413       DEALLOCATE (zhave)
3414
3415    END SUBROUTINE init_soil_3_real
3416
3417 END MODULE module_soil_pre
3418
3419 #endif
3420