da_analysis_stats.inc

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