da_analysis_stats.inc

References to this file elsewhere.
1 subroutine da_analysis_stats  (stats_unit, xp, xa)
2    
3    !------------------------------------------------------------------------
4    ! Purpose: Calculate min, max, mean and RMS of input 1d field.
5    !------------------------------------------------------------------------
6 
7    implicit none
8 
9    integer,           intent (in) :: stats_unit ! Output unit for stats.
10    type (xpose_type), intent (in) :: xp         ! Dimensions and xpose buffers.
11    type (x_type),     intent (in) :: xa         ! gridded analysis increments.
12    
13    integer :: i, j, k
14    integer :: ij_g, ijk_g   ! ij, ijk for global domain. 
15    integer :: is, ie        ! i range of processor subdomain.
16    integer :: js, je        ! j range of processor subdomain.
17    integer :: ks, ke, kdim  ! k range 
18 
19    real    :: um, vm, tm, pm, qm ! On local domain.
20    real    :: rij_g, rijk_g      ! On global domain.
21 
22    type (maxmin_field_type) :: max_u(xp%kts:xp%kte), max_v(xp%kts:xp%kte), &
23                                max_t(xp%kts:xp%kte), max_p(xp%kts:xp%kte), &
24                                max_q(xp%kts:xp%kte), &
25                                min_u(xp%kts:xp%kte), min_v(xp%kts:xp%kte), &
26                                min_t(xp%kts:xp%kte), min_p(xp%kts:xp%kte), &
27                                min_q(xp%kts:xp%kte)
28  
29    real                        uv(xp%kts:xp%kte), vv(xp%kts:xp%kte), &
30                                tv(xp%kts:xp%kte), pv(xp%kts:xp%kte), &
31                                qv(xp%kts:xp%kte)
32 
33    call da_trace_entry("da_analysis_stats")
34 
35    is = xp%its
36    ie = xp%ite
37    js = xp%jts
38    je = xp%jte
39    ks = xp%kts
40    ke = xp%kte
41 
42    kdim = ke-ks+1
43 
44    ij_g = mix * mjy
45    ijk_g = ij_g * mkz
46    
47    rij_g  = 1.0/real(ij_g)
48    rijk_g = 1.0/real(ijk_g)
49   
50    if (rootproc) then
51       write(unit=stats_unit, fmt='(/a/)') &
52          ' Minimum of gridded analysis increments'
53 
54       write(unit=stats_unit, fmt='(6a/)') &
55          ' Lvl         ', &
56          'u     i    j          ', &
57          'v     i    j          ', &
58          't     i    j          ', &
59          'p     i    j          ', &
60          'q     i    j          '
61    end if
62 
63    call da_maxmin_in_field(xa%u(is:ie,js:je,ks:ke), max_u, min_u)
64    call da_proc_maxmin_combine(kdim, max_u, min_u)
65    call da_maxmin_in_field(xa%v(is:ie,js:je,ks:ke), max_v, min_v)
66    call da_proc_maxmin_combine(kdim, max_v, min_v)
67    call da_maxmin_in_field(xa%t(is:ie,js:je,ks:ke), max_t, min_t)
68    call da_proc_maxmin_combine(kdim, max_t, min_t)
69    call da_maxmin_in_field(xa%p(is:ie,js:je,ks:ke), max_p, min_p)
70    call da_proc_maxmin_combine(kdim, max_p, min_p)
71    call da_maxmin_in_field(xa%q(is:ie,js:je,ks:ke), max_q, min_q)
72    call da_proc_maxmin_combine(kdim, max_q, min_q)
73 
74    um = 999999.0
75    vm = 999999.0
76    tm = 999999.0
77    pm = 999999.0
78    qm = 999999.0
79 
80    do k = ks, ke   
81       if (rootproc) then
82          write(unit=stats_unit, fmt='(i4,4(f12.4,2i5),e12.4,2i5)') k, &
83             min_u(k), min_v(k), min_t(k), min_p(k), min_q(k)
84       end if
85 
86       um=minval(min_u(:)%value)
87       vm=minval(min_v(:)%value)
88       tm=minval(min_t(:)%value)
89       pm=minval(min_p(:)%value)
90       qm=minval(min_q(:)%value)
91    end do
92 
93    if (rootproc) then
94       write(unit=stats_unit, fmt='(a,4(f12.4,10x),e12.4)') ' ALL', &
95          um, vm, tm, pm, qm
96    
97   
98       write(unit=stats_unit, fmt='(/a/)') &
99          ' Maximum of gridded analysis increments'
100 
101       write(unit=stats_unit, fmt='(6a/)') &
102          ' Lvl         ', &
103          'u     i    j          ', &
104          'v     i    j          ', &
105          't     i    j          ', &
106          'p     i    j          ', &
107          'q     i    j          '
108    end if
109 
110    do k = ks, ke
111       if (rootproc) then
112          write(unit=stats_unit, fmt='(i4,4(f12.4,2i5),e12.4,2i5)') k, &
113             max_u(k), max_v(k), max_t(k), max_p(k), max_q(k)
114       end if
115 
116       um=maxval(max_u(:)%value)
117       vm=maxval(max_v(:)%value)
118       tm=maxval(max_t(:)%value)
119       pm=maxval(max_p(:)%value)
120       qm=maxval(max_q(:)%value)
121    end do
122 
123    if (rootproc) then
124       write(unit=stats_unit, fmt='(a,4(f12.4,10x),e12.4)') ' ALL', &
125          um, vm, tm, pm, qm
126    
127       write(unit=stats_unit, fmt='(/a/)') ' Mean of gridded analysis increments'
128 
129       write(unit=stats_unit, fmt='(a/)') &
130          ' Lvl        u           v           t           p           q'
131    end if
132 
133    um = 0.0
134    vm = 0.0
135    tm = 0.0
136    pm = 0.0
137    qm = 0.0
138 
139    do k = ks, ke
140       uv(k) = sum(xa%u(is:ie,js:je,k))
141       vv(k) = sum(xa%v(is:ie,js:je,k))
142       tv(k) = sum(xa%t(is:ie,js:je,k))
143       pv(k) = sum(xa%p(is:ie,js:je,k))
144       qv(k) = sum(xa%q(is:ie,js:je,k))
145    end do
146 
147    call da_proc_sum_real  (uv)
148    call da_proc_sum_real  (vv)
149    call da_proc_sum_real  (tv)
150    call da_proc_sum_real  (pv)
151    call da_proc_sum_real  (qv)
152 
153    if (rootproc) then
154       do k = ks, ke
155          write(unit=stats_unit, fmt='(i4,4f12.4,4e12.4)') k, &
156             uv(k)*rij_g, vv(k)*rij_g, tv(k)*rij_g, &
157             pv(k)*rij_g, qv(k)*rij_g
158 
159          um=um+uv(k)
160          vm=vm+vv(k)
161          tm=tm+tv(k)
162          pm=pm+pv(k)
163          qm=qm+qv(k)
164       end do
165    end if
166 
167    if (rootproc) then
168       write(unit=stats_unit, fmt='(a,4f12.4,4e12.4)') ' ALL', &
169          um*rijk_g, vm*rijk_g, tm*rijk_g, pm*rijk_g, qm*rijk_g
170 
171       write(unit=stats_unit, fmt='(/a/)') ' RMSE of gridded analysis increments'
172 
173       write(unit=stats_unit, fmt='(a/)') &
174          ' Lvl        u           v           t           p           q'
175    end if
176 
177    um = 0.0
178    vm = 0.0
179    tm = 0.0
180    pm = 0.0
181    qm = 0.0
182    uv = 0.0
183    vv = 0.0
184    tv = 0.0
185    pv = 0.0
186    qv = 0.0
187 
188    do k = ks, ke
189       do j=js,je
190          do i=is,ie
191             uv(k) = uv(k) + xa%u(i,j,k) * xa%u(i,j,k)
192             vv(k) = vv(k) + xa%v(i,j,k) * xa%v(i,j,k)
193             tv(k) = tv(k) + xa%t(i,j,k) * xa%t(i,j,k)
194             pv(k) = pv(k) + xa%p(i,j,k) * xa%p(i,j,k)
195             qv(k) = qv(k) + xa%q(i,j,k) * xa%q(i,j,k)
196          end do
197       end do
198    end do
199 
200    call da_proc_sum_real  (uv)
201    call da_proc_sum_real  (vv)
202    call da_proc_sum_real  (tv)
203    call da_proc_sum_real  (pv)
204    call da_proc_sum_real  (qv)
205 
206    if (rootproc) then
207       do k = ks, ke
208          write(unit=stats_unit, fmt='(i4,4f12.4,4e12.4)') k, &
209             sqrt(uv(k)*rij_g), &
210             sqrt(vv(k)*rij_g), &
211             sqrt(tv(k)*rij_g), &
212             sqrt(pv(k)*rij_g), &
213             sqrt(qv(k)*rij_g)
214 
215          um=um+uv(k)
216          vm=vm+vv(k)
217          tm=tm+tv(k)
218          pm=pm+pv(k)
219          qm=qm+qv(k)
220       end do
221    end if
222 
223    if (rootproc) then
224       write(unit=stats_unit, fmt='(a,4f12.4,4e12.4)') ' ALL', &
225          sqrt(um*rijk_g), sqrt(vm*rijk_g), sqrt(tm*rijk_g), &
226          sqrt(pm*rijk_g), sqrt(qm*rijk_g)
227    end if
228 
229    call da_trace_exit("da_analysis_stats")
230 
231 contains
232 
233 #include "da_maxmin_in_field.inc"
234    
235 end subroutine da_analysis_stats
236 
237