da_detsurtyp.inc

References to this file elsewhere.
1 subroutine da_detsurtyp (snow,ice,landsea, vegtyp, soiltyp, & 
2                          is, ie, js, je, &
3                          i, j, dx, dy, dxm, dym, &
4                          isflg,ob_vegtyp,ob_soiltyp, &
5                          seap, icep, lndp, snop)
6 
7    !---------------------------------------------------------------------------
8    ! Purpose: determine surface type at obs locations.
9    !
10    ! METHOD: using the background information 
11    !      1. get the background landmask/snow/sea-ice
12    !      2. calculate percentage of sea/ice/land/snow
13    !      3. determine surface type at obs location
14    !      4. find nearest grid vegtype and soil type 
15    !---------------------------------------------------------------------------
16 
17    implicit none
18 
19    integer, intent(in)  :: is, ie, js, je
20    integer, intent(in)  :: i, j
21    real   , intent(in)  :: dx, dxm, dy, dym
22    real   , intent(in)  :: snow(is:ie,js:je)    ! Input variable
23    real   , intent(in)  :: ice(is:ie,js:je)     ! Input variable
24    real   , intent(in)  :: landsea(is:ie,js:je) ! Input variable
25    integer, intent(in)  :: vegtyp(is:ie,js:je)  ! Input variable
26    integer, intent(in)  :: soiltyp(is:ie,js:je) ! Input variable
27    integer, intent(out) :: isflg                ! Output variable 
28    real,    intent(out) :: ob_vegtyp            ! Output variable
29    real,    intent(out) :: ob_soiltyp            ! Output variable
30    real,    intent(out) :: seap, icep, lndp, snop ! percentage of surface type
31    
32    !     isflg    - surface flag
33    !                0 sea
34    !                1 sea ice
35    !                2 land
36    !                3 snow
37    !                4 mixed predominately sea
38    !                5 mixed predominately sea ice
39    !                6 mixed predominately land
40    !                7 mixed predominately snow
41 
42    !  local variables
43    integer   ::  n, xbflag(4)   ! surface type at xb location
44                                ! 0:sea  1:sea-ice 2:land  3:snow
45    real      ::  w(4),minw      ! weight at 4 xb locations
46 
47    ! 1.0 determine surface type of xb at 4 location around obs
48    !-------------------------------------------------------
49    if ( nint(landsea(i  ,j  )) == 0 ) xbflag(1) = 0 ! sea
50    if ( nint(landsea(i+1,j  )) == 0 ) xbflag(2) = 0 ! sea
51    if ( nint(landsea(i  ,j+1)) == 0 ) xbflag(3) = 0 ! sea
52    if ( nint(landsea(i+1,j+1)) == 0 ) xbflag(4) = 0 ! sea
53 
54    if ( nint(landsea(i  ,j  )) == 1 ) xbflag(1) = 2 ! land/snow/sea-ice
55    if ( nint(landsea(i+1,j  )) == 1 ) xbflag(2) = 2 ! land/snow/sea-ice
56    if ( nint(landsea(i  ,j+1)) == 1 ) xbflag(3) = 2 ! land/snow/sea-ice
57    if ( nint(landsea(i+1,j+1)) == 1 ) xbflag(4) = 2 ! land/snow/sea-ice
58 
59    if ( nint(snow(i  ,j  )) == 1 ) xbflag(1) = 3 ! snow
60    if ( nint(snow(i+1,j  )) == 1 ) xbflag(2) = 3 ! snow
61    if ( nint(snow(i  ,j+1)) == 1 ) xbflag(3) = 3 ! snow
62    if ( nint(snow(i+1,j+1)) == 1 ) xbflag(4) = 3 ! snow
63 
64    if ( nint(ice(i  ,j  )) == 1 ) xbflag(1) = 1 ! sea-ice
65    if ( nint(ice(i+1,j  )) == 1 ) xbflag(2) = 1 ! sea-ice
66    if ( nint(ice(i  ,j+1)) == 1 ) xbflag(3) = 1 ! sea-ice
67    if ( nint(ice(i+1,j+1)) == 1 ) xbflag(4) = 1 ! sea-ice 
68 
69    ! 2.0 determine surface type percentage at obs location
70    !------------------------------------------------------
71    !  (i,j+1) 
72    !    -----------------(i+1,j+1)
73    !    |   w2          |
74    !    |dym        w1  |
75    !    |      obs      |
76    !    |------- *      |
77    !    |dy  w4  |  w3  |
78    !    |   dx   | dxm  |
79    !    |----------------
80    !   (i,j)            (i+1,j)
81    !
82    !--------------------------------------------------------
83 
84    seap = 0.
85    icep = 0.
86    snop = 0.
87    lndp = 0.
88    w(1) = dym*dxm   ! weight for point (i,j)
89    w(2) = dym*dx    ! weight for point (i+1,j)
90    w(3) = dy *dxm   ! weight for point (i,j+1)
91    w(4) = dy *dx    ! weight for point (i+1,j+1)
92 
93    do n = 1, 4
94       if (xbflag(n) == 0) seap = seap+w(n)
95       if (xbflag(n) == 1) icep = icep+w(n)
96       if (xbflag(n) == 2) lndp = lndp+w(n)
97       if (xbflag(n) == 3) snop = snop+w(n)
98    end do
99 
100    ! fo2d   = dym*(dxm*fm2d(i,j  ) + dx*fm2d(i+1,j  )) &
101    !         + dy *(dxm*fm2d(i,j+1) + dx*fm2d(i+1,j+1))
102 
103    ! 3.0 determine final surface flag at obs location
104    !-----------------------------------------
105    if (seap >= 0.99) isflg = 0
106    if (icep >= 0.99) isflg = 1
107    if (lndp >= 0.99) isflg = 2
108    if (snop >= 0.99) isflg = 3
109    if ( .not. (seap >= 0.99) .and. &
110         .not. (icep >= 0.99) .and. &
111         .not. (lndp >= 0.99) .and. &
112         .not. (snop >= 0.99) ) then
113       if (seap > lndp) then
114          if (seap > icep) then
115             if (seap > snop) then
116                isflg = 4
117             else
118                isflg = 7
119             end if
120          else
121             if (icep > snop) then
122                isflg = 5
123             else
124                isflg = 7
125             end if
126          end if
127       else
128          if (lndp > icep) then
129             if (lndp > snop) then
130                isflg = 6
131             else
132                isflg = 7
133             end if
134          else
135             if (icep > snop) then
136                isflg = 5
137             else
138                isflg = 7
139             end if
140          end if
141       end if 
142    end if
143    
144    ! 4.0 find nearest grid vegtype and soil type
145    !         at obs location
146    !-----------------------------------------
147    minw=min(w(1),w(2),w(3),w(4))
148    if (minw == w(1)) then
149       ob_vegtyp  = float(vegtyp (i+1,j+1))
150       ob_soiltyp = float(soiltyp(i+1,j+1))
151    else if (minw == w(2)) then
152       ob_vegtyp  = float(vegtyp (i,j+1))
153       ob_soiltyp = float(soiltyp(i,j+1))
154    else if (minw == w(3)) then
155       ob_vegtyp  = float(vegtyp (i+1,j))
156       ob_soiltyp = float(soiltyp(i+1,j))
157    else if (minw == w(4)) then
158       ob_vegtyp  = float(vegtyp (i,j))
159       ob_soiltyp = float(soiltyp(i,j))
160    end if
161 
162    return
163 
164 end subroutine da_detsurtyp
165