hydro/CPL/WRF_cpl/module_wrf_HYDRO.F

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !  <brief list of changes to this source file>
   6 !
   7 !  Usage:
   8 !  Parameters: <Specify typical arguments passed>
   9 !  Input Files:
  10 !        <list file names and briefly describe the data they include>
  11 !  Output Files:
  12 !        <list file names and briefly describe the information they include>
  13 !
  14 !  Condition codes:
  15 !        <list exit condition or error codes returned >
  16 !        If appropriate, descriptive troubleshooting instructions or
  17 !        likely causes for failures could be mentioned here with the
  18 !        appropriate error code
  19 !
  20 !  User controllable options: <if applicable>
  21
  22 module module_WRF_HYDRO
  23
  24 #ifdef MPP_LAND
  25     use mpi
  26     use module_mpp_land, only: global_nx, global_ny, decompose_data_real, &
  27                  write_io_real, my_id, mpp_land_bcast_real1, IO_id, &
  28                 mpp_land_bcast_real, mpp_land_bcast_int1, mpp_land_init
  29     use module_CPL_LAND, only: CPL_LAND_INIT, cpl_outdate, HYDRO_COMM_WORLD
  30     use module_hydro_stop, only: HYDRO_stop
  31 #endif
  32     use module_HYDRO_drv, only: HYDRO_ini, HYDRO_exe
  33
  34     use module_rt_data, only:  rt_domain
  35     use module_gw_gw2d_data, only:  gw2d
  36     use module_CPL_LAND, only: CPL_LAND_INIT, cpl_outdate
  37     use config_base, only: nlst
  38     USE module_domain, ONLY : domain, domain_clock_get
  39     USE module_configure, ONLY : grid_config_rec_type
  40     !yw USE module_configure, only : config_flags
  41     USE module_configure, only: model_config_rec
  42
  43
  44     implicit none
  45
  46     !yw   added for check soil moisture and soiltype
  47     integer ::  checkSOIL_flag
  48
  49 #ifndef MPP_LAND
  50     character(len=19) :: cpl_outdate
  51 #endif
  52 !
  53 ! added to consider the adaptive time step from WRF model.
  54     real    :: dtrt_ter0  , dtrt_ch0
  55     integer ::  mm0
  56
  57
  58
  59
  60 CONTAINS
  61
  62 !wrf_cpl_HYDRO will not call the off-line lsm
  63 !ywGW subroutine wrf_cpl_HYDRO(HYDRO_dt,grid, config_flags, its,ite,jts,jte)
  64     subroutine wrf_cpl_HYDRO(HYDRO_dt,grid,its,ite,jts,jte)
  65
  66        implicit none
  67        TYPE ( domain ), INTENT(INOUT) :: grid
  68 !ywGW       TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags
  69        integer its, ite, jts, jte, ij
  70        real :: HYDRO_dt
  71
  72
  73         integer k, ix,jx, mm, nn
  74
  75         integer ::  did
  76
  77         integer ntime
  78
  79         integer :: i,j
  80
  81         integer :: ierr
  82
  83 !output flux and state variable
  84
  85         did = 1
  86
  87
  88         ix = ite - its + 1
  89         jx = jte - jts + 1
  90
  91         if(HYDRO_dt .le. 0) then
  92              write(6,*) "WARNING: HYDRO_dt <= 0 from land input. set it to be 1 seconds."
  93              HYDRO_dt = 1
  94         endif
  95
  96         ntime = 1
  97
  98
  99             nlst(did)%dt = HYDRO_dt
 100
 101
 102         if(.not. RT_DOMAIN(did)%initialized) then
 103            !yw nlst_rt(did)%nsoil = config_flags%num_soil_layers
 104            !nlst_rt(did)%nsoil = model_config_rec%num_metgrid_soil_levels
 105            nlst(did)%nsoil = grid%num_soil_layers
 106
 107
 108 #ifdef MPP_LAND
 109            call MPI_COMM_DUP(MPI_COMM_WORLD, HYDRO_COMM_WORLD, ierr)
 110            call MPP_LAND_INIT(grid%e_we - grid%s_we - 1, grid%e_sn - grid%s_sn - 1)
 111
 112            call mpp_land_bcast_int1 (nlst(did)%nsoil)
 113 #endif
 114            allocate(nlst(did)%zsoil8(nlst(did)%nsoil))
 115            if(grid%zs(1) <  0) then
 116               nlst(did)%zsoil8(1:nlst(did)%nsoil) = grid%zs(1:nlst(did)%nsoil)
 117            else
 118               nlst(did)%zsoil8(1:nlst(did)%nsoil) = -1*grid%zs(1:nlst(did)%nsoil)
 119            endif
 120
 121             CALL domain_clock_get( grid, current_timestr=cpl_outdate)
 122             nlst(did)%startdate(1:19) = cpl_outdate(1:19)
 123             nlst(did)%olddate(1:19) = cpl_outdate(1:19)
 124
 125
 126 #ifdef MPP_LAND
 127             call CPL_LAND_INIT(its,ite,jts,jte)
 128 #endif
 129
 130 #ifdef HYDRO_D
 131                write(6,*) "sf_surface_physics is ", grid%sf_surface_physics
 132 #endif
 133
 134            if(grid%sf_surface_physics .eq. 5) then
 135                 ! clm4
 136                call HYDRO_ini(ntime,did=did,ix0=1,jx0=1)
 137            else
 138                call HYDRO_ini(ntime,did,ix0=ix,jx0=jx,vegtyp=grid%IVGTYP(its:ite,jts:jte),soltyp=grid%isltyp(its:ite,jts:jte))
 139            endif
 140
 141
 142
 143             if(nlst(did)%sys_cpl .ne. 2) then
 144                call hydro_stop("In module_wrf_HYDRO.F wrf_cpl_HYDRO() - "// &
 145                                "sys_cpl should be 2.  Check hydro.namelist file.")
 146             endif
 147
 148
 149             nlst(did)%startdate(1:19) = cpl_outdate(1:19)
 150             nlst(did)%olddate(1:19) = cpl_outdate(1:19)
 151
 152             nlst(did)%dt = HYDRO_dt
 153             if(nlst(did)%dtrt_ter .ge. HYDRO_dt) then
 154                nlst(did)%dtrt_ter = HYDRO_dt
 155                mm0 = 1
 156             else
 157                mm = HYDRO_dt/nlst(did)%dtrt_ter
 158                if(mm*nlst(did)%dtrt_ter.lt. HYDRO_dt) nlst(did)%dtrt_ter = HYDRO_dt/mm
 159                mm0 = mm
 160             endif
 161
 162             dtrt_ter0 = nlst(did)%dtrt_ter
 163
 164             if(nlst(did)%dtrt_ch .ge. HYDRO_dt) then
 165                nlst(did)%dtrt_ch = HYDRO_dt
 166                mm0 = 1
 167             else
 168                mm = HYDRO_dt/nlst(did)%dtrt_ch
 169                if(mm*nlst(did)%dtrt_ch.lt. HYDRO_dt) nlst(did)%dtrt_ch = HYDRO_dt/mm
 170                mm0 = mm
 171             endif
 172
 173             dtrt_ch0 = nlst(did)%dtrt_ch
 174         endif
 175
 176             if((mm0*nlst(did)%dtrt_ter) .ne. HYDRO_dt) then   ! WRF model time step changed.
 177                if(dtrt_ter0 .ge. HYDRO_dt) then
 178                   nlst(did)%dtrt_ter = HYDRO_dt
 179                   mm0 = 1
 180                else
 181                   mm = HYDRO_dt/dtrt_ter0
 182                   if(mm*dtrt_ter0 .lt. HYDRO_dt) nlst(did)%dtrt_ter = HYDRO_dt/mm
 183                   mm0 = mm
 184                endif
 185             endif
 186
 187             if((mm0*nlst(did)%dtrt_ch) .ne. HYDRO_dt) then   ! WRF model time step changed.
 188                if(dtrt_ch0 .ge. HYDRO_dt) then
 189                   nlst(did)%dtrt_ch = HYDRO_dt
 190                   mm0 = 1
 191                else
 192                   mm = HYDRO_dt/dtrt_ch0
 193                   if(mm*dtrt_ch0 .lt. HYDRO_dt) nlst(did)%dtrt_ch = HYDRO_dt/mm
 194                   mm0 = mm
 195                endif
 196             endif
 197
 198 #ifdef HYDRO_D
 199         write(6,*) "mm, nlst(did)%dt = ",mm, nlst(did)%dt
 200 #endif
 201
 202         if(nlst(did)%rtFlag .eq. 0) return
 203
 204
 205         nn = nlst(did)%nsoil
 206
 207         ! get the data from WRF
 208
 209
 210
 211        if((.not. RT_DOMAIN(did)%initialized) .and. (nlst(did)%rst_typ .eq. 1) ) then
 212 #ifdef HYDRO_D
 213            write(6,*) "restart initial data from offline file"
 214 #endif
 215        else
 216             do k = 1, nlst(did)%nsoil
 217                 RT_DOMAIN(did)%STC(:,:,k) = grid%TSLB(its:ite,k,jts:jte)
 218                 RT_DOMAIN(did)%smc(:,:,k) = grid%smois(its:ite,k,jts:jte)
 219                 RT_DOMAIN(did)%sh2ox(:,:,k) = grid%sh2o(its:ite,k,jts:jte)
 220             end do
 221             rt_domain(did)%infxsrt = grid%infxsrt(its:ite,jts:jte)
 222             rt_domain(did)%soldrain = grid%soldrain(its:ite,jts:jte)
 223        endif
 224
 225             call HYDRO_exe(did)
 226
 227
 228 ! add for update the WRF state variable.
 229             do k = 1, nlst(did)%nsoil
 230                 ! grid%TSLB(its:ite,k,jts:jte) = RT_DOMAIN(did)%STC(:,:,k)
 231                 grid%smois(its:ite,k,jts:jte) = RT_DOMAIN(did)%smc(:,:,k)
 232                 grid%sh2o(its:ite,k,jts:jte) = RT_DOMAIN(did)%sh2ox(:,:,k)
 233             end do
 234
 235 ! update WRF variable after running routing model.
 236             grid%sfcheadrt(its:ite,jts:jte) = rt_domain(did)%overland%control%surface_water_head_lsm
 237
 238 ! provide groundwater soil flux to WRF for fully coupled simulations (FERSCH 09/2014)
 239             if(nlst(did)%GWBASESWCRT .eq. 3 ) then
 240 !Wei Yu: comment the following two lines. Not ready for WRF3.7 release
 241 !yw             grid%qsgw(its:ite,jts:jte) = gw2d(did)%qsgw
 242 !yw             config_flags%gwsoilcpl = nlst_rt(did)%gwsoilcpl
 243             end if
 244
 245 !yw not sure for the following
 246 !           grid%xice(its:ite,jts:jte) = rt_domain(did)%sice
 247
 248             RT_DOMAIN(did)%initialized = .true.
 249      end subroutine wrf_cpl_HYDRO
 250
 251
 252
 253
 254
 255 !program drive rtland
 256 ! This subroutine will be used if the 4-layer Noah lsm is not used.
 257       subroutine wrf2lsm (z1,v1,kk1,z,vout,ix,jx,kk,vegtyp)
 258 !  input: z1,v1,kk1,z,ix,jx,kk
 259 !  output: vout
 260 !  interpolate based on soil layer: z1 and z
 261 !  z :  soil layer of output variable.
 262 !  z1: array of soil layers of input variable.
 263          implicit none
 264          integer:: i,j,k
 265          integer:: kk1, ix,jx,kk, vegtyp(ix,jx)
 266          real :: z1(kk1), z(kk), v1(ix,kk1,jx),vout(ix,jx,kk)
 267
 268
 269          do j = 1, jx
 270             do i = 1, ix
 271                 do k = 1, kk
 272                   call interpLayer(Z1,v1(i,1:kk1,j),kk1,Z(k),vout(i,j,k))
 273                 end do
 274             end do
 275          end do
 276       end subroutine wrf2lsm
 277
 278 ! This subroutine will be used if the 4-layer Noah lsm is not used.
 279       subroutine lsm2wrf (z1,v1,kk1,z,vout,ix,jx,kk,vegtyp)
 280 !  input: z1,v1,kk1,z,ix,jx,kk
 281 !  output: vout
 282 !  interpolate based on soil layer: z1 and z
 283 !  z :  soil layer of output variable.
 284 !  z1: array of soil layers of input variable.
 285          implicit none
 286          integer:: i,j,k
 287          integer:: kk1, ix,jx,kk, vegtyp(ix,jx)
 288          real :: z1(kk1), z(kk), v1(ix,jx,kk1),vout(ix,kk,jx)
 289
 290
 291          do j = 1, jx
 292             do i = 1, ix
 293                  do k = 1, kk
 294                     call interpLayer(Z1,v1(i,j,1:kk1),kk1,Z(k),vout(i,k,j))
 295                  end do
 296             end do
 297          end do
 298       end subroutine lsm2wrf
 299
 300       subroutine interpLayer(inZ,inV,inK,outZ,outV)
 301          implicit none
 302          integer:: k, k1, k2
 303          integer :: inK
 304          real:: inV(inK),inZ(inK)
 305          real:: outV, outZ, w1, w2
 306
 307          if(outZ .le. inZ(1)) then
 308              w1 = (inZ(2)-outZ)/(inZ(2)-inZ(1))
 309              w2 = (inZ(1)-outZ)/(inZ(2)-inZ(1))
 310              outV = inV(1)*w1-inV(2)*w2
 311              return
 312          elseif(outZ .ge. inZ(inK)) then
 313              w1 = (outZ-inZ(inK-1))/(inZ(inK)-inZ(inK-1))
 314              w2 = (outZ-inZ(inK))  /(inZ(inK)-inZ(inK-1))
 315              outV = inV(inK)*w1 -inV(inK-1)* w2
 316              return
 317          else
 318             do k = 2, inK
 319              if((inZ(k) .ge. outZ).and.(inZ(k-1) .le. outZ) ) then
 320                 k1  = k-1
 321                 k2 = k
 322                 w1 = (outZ-inZ(k1))/(inZ(k2)-inZ(k1))
 323                 w2 = (inZ(k2)-outZ)/(inZ(k2)-inZ(k1))
 324                 outV = inV(k2)*w1 + inV(k1)*w2
 325                 return
 326              end if
 327             end do
 328          endif
 329       end subroutine interpLayer
 330
 331       subroutine lsm_wrf_input(did,vegtyp,soltyp,ix,jx)
 332          implicit none
 333          integer did, leng
 334          parameter(leng=100)
 335          integer :: i,j, nn, ix,jx
 336          integer, dimension(ix,jx) :: soltyp, vegtyp
 337          real, dimension(leng) :: xdum1, MAXSMC,refsmc,wltsmc
 338
 339
 340          where(soltyp == 14) VEGTYP = 16
 341          where(VEGTYP == 16 ) soltyp = 14
 342
 343          RT_DOMAIN(did)%VEGTYP = vegtyp
 344
 345 !      input OV_ROUGH from OVROUGH.TBL
 346 #ifdef MPP_LAND
 347        if(my_id .eq. IO_id) then
 348 #endif
 349
 350 #ifndef NCEP_WCOSS
 351        open(71,file="HYDRO.TBL", form="formatted")
 352 !read OV_ROUGH first
 353           read(71,*) nn
 354           read(71,*)
 355           do i = 1, nn
 356              read(71,*) RT_DOMAIN(did)%OV_ROUGH(i)
 357           end do
 358 !read parameter for LKSAT
 359           read(71,*) nn
 360           read(71,*)
 361           do i = 1, nn
 362              read(71,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i)
 363           end do
 364        close(71)
 365
 366 #else
 367       open(13, form="formatted")
 368 !read OV_ROUGH first
 369           read(13,*) nn
 370           read(13,*)
 371           do i = 1, nn
 372              read(13,*) RT_DOMAIN(did)%OV_ROUGH(i)
 373           end do
 374 !read parameter for LKSAT
 375           read(13,*) nn
 376           read(13,*)
 377           do i = 1, nn
 378              read(13,*) xdum1(i), MAXSMC(i),refsmc(i),wltsmc(i)
 379           end do
 380        close(13)
 381 #endif
 382 #ifdef MPP_LAND
 383        endif
 384        call mpp_land_bcast_real(leng,RT_DOMAIN(did)%OV_ROUGH)
 385        call mpp_land_bcast_real(leng,xdum1)
 386        call mpp_land_bcast_real(leng,MAXSMC)
 387        call mpp_land_bcast_real(leng,refsmc)
 388        call mpp_land_bcast_real(leng,wltsmc)
 389 #endif
 390
 391        rt_domain(did)%lksat = 0.0
 392        do j = 1, RT_DOMAIN(did)%jx
 393              do i = 1, RT_DOMAIN(did)%ix
 394                 rt_domain(did)%lksat(i,j) = xdum1(soltyp(i,j) ) * 1000.0
 395                 IF(rt_domain(did)%VEGTYP(i,j) == 1 ) THEN   ! urban
 396                     rt_domain(did)%SMCMAX1(i,j) = 0.45
 397                     rt_domain(did)%SMCREF1(i,j) = 0.42
 398                     rt_domain(did)%SMCWLT1(i,j) = 0.40
 399                 else
 400                     rt_domain(did)%SMCMAX1(i,j) = MAXSMC(soltyp(I,J))
 401                     rt_domain(did)%SMCREF1(i,j) = refsmc(soltyp(I,J))
 402                     rt_domain(did)%SMCWLT1(i,j) = wltsmc(soltyp(I,J))
 403                 ENDIF
 404              end do
 405        end do
 406
 407
 408       end subroutine lsm_wrf_input
 409
 410       subroutine  checkSoil(did)
 411           implicit none
 412           integer :: did
 413           where(rt_domain(did)%smc(:,:,1) <=0) RT_DOMAIN(did)%VEGTYP = 16
 414           where(rt_domain(did)%sh2ox(:,:,1) <=0) RT_DOMAIN(did)%VEGTYP = 16
 415           where(rt_domain(did)%smc(:,:,1) >=100) RT_DOMAIN(did)%VEGTYP = 16
 416           where(rt_domain(did)%sh2ox(:,:,1) >=100) RT_DOMAIN(did)%VEGTYP = 16
 417       end subroutine checkSoil
 418
 419 end module module_wrf_HYDRO