da_detsurtyp.inc

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