hydro/Routing/module_HYDRO_utils.F

   1 !  Program Name:
   2 !  Author(s)/Contact(s):
   3 !  Abstract:
   4 !  History Log:
   5 !
   6 !  Usage:
   7 !  Parameters: <Specify typical arguments passed>
   8 !  Input Files:
   9 !        <list file names and briefly describe the data they include>
  10 !  Output Files:
  11 !        <list file names and briefly describe the information they include>
  12 !
  13 !  Condition codes:
  14 !        <list exit condition or error codes returned >
  15 !        If appropriate, descriptive troubleshooting instructions or
  16 !        likely causes for failures could be mentioned here with the
  17 !        appropriate error code
  18 !
  19 !  User controllable options: <if applicable>
  20
  21 module module_HYDRO_utils
  22   use module_RT_data, only: rt_domain
  23   use config_base, only: nlst
  24   use iso_fortran_env, only: int64
  25 #ifdef MPP_LAND
  26   use module_mpp_land, only: global_nx, global_ny, my_id, IO_id, &
  27        decompose_data_real, write_io_real, MPP_LAND_COM_REAL, &
  28        write_io_int, mpp_land_bcast_real, global_rt_nx, global_rt_ny, &
  29        decompose_rt_real, write_io_rt_real
  30   use MODULE_mpp_GWBUCKET, only: gw_decompose_real
  31 #endif
  32
  33
  34   implicit none
  35   logical lr_dist_flag    !land routing distance calculated or not.
  36
  37 contains
  38
  39   integer function get2d_real(var_name,out_buff,ix,jx,fileName)
  40     implicit none
  41 #         include "netcdf.inc"
  42     integer :: ivar, iret,varid,ncid,ix,jx
  43     real out_buff(ix,jx)
  44     character(len=*), intent(in) :: var_name
  45     character(len=*), intent(in) :: fileName
  46     get2d_real = -1
  47
  48     iret = nf_open(trim(fileName), NF_NOWRITE, ncid)
  49     if (iret .ne. 0) then
  50 #ifdef HYDRO_D
  51        print*,"Failed to open the netcdf file: ",trim(fileName)
  52 #endif
  53        out_buff = -9999.
  54        return
  55     endif
  56     ivar = nf_inq_varid(ncid,trim(var_name),  varid)
  57     if(ivar .ne. 0) then
  58        ivar = nf_inq_varid(ncid,trim(var_name//"_M"),  varid)
  59        if(ivar .ne. 0) then
  60 #ifdef HYDRO_D
  61           write(6,*) "Read Error: could not find ",var_name
  62 #endif
  63           return
  64        endif
  65     end if
  66     iret = nf_get_var_real(ncid, varid, out_buff)
  67     iret = nf_close(ncid)
  68     get2d_real =  ivar
  69   end function get2d_real
  70
  71
  72 ! this module create the distance dx, dy and diagnoal for routing
  73 ! 8 direction as the slop:
  74 ! 1: i,j+1
  75 ! 2: i+1, j+1
  76 ! 3: i+1, j
  77 ! 4: i+1, j-1
  78 ! 5: i, j-1
  79 ! 6: i-1, j-1
  80 ! 7: i-1, j
  81 ! 8: i-1, j+1
  82    real function get_dy(i,j,v,ix,jx)
  83       ! south north
  84        integer :: i,j,ix,jx
  85        real,dimension(ix,jx,9) :: v
  86        if( v(i,j,1) .le. 0) then
  87           get_dy = v(i,j,5)
  88        else if( v(i,j,5) .le. 0) then
  89           get_dy = v(i,j,1)
  90        else
  91           get_dy = (v(i,j,1) + v(i,j,5) ) / 2
  92        endif
  93        return
  94    end function get_dy
  95
  96    real function get_dx(i,j,v,ix,jx)
  97       ! east-west
  98        integer :: i,j, ix,jx
  99        real,dimension(ix,jx,9) :: v
 100        if( v(i,j,3) .le. 0) then
 101           get_dx = v(i,j,7)
 102        else if( v(i,j,7) .le. 0) then
 103           get_dx = v(i,j,3)
 104        else
 105           get_dx = (v(i,j,3) + v(i,j,7) ) / 2
 106        endif
 107        return
 108    end function get_dx
 109
 110    real function get_ll_d(lat1_in, lat2_in, lon1_in, lon2_in)
 111      implicit none
 112      real:: lat1, lat2, lon1, lon2
 113      real:: lat1_in, lat2_in, lon1_in, lon2_in
 114      real::  r, pai, a,c, dlat, dlon, b1,b2
 115      pai = 3.14159
 116      lat1 = lat1_in * pai/180
 117      lat2 = lat2_in * pai/180
 118      lon1 = lon1_in * pai/180
 119      lon2 = lon2_in * pai/180
 120      r = 6378.1*1000
 121      dlat = lat2 -lat1
 122      dlon = lon2 -lon1
 123      a = sin(dlat/2)*sin(dlat/2) + cos(lat1)*cos(lat2)*sin(dlon/2)*sin(dlon/2)
 124      b1 = sqrt(a)
 125      b2 = sqrt(1-a)
 126      c = 2.0*atan2(b1,b2)
 127      get_ll_d = R*c
 128      return
 129
 130    end function get_ll_d
 131
 132    real function get_ll_d_tmp(lat1_in, lat2_in, lon1_in, lon2_in)
 133      implicit none
 134      real:: lat1, lat2, lon1, lon2
 135      real:: lat1_in, lat2_in, lon1_in, lon2_in
 136      real::  r, pai
 137      pai = 3.14159
 138      lat1 = lat1_in * pai/180
 139      lat2 = lat2_in * pai/180
 140      lon1 = lon1_in * pai/180
 141      lon2 = lon2_in * pai/180
 142      r = 6371*1000
 143      get_ll_d_tmp = acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon2-lon1))*r
 144      return
 145
 146    end function get_ll_d_tmp
 147
 148    subroutine get_rt_dxdy_ll(did)
 149 !   use the land lat and lon to derive the routing distrt
 150       implicit none
 151       integer:: did, k
 152       integer iret
 153 !     external get2d_real
 154 !     real get2d_real
 155 #ifdef MPP_LAND
 156       real, dimension(global_rt_nx,global_rt_ny):: latrt, lonrt
 157       real, dimension(global_rt_nx,global_rt_ny,9):: dist
 158       if(my_id .eq. IO_id) then
 159  ! read the lat and lon.
 160          iret =  get2d_real("LONGITUDE",lonrt,global_rt_nx,global_rt_ny,&
 161                      trim(nlst(did)%GEO_FINEGRID_FLNM ))
 162          iret =  get2d_real("LATITUDE",latrt,global_rt_nx,global_rt_ny,&
 163                      trim(nlst(did)%GEO_FINEGRID_FLNM ))
 164          call get_dist_ll(dist,latrt,lonrt,global_rt_nx,global_rt_ny)
 165       end if
 166      do k = 1 , 9
 167         call decompose_RT_real(dist(:,:,k),rt_domain(did)%overland%properties%distance_to_neighbor(:,:,k), &
 168                 global_rt_nx,global_rt_ny,rt_domain(did)%ixrt,rt_domain(did)%jxrt)
 169      end do
 170 #else
 171       real, dimension(rt_domain(did)%ixrt,rt_domain(did)%jxrt):: latrt, lonrt
 172  ! read the lat and lon.
 173          iret =  get2d_real("LONGITUDE",lonrt,rt_domain(did)%ixrt,rt_domain(did)%jxrt,&
 174                      trim(nlst(did)%GEO_FINEGRID_FLNM ))
 175          iret =  get2d_real("LATITUDE",latrt,rt_domain(did)%ixrt,rt_domain(did)%jxrt,&
 176                      trim(nlst(did)%GEO_FINEGRID_FLNM ))
 177          call get_dist_ll(rt_domain(did)%overland%properties%distance_to_neighbor,latrt,lonrt,rt_domain(did)%ixrt,rt_domain(did)%jxrt)
 178 #endif
 179
 180    end subroutine get_rt_dxdy_ll
 181
 182 !  get dx and dy of lat and lon
 183    subroutine get_dist_ll(dist,lat,lon,ix,jx)
 184       implicit none
 185       integer:: ix,jx
 186       real, dimension(ix,jx,9):: dist
 187       real, dimension(ix,jx):: lat, lon
 188       integer:: i,j
 189       real x,y
 190       dist = -1
 191       do j = 1, jx
 192         do i = 1, ix
 193           if(j .lt. jx) dist(i,j,1) = &
 194              get_ll_d(lat(i,j), lat(i,j+1), lon(i,j), lon(i,j+1))
 195           if(j .lt. jx .and. i .lt. ix) dist(i,j,2) =  &
 196              get_ll_d(lat(i,j), lat(i+1,j+1), lon(i,j), lon(i+1,j+1))
 197           if(i .lt. ix) dist(i,j,3) = &
 198              get_ll_d(lat(i,j), lat(i+1,j), lon(i,j), lon(i+1,j))
 199           if(j .gt. 1 .and. i .lt. ix) dist(i,j,4) = &
 200              get_ll_d(lat(i,j), lat(i+1,j-1), lon(i,j), lon(i+1,j-1))
 201           if(j .gt. 1 ) dist(i,j,5) = &
 202              get_ll_d(lat(i,j), lat(i,j-1), lon(i,j), lon(i,j-1))
 203           if(j .gt. 1 .and. i .gt. 1) dist(i,j,6) = &
 204              get_ll_d(lat(i,j), lat(i-1,j-1), lon(i,j), lon(i-1,j-1))
 205           if(i .gt. 1) dist(i,j,7) = &
 206              get_ll_d(lat(i,j), lat(i-1,j), lon(i,j), lon(i-1,j))
 207           if(j .lt. jx .and. i .gt. 1) dist(i,j,8) = &
 208              get_ll_d(lat(i,j), lat(i-1,j+1), lon(i,j), lon(i-1,j+1))
 209         end do
 210       end do
 211       do j = 1, jx
 212         do i = 1, ix
 213             if(j.eq.1) then
 214                y =  get_ll_d(lat(i,j), lat(i,j+1), lon(i,j), lon(i,j+1))
 215             else if(j.eq.jx) then
 216                y =  get_ll_d(lat(i,j-1), lat(i,j), lon(i,j-1), lon(i,j))
 217             else
 218                y =  get_ll_d(lat(i,j-1), lat(i,j+1), lon(i,j-1), lon(i,j+1))/2.0
 219             endif
 220
 221             if(i.eq.ix) then
 222                 x =  get_ll_d(lat(i,j), lat(i-1,j), lon(i,j), lon(i-1,j))
 223             else if(i.eq.1) then
 224                 x =  get_ll_d(lat(i,j), lat(i+1,j), lon(i,j), lon(i+1,j))
 225             else
 226                 x =  get_ll_d(lat(i-1,j), lat(i+1,j), lon(i-1,j), lon(i+1,j))/2.0
 227             endif
 228             dist(i,j,9) = x * y
 229         end do
 230       end do
 231 #ifdef HYDRO_D
 232       write(6,*) "finished get_dist_ll"
 233 #endif
 234    end subroutine get_dist_ll
 235
 236 !  get dx and dy of map projected
 237    subroutine get_dxdy_mp(dist,ix,jx,dx,dy)
 238       implicit none
 239       integer:: ix,jx
 240       real :: dx,dy
 241       integer:: i,j
 242       real :: v1
 243       ! out variable
 244       real, dimension(ix,jx,9)::dist
 245       dist = -1
 246       v1 = sqrt(dx*dx + dy*dy)
 247       do j = 1, jx
 248         do i = 1, ix
 249           if(j .lt. jx) dist(i,j,1) = dy
 250           if(j .lt. jx .and. i .lt. ix) dist(i,j,2) = v1
 251           if(i .lt. ix) dist(i,j,3) = dx
 252           if(j .gt. 1 .and. i .lt. ix) dist(i,j,4) = v1
 253           if(j .gt. 1 ) dist(i,j,5) = dy
 254           if(j .gt. 1 .and. i .gt. 1) dist(i,j,6) = v1
 255           if(i .gt. 1) dist(i,j,7) = dx
 256           if(j .lt. jx .and. i .gt. 1) dist(i,j,8) = v1
 257           dist(i,j,9) = dx * dy
 258         end do
 259       end do
 260 #ifdef HYDRO_D
 261       write(6,*) "finished get_dxdy_mp "
 262 #endif
 263    end subroutine get_dxdy_mp
 264
 265    subroutine get_dist_lsm(did)
 266      integer did
 267 #ifdef MPP_LAND
 268      integer ix,jx,ixrt,jxrt, k
 269      real , dimension(global_nx,global_ny):: latitude,longitude
 270      real, dimension(global_nx,global_ny,9):: dist
 271      if(nlst(did)%dxrt0 .lt. 0) then
 272            ! lat and lon grid
 273           call write_io_real(rt_domain(did)%lat_lsm,latitude)
 274           call write_io_real(rt_domain(did)%lon_lsm,longitude)
 275           if(my_id.eq.IO_id) then
 276                call get_dist_ll(dist,latitude,longitude,  &
 277                          global_nx,global_ny)
 278           endif
 279
 280      else
 281            ! mapp projected grid.
 282           if(my_id.eq.IO_id) then
 283               call get_dxdy_mp(dist,global_nx,global_ny, &
 284                  nlst(did)%dxrt0*nlst(did)%AGGFACTRT,nlst(did)%dxrt0*nlst(did)%AGGFACTRT)
 285           endif
 286      endif
 287      do k = 1 , 9
 288         call decompose_data_real(dist(:,:,k),rt_domain(did)%dist_lsm(:,:,k))
 289      end do
 290 #else
 291      if(nlst(did)%dxrt0 .lt. 0) then
 292         ! lat and lon grid
 293         call get_dist_ll(rt_domain(did)%dist_lsm,rt_domain(did)%lat_lsm,rt_domain(did)%lon_lsm,  &
 294                       rt_domain(did)%ix,rt_domain(did)%jx)
 295      else
 296         ! mapp projected grid.
 297         call get_dxdy_mp(rt_domain(did)%dist_lsm,rt_domain(did)%ix,rt_domain(did)%jx, &
 298               nlst(did)%dxrt0*nlst(did)%AGGFACTRT,nlst(did)%dxrt0*nlst(did)%AGGFACTRT)
 299      endif
 300 #endif
 301
 302
 303    end subroutine get_dist_lsm
 304
 305    subroutine get_dist_lrt(did)
 306      integer did, k
 307
 308 !     real :: tmp_dist(global_rt_nx, global_rt_ny,9)
 309
 310 ! calculate the distance for land routing from the lat /lon of land surface model
 311      if(nlst(did)%dxrt0 .lt. 0) then
 312         ! using lat and lon grid when channel routing is off
 313         call get_rt_dxdy_ll(did)
 314      else
 315         ! mapp projected grid.
 316          call get_dxdy_mp(rt_domain(did)%overland%properties%distance_to_neighbor,rt_domain(did)%ixrt,rt_domain(did)%jxrt, &
 317               nlst(did)%dxrt0,nlst(did)%dxrt0)
 318 #ifdef MPP_LAND
 319         do k = 1, 9
 320            call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%distance_to_neighbor(:,:,k),rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
 321         end do
 322 #endif
 323      endif
 324
 325
 326    end subroutine get_dist_lrt
 327
 328 !   subroutine get_dist_crt(did)
 329 !      integer did, k
 330 ! calculate the distance from channel routing
 331 !     if(nlst_rt(did)%dxrt0 .lt. 0) then
 332 !        ! lat and lon grid
 333 !        if(rt_domain(did)%overland%properties%distance_to_neighbor(1,1,9) .eq. -999)   &
 334 !           call get_dist_ll(rt_domain(did)%overland%properties%distance_to_neighbor,rt_domain(did)%latval,rt_domain(did)%lonval,  &
 335 !                      rt_domain(did)%ixrt,rt_domain(did)%jxrt)
 336 !     else
 337 !        ! mapp projected grid.
 338 !        if(rt_domain(did)%overland%properties%distance_to_neighbor(1,1,9) .eq. -999)   &
 339 !           call get_dxdy_mp(rt_domain(did)%overland%properties%distance_to_neighbor,rt_domain(did)%ixrt,rt_domain(did)%jxrt, &
 340 !              nlst_rt(did)%dxrt0,nlst_rt(did)%dxrt0)
 341 !     endif
 342 !#ifdef MPP_LAND
 343 !     do k = 1, 9
 344 !       call MPP_LAND_COM_REAL(rt_domain(did)%overland%properties%distance_to_neighbor(:,:,k),rt_domain(did)%IXRT,rt_domain(did)%JXRT,99)
 345 !     end do
 346 !#endif
 347 !   end subroutine get_dist_crt
 348
 349    subroutine get_basn_area(did)
 350       implicit none
 351       integer :: did, ix,jx, k
 352       real :: basns_area(rt_domain(did)%gnumbasns)
 353 #ifdef MPP_LAND
 354       integer :: mask(global_nx, global_ny)
 355       real :: dist_lsm(global_nx, global_ny,9)
 356 #else
 357       integer :: mask(rt_domain(did)%ix, rt_domain(did)%jx)
 358       real :: dist_lsm(rt_domain(did)%ix, rt_domain(did)%jx,9)
 359 #endif
 360 #ifdef MPP_LAND
 361       ix = global_nx
 362       jx = global_ny
 363       call write_IO_int(rt_domain(did)%GWSUBBASMSK,mask)
 364       do k = 1,  9
 365          call write_IO_real(rt_domain(did)%dist_lsm(:,:,k),dist_lsm(:,:,k))
 366       end do
 367 #else
 368       ix = rt_domain(did)%ix
 369       jx = rt_domain(did)%jx
 370       mask = rt_domain(did)%GWSUBBASMSK
 371       dist_lsm = rt_domain(did)%dist_lsm
 372 #endif
 373
 374 #ifdef MPP_LAND
 375       if(my_id .eq. IO_id) then
 376          call get_area_g(basns_area,mask, rt_domain(did)%gnumbasns,ix,jx,dist_lsm)
 377       end if
 378 !      call mpp_land_bcast_real(rt_domain(did)%numbasns,rt_domain(did)%basns_area)
 379
 380       call gw_decompose_real(rt_domain(did)%gnumbasns, rt_domain(did)%numbasns, &
 381            rt_domain(did)%basnsInd, basns_area,rt_domain(did)%basns_area)
 382 #else
 383       call get_area_g(rt_domain(did)%basns_area,mask, rt_domain(did)%gnumbasns,ix,jx,dist_lsm)
 384 #endif
 385    end subroutine get_basn_area
 386
 387    subroutine get_area_g(basns_area,GWSUBBASMSK, numbasns,ix,jx,dist)
 388       integer :: i,j, n, ix,jx, numbasns
 389       integer :: count(numbasns)
 390       real :: basns_area(numbasns) , dist(ix,jx,9)
 391       integer :: GWSUBBASMSK(ix,jx)
 392       basns_area = 0
 393       count = 0
 394       do  j = 1, jx
 395         do  i = 1, ix
 396            n = GWSUBBASMSK(i,j)
 397            if(n .gt. 0) then
 398               basns_area(n) = basns_area(n)+dist(i,j,9)
 399               count(n) = count(n) + 1
 400            endif
 401         end do
 402       end do
 403       do i = 1, numbasns
 404          if(count(i) .gt. 0) then
 405              basns_area(i) = basns_area(i) / count(i)
 406          end if
 407       end do
 408    end subroutine get_area_g
 409
 410     subroutine get_area_g8(basns_area,GWSUBBASMSK, numbasns,ix,jx,dist)
 411         integer :: i,j, n, ix,jx, numbasns
 412         integer :: count(numbasns)
 413         real :: basns_area(numbasns) , dist(ix,jx,9)
 414         integer(kind=int64) :: GWSUBBASMSK(ix,jx)
 415         basns_area = 0
 416         count = 0
 417         do  j = 1, jx
 418             do  i = 1, ix
 419                 n = GWSUBBASMSK(i,j)
 420                 if(n .gt. 0) then
 421                     basns_area(n) = basns_area(n)+dist(i,j,9)
 422                     count(n) = count(n) + 1
 423                 endif
 424             end do
 425         end do
 426         do i = 1, numbasns
 427             if(count(i) .gt. 0) then
 428                 basns_area(i) = basns_area(i) / count(i)
 429             end if
 430         end do
 431     end subroutine get_area_g8
 432
 433    subroutine get_node_area(did)
 434        integer :: did
 435
 436        call get_area_g8(rt_domain(did)%node_area, rt_domain(did)%CH_NETLNK, &
 437          rt_domain(did)%NLINKS,rt_domain(did)%ixrt, rt_domain(did)%jxrt, rt_domain(did)%overland%properties%distance_to_neighbor)
 438    end subroutine get_node_area
 439
 440
 441 end module module_HYDRO_utils