FV3 Bundle
fv3jedi_increment_mod.f90
Go to the documentation of this file.
1 ! (C) Copyright 2017-2018 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 !> Handle increment for the FV3JEDI model
7 
8 module fv3jedi_increment_mod
9 
10 use iso_c_binding
11 use config_mod
12 use datetime_mod
13 
14 use fv3jedi_constants_mod, only: rad2deg, constoz, cp, alhl, rgas
15 use fv3jedi_geom_mod, only: fv3jedi_geom
16 use fv3jedi_increment_utils_mod, only: fv3jedi_increment
17 use fv3jedi_increment_io_mod
18 use fv3jedi_kinds_mod, only: kind_real
19 use fv3jedi_state_utils_mod, only: fv3jedi_state
20 use fv3jedi_vars_mod, only: fv3jedi_vars
21 use fv3jedi_getvalues_mod, only: getvalues_tl, getvalues_ad
22 
23 implicit none
24 private
25 
26 public :: fv3jedi_increment
27 
28 public :: create, delete, zeros, random, copy, &
29  self_add, self_schur, self_sub, self_mul, axpy_inc, axpy_state, &
30  dot_prod, add_incr, diff_incr, &
31  read_file, write_file, gpnorm, incrms, &
32  change_resol, getvalues_tl, getvalues_ad, &
33  ug_coord, increment_to_ug, increment_from_ug, dirac, jnormgrad
34 
35 ! ------------------------------------------------------------------------------
36 
37 contains
38 
39 ! ------------------------------------------------------------------------------
40 
41 subroutine create(self, geom, vars)
42 
43 implicit none
44 type(fv3jedi_increment), intent(inout) :: self
45 type(fv3jedi_geom), target, intent(in) :: geom
46 type(fv3jedi_vars), intent(in) :: vars
47 
48 integer :: var
49 
50 ! Copy the variable names
51 self%vars%nv = vars%nv
52 allocate(self%vars%fldnames(self%vars%nv))
53 self%vars%fldnames = vars%fldnames
54 
55 ! Allocate variables based on names
56 do var = 1, self%vars%nv
57 
58  select case (trim(self%vars%fldnames(var)))
59 
60  case("ua")
61  if (.not.allocated( self%ua)) allocate ( self%ua(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
62  case("va")
63  if (.not.allocated( self%va)) allocate ( self%va(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
64  case("t")
65  if (.not.allocated( self%t)) allocate ( self%t(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
66  case("ps")
67  if (.not.allocated( self%ps)) allocate ( self%ps(geom%isc:geom%iec, geom%jsc:geom%jec ))
68  case("q")
69  if (.not.allocated( self%q)) allocate ( self%q(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
70  case("qi")
71  if (.not.allocated( self%qi)) allocate ( self%qi(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
72  case("ql")
73  if (.not.allocated( self%ql)) allocate ( self%ql(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
74  case("o3")
75  if (.not.allocated( self%o3)) allocate ( self%o3(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
76  case("psi")
77  if (.not.allocated( self%psi)) allocate ( self%psi(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
78  case("chi")
79  if (.not.allocated( self%chi)) allocate ( self%chi(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
80  case("tv")
81  if (.not.allocated( self%tv)) allocate ( self%tv(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
82  case("qc")
83  if (.not.allocated( self%qc)) allocate ( self%qc(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
84  case("qic")
85  if (.not.allocated( self%qic)) allocate ( self%qic(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
86  case("qlc")
87  if (.not.allocated( self%qlc)) allocate ( self%qlc(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
88  case("o3c")
89  if (.not.allocated( self%o3c)) allocate ( self%o3c(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
90  case("w")
91  if (.not.allocated( self%w)) allocate ( self%w(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
92  case("delz")
93  if (.not.allocated(self%delz)) allocate (self%delz(geom%isc:geom%iec, geom%jsc:geom%jec , geom%npz))
94  case("ud")
95  case("vd")
96  case("delp")
97  case default
98  call abor1_ftn("Increment: unknown variable "//trim(self%vars%fldnames(var)))
99 
100  end select
101 
102 enddo
103 
104 self%hydrostatic = .true.
105 if (allocated(self%w).and.allocated(self%delz).and.allocated(self%delp)) self%hydrostatic = .false.
106 
107 if (allocated(self%ps) .and. allocated(self%delp)) then
108  call abor1_ftn("Increment: Ps and delp are both allocated, only one can be used")
109 endif
110 
111 ! Initialize all domain arrays to zero
112 call zeros(self)
113 
114 ! For convenience
115 self%isc = geom%isc
116 self%iec = geom%iec
117 self%jsc = geom%jsc
118 self%jec = geom%jec
119 self%isd = geom%isd
120 self%ied = geom%ied
121 self%jsd = geom%jsd
122 self%jed = geom%jed
123 self%npx = geom%npx
124 self%npy = geom%npy
125 self%npz = geom%npz
126 
127 end subroutine create
128 
129 ! ------------------------------------------------------------------------------
130 
131 subroutine delete(self)
132 implicit none
133 type(fv3jedi_increment), intent(inout) :: self
134 
135 if (allocated(self%ua )) deallocate (self%ua )
136 if (allocated(self%va )) deallocate (self%va )
137 if (allocated(self%t )) deallocate (self%t )
138 if (allocated(self%ps )) deallocate (self%ps )
139 if (allocated(self%q )) deallocate (self%q )
140 if (allocated(self%qi )) deallocate (self%qi )
141 if (allocated(self%ql )) deallocate (self%ql )
142 if (allocated(self%o3 )) deallocate (self%o3 )
143 
144 if (allocated(self%psi )) deallocate(self%psi )
145 if (allocated(self%chi )) deallocate(self%chi )
146 if (allocated(self%tv )) deallocate(self%tv )
147 if (allocated(self%qc )) deallocate(self%qc )
148 if (allocated(self%qic )) deallocate(self%qic )
149 if (allocated(self%qlc )) deallocate(self%qlc )
150 if (allocated(self%o3c )) deallocate(self%o3c )
151 
152 if (allocated(self%w )) deallocate (self%w )
153 if (allocated(self%delz)) deallocate (self%delz)
154 if (allocated(self%delp)) deallocate (self%delp)
155 
156 end subroutine delete
157 
158 ! ------------------------------------------------------------------------------
159 
160 subroutine zeros(self)
161 implicit none
162 type(fv3jedi_increment), intent(inout) :: self
163 
164 !Zero out the entire domain
165 
166 if(allocated(self%ua )) self%ua = 0.0_kind_real
167 if(allocated(self%va )) self%va = 0.0_kind_real
168 if(allocated(self%t )) self%t = 0.0_kind_real
169 if(allocated(self%ps )) self%ps = 0.0_kind_real
170 if(allocated(self%q )) self%q = 0.0_kind_real
171 if(allocated(self%qi )) self%qi = 0.0_kind_real
172 if(allocated(self%ql )) self%ql = 0.0_kind_real
173 if(allocated(self%o3 )) self%o3 = 0.0_kind_real
174 
175 if(allocated(self%psi)) self%psi = 0.0_kind_real
176 if(allocated(self%chi)) self%chi = 0.0_kind_real
177 if(allocated(self%tv )) self%tv = 0.0_kind_real
178 if(allocated(self%qc )) self%qc = 0.0_kind_real
179 if(allocated(self%qic)) self%qic = 0.0_kind_real
180 if(allocated(self%qlc)) self%qlc = 0.0_kind_real
181 if(allocated(self%o3c)) self%o3c = 0.0_kind_real
182 
183 if(allocated(self%w )) self%w = 0.0_kind_real
184 if(allocated(self%delz)) self%delz = 0.0_kind_real
185 if(allocated(self%delp)) self%delp = 0.0_kind_real
186 
187 end subroutine zeros
188 
189 ! ------------------------------------------------------------------------------
190 
191 subroutine ones(self)
192 implicit none
193 type(fv3jedi_increment), intent(inout) :: self
194 
195 call zeros(self)
196 
197 if(allocated(self%ua )) self%ua = 1.0_kind_real
198 if(allocated(self%va )) self%va = 1.0_kind_real
199 if(allocated(self%t )) self%t = 1.0_kind_real
200 if(allocated(self%ps )) self%ps = 1.0_kind_real
201 if(allocated(self%q )) self%q = 1.0_kind_real
202 if(allocated(self%qi )) self%qi = 1.0_kind_real
203 if(allocated(self%ql )) self%ql = 1.0_kind_real
204 if(allocated(self%o3 )) self%o3 = 1.0_kind_real
205 
206 if(allocated(self%psi)) self%psi = 1.0_kind_real
207 if(allocated(self%chi)) self%chi = 1.0_kind_real
208 if(allocated(self%tv )) self%tv = 1.0_kind_real
209 if(allocated(self%qc )) self%qc = 1.0_kind_real
210 if(allocated(self%qic)) self%qic = 1.0_kind_real
211 if(allocated(self%qlc)) self%qlc = 1.0_kind_real
212 if(allocated(self%o3c)) self%o3c = 1.0_kind_real
213 
214 if(allocated(self%w )) self%w = 1.0_kind_real
215 if(allocated(self%delz)) self%delz = 1.0_kind_real
216 if(allocated(self%delp)) self%delp = 1.0_kind_real
217 
218 end subroutine ones
219 
220 ! ------------------------------------------------------------------------------
221 
222 subroutine random(self)
223 use random_vectors_mod
224 implicit none
225 type(fv3jedi_increment), intent(inout) :: self
226 integer :: nq
227 
228 if(allocated(self%ua )) call random_vector(self%ua )
229 if(allocated(self%va )) call random_vector(self%va )
230 if(allocated(self%t )) call random_vector(self%t )
231 if(allocated(self%ps )) call random_vector(self%ps )
232 if(allocated(self%q )) call random_vector(self%q )
233 if(allocated(self%qi )) call random_vector(self%qi )
234 if(allocated(self%ql )) call random_vector(self%ql )
235 if(allocated(self%o3 )) call random_vector(self%o3 )
236 
237 if(allocated(self%psi )) call random_vector(self%psi )
238 if(allocated(self%chi )) call random_vector(self%chi )
239 if(allocated(self%tv )) call random_vector(self%tv )
240 if(allocated(self%qc )) call random_vector(self%qc )
241 if(allocated(self%qic )) call random_vector(self%qic )
242 if(allocated(self%qlc )) call random_vector(self%qlc )
243 if(allocated(self%o3c )) call random_vector(self%o3c )
244 
245 if(allocated(self%w )) call random_vector(self%w )
246 if(allocated(self%delz)) call random_vector(self%delz)
247 if(allocated(self%delp)) call random_vector(self%delp)
248 
249 end subroutine random
250 
251 ! ------------------------------------------------------------------------------
252 
253 subroutine copy(self,rhs)
254 implicit none
255 type(fv3jedi_increment), intent(inout) :: self
256 type(fv3jedi_increment), intent(in) :: rhs
257 
258 self%isc = rhs%isc
259 self%iec = rhs%iec
260 self%jsc = rhs%jsc
261 self%jec = rhs%jec
262 self%isd = rhs%isd
263 self%ied = rhs%ied
264 self%jsd = rhs%jsd
265 self%jed = rhs%jed
266 self%npx = rhs%npx
267 self%npy = rhs%npy
268 self%npz = rhs%npz
269 self%hydrostatic = rhs%hydrostatic
270 self%calendar_type = rhs%calendar_type
271 self%date = rhs%date
272 self%date_init = rhs%date_init
273 
274 if(allocated(self%ua )) self%ua = rhs%ua
275 if(allocated(self%va )) self%va = rhs%va
276 if(allocated(self%t )) self%t = rhs%t
277 if(allocated(self%ps )) self%ps = rhs%ps
278 if(allocated(self%q )) self%q = rhs%q
279 if(allocated(self%qi )) self%qi = rhs%qi
280 if(allocated(self%ql )) self%ql = rhs%ql
281 if(allocated(self%o3 )) self%o3 = rhs%o3
282 
283 if(allocated(self%psi )) self%psi = rhs%psi
284 if(allocated(self%chi )) self%chi = rhs%chi
285 if(allocated(self%tv )) self%tv = rhs%tv
286 if(allocated(self%qc )) self%qc = rhs%qc
287 if(allocated(self%qic )) self%qic = rhs%qic
288 if(allocated(self%qlc )) self%qlc = rhs%qlc
289 if(allocated(self%o3c )) self%o3c = rhs%o3c
290 
291 if(allocated(self%w )) self%w = rhs%w
292 if(allocated(self%delz)) self%delz = rhs%delz
293 if(allocated(self%delp)) self%delp = rhs%delp
294 
295 return
296 end subroutine copy
297 
298 ! ------------------------------------------------------------------------------
299 
300 subroutine self_add(self,rhs)
301 implicit none
302 type(fv3jedi_increment), intent(inout) :: self
303 type(fv3jedi_increment), intent(in) :: rhs
304 
305 if(allocated(self%ua )) self%ua = self%ua + rhs%ua
306 if(allocated(self%va )) self%va = self%va + rhs%va
307 if(allocated(self%t )) self%t = self%t + rhs%t
308 if(allocated(self%q )) self%q = self%q + rhs%q
309 if(allocated(self%ps )) self%ps = self%ps + rhs%ps
310 if(allocated(self%qi )) self%qi = self%qi + rhs%qi
311 if(allocated(self%ql )) self%ql = self%ql + rhs%ql
312 if(allocated(self%o3 )) self%o3 = self%o3 + rhs%o3
313 
314 if(allocated(self%psi )) self%psi = self%psi + rhs%psi
315 if(allocated(self%chi )) self%chi = self%chi + rhs%chi
316 if(allocated(self%tv )) self%tv = self%tv + rhs%tv
317 if(allocated(self%qc )) self%qc = self%qc + rhs%qc
318 if(allocated(self%qic )) self%qic = self%qic + rhs%qic
319 if(allocated(self%qlc )) self%qlc = self%qlc + rhs%qlc
320 if(allocated(self%o3c )) self%o3c = self%o3c + rhs%o3c
321 
322 if(allocated(self%w )) self%w = self%w + rhs%w
323 if(allocated(self%delz)) self%delz = self%delz + rhs%delz
324 if(allocated(self%delp)) self%delp = self%delp + rhs%delp
325 
326 return
327 end subroutine self_add
328 
329 ! ------------------------------------------------------------------------------
330 
331 subroutine self_schur(self,rhs)
332 implicit none
333 type(fv3jedi_increment), intent(inout) :: self
334 type(fv3jedi_increment), intent(in) :: rhs
335 
336 if(allocated(self%ua )) self%ua = self%ua * rhs%ua
337 if(allocated(self%va )) self%va = self%va * rhs%va
338 if(allocated(self%t )) self%t = self%t * rhs%t
339 if(allocated(self%ps )) self%ps = self%ps * rhs%ps
340 if(allocated(self%q )) self%q = self%q * rhs%q
341 if(allocated(self%qi )) self%qi = self%qi * rhs%qi
342 if(allocated(self%ql )) self%ql = self%ql * rhs%ql
343 if(allocated(self%o3 )) self%o3 = self%o3 * rhs%o3
344 
345 if(allocated(self%psi )) self%psi = self%psi * rhs%psi
346 if(allocated(self%chi )) self%chi = self%chi * rhs%chi
347 if(allocated(self%tv )) self%tv = self%tv * rhs%tv
348 if(allocated(self%qc )) self%qc = self%qc * rhs%qc
349 if(allocated(self%qic )) self%qic = self%qic * rhs%qic
350 if(allocated(self%qlc )) self%qlc = self%qlc * rhs%qlc
351 if(allocated(self%o3c )) self%o3c = self%o3c * rhs%o3c
352 
353 if(allocated(self%w )) self%w = self%w * rhs%w
354 if(allocated(self%delz)) self%delz = self%delz * rhs%delz
355 if(allocated(self%delp)) self%delp = self%delp * rhs%delp
356 
357 return
358 end subroutine self_schur
359 
360 ! ------------------------------------------------------------------------------
361 
362 subroutine self_sub(self,rhs)
363 implicit none
364 type(fv3jedi_increment), intent(inout) :: self
365 type(fv3jedi_increment), intent(in) :: rhs
366 
367 if(allocated(self%ua )) self%ua = self%ua - rhs%ua
368 if(allocated(self%va )) self%va = self%va - rhs%va
369 if(allocated(self%t )) self%t = self%t - rhs%t
370 if(allocated(self%ps )) self%ps = self%ps - rhs%ps
371 if(allocated(self%q )) self%q = self%q - rhs%q
372 if(allocated(self%qi )) self%qi = self%qi - rhs%qi
373 if(allocated(self%ql )) self%ql = self%ql - rhs%ql
374 if(allocated(self%o3 )) self%o3 = self%o3 - rhs%o3
375 
376 if(allocated(self%psi )) self%psi = self%psi - rhs%psi
377 if(allocated(self%chi )) self%chi = self%chi - rhs%chi
378 if(allocated(self%tv )) self%tv = self%tv - rhs%tv
379 if(allocated(self%qc )) self%qc = self%qc - rhs%qc
380 if(allocated(self%qic )) self%qic = self%qic - rhs%qic
381 if(allocated(self%qlc )) self%qlc = self%qlc - rhs%qlc
382 if(allocated(self%o3c )) self%o3c = self%o3c - rhs%o3c
383 
384 if(allocated(self%w )) self%w = self%w - rhs%w
385 if(allocated(self%delz)) self%delz = self%delz - rhs%delz
386 if(allocated(self%delp)) self%delp = self%delp - rhs%delp
387 
388 return
389 end subroutine self_sub
390 
391 ! ------------------------------------------------------------------------------
392 
393 subroutine self_mul(self,zz)
394 implicit none
395 type(fv3jedi_increment), intent(inout) :: self
396 real(kind=kind_real), intent(in) :: zz
397 
398 if(allocated(self%ua )) self%ua = zz * self%ua
399 if(allocated(self%va )) self%va = zz * self%va
400 if(allocated(self%t )) self%t = zz * self%t
401 if(allocated(self%ps )) self%ps = zz * self%ps
402 if(allocated(self%q )) self%q = zz * self%q
403 if(allocated(self%qi )) self%qi = zz * self%qi
404 if(allocated(self%ql )) self%ql = zz * self%ql
405 if(allocated(self%o3 )) self%o3 = zz * self%o3
406 
407 if(allocated(self%psi )) self%psi = zz * self%psi
408 if(allocated(self%chi )) self%chi = zz * self%chi
409 if(allocated(self%tv )) self%tv = zz * self%tv
410 if(allocated(self%qc )) self%qc = zz * self%qc
411 if(allocated(self%qic )) self%qic = zz * self%qic
412 if(allocated(self%qlc )) self%qlc = zz * self%qlc
413 if(allocated(self%o3c )) self%o3c = zz * self%o3c
414 
415 if(allocated(self%w )) self%w = zz * self%w
416 if(allocated(self%delz)) self%delz = zz * self%delz
417 if(allocated(self%delp)) self%delp = zz * self%delp
418 
419 return
420 end subroutine self_mul
421 
422 ! ------------------------------------------------------------------------------
423 
424 subroutine axpy_inc(self,zz,rhs)
425 implicit none
426 type(fv3jedi_increment), intent(inout) :: self
427 real(kind=kind_real), intent(in) :: zz
428 type(fv3jedi_increment), intent(in) :: rhs
429 
430 if(allocated(self%ua )) self%ua = self%ua + zz * rhs%ua
431 if(allocated(self%va )) self%va = self%va + zz * rhs%va
432 if(allocated(self%t )) self%t = self%t + zz * rhs%t
433 if(allocated(self%ps )) self%ps = self%ps + zz * rhs%ps
434 if(allocated(self%q )) self%q = self%q + zz * rhs%q
435 if(allocated(self%qi )) self%qi = self%qi + zz * rhs%qi
436 if(allocated(self%ql )) self%ql = self%ql + zz * rhs%ql
437 if(allocated(self%o3 )) self%o3 = self%o3 + zz * rhs%o3
438 
439 if(allocated(self%psi )) self%psi = self%psi + zz * rhs%psi
440 if(allocated(self%chi )) self%chi = self%chi + zz * rhs%chi
441 if(allocated(self%tv )) self%tv = self%tv + zz * rhs%tv
442 if(allocated(self%qc )) self%qc = self%qc + zz * rhs%qc
443 if(allocated(self%qic )) self%qic = self%qic + zz * rhs%qic
444 if(allocated(self%qlc )) self%qlc = self%qlc + zz * rhs%qlc
445 if(allocated(self%o3c )) self%o3c = self%o3c + zz * rhs%o3c
446 
447 if(allocated(self%w )) self%w = self%w + zz * rhs%w
448 if(allocated(self%delz)) self%delz = self%delz + zz * rhs%delz
449 if(allocated(self%delp)) self%delp = self%delp + zz * rhs%delp
450 
451 return
452 end subroutine axpy_inc
453 
454 ! ------------------------------------------------------------------------------
455 
456 subroutine axpy_state(self,zz,rhs)
457 implicit none
458 type(fv3jedi_increment), intent(inout) :: self
459 real(kind=kind_real), intent(in) :: zz
460 type(fv3jedi_state), intent(in) :: rhs
461 
462 real(kind=kind_real), allocatable :: rhs_ps(:,:)
463 
464 if(allocated(self%ua )) self%ua = self%ua + zz * rhs%ua
465 if(allocated(self%va )) self%va = self%va + zz * rhs%va
466 if(allocated(self%t )) self%t = self%t + zz * rhs%t
467 
468 if(allocated(self%ps))then
469  allocate(rhs_ps(rhs%isc:rhs%iec,rhs%jsc:rhs%jec))
470  rhs_ps = sum(rhs%delp,3)
471  self%ps = self%ps + zz * rhs_ps
472  deallocate(rhs_ps)
473 endif
474 
475 if(allocated(self%delp)) self%delp = self%delp + zz * rhs%delp
476 
477 if(allocated(self%q )) self%q = self%q + zz * rhs%q
478 if(allocated(self%qi )) self%qi = self%qi + zz * rhs%qi
479 if(allocated(self%ql )) self%ql = self%ql + zz * rhs%ql
480 if(allocated(self%o3 )) self%o3 = self%o3 + zz * rhs%o3
481 
482 if(allocated(self%w )) self%w = self%w + zz * rhs%w
483 if(allocated(self%delz)) self%delz = self%delz + zz * rhs%delz
484 
485 return
486 end subroutine axpy_state
487 
488 ! ------------------------------------------------------------------------------
489 
490 subroutine dot_prod(inc1,inc2,zprod)
491 
492 use fckit_mpi_module, only: fckit_mpi_comm, fckit_mpi_sum
493 
494 implicit none
495 type(fv3jedi_increment), intent(in) :: inc1, inc2
496 real(kind=kind_real), intent(inout) :: zprod
497 real(kind=kind_real) :: zp
498 integer :: i,j,k
499 integer :: ierr
500 type(fckit_mpi_comm) :: f_comm
501 
502 f_comm = fckit_mpi_comm()
503 
504 zp=0.0_kind_real
505 
506 !ua
507 if (allocated(inc1%ua)) then
508  do k = 1,inc1%npz
509  do j = inc1%jsc,inc1%jec
510  do i = inc1%isc,inc1%iec
511  zp = zp + inc1%ua(i,j,k) * inc2%ua(i,j,k)
512  enddo
513  enddo
514  enddo
515 endif
516 
517 !va
518 if (allocated(inc1%va)) then
519  do k = 1,inc1%npz
520  do j = inc1%jsc,inc1%jec
521  do i = inc1%isc,inc1%iec
522  zp = zp + inc1%va(i,j,k) * inc2%va(i,j,k)
523  enddo
524  enddo
525  enddo
526 endif
527 
528 !t
529 if (allocated(inc1%t)) then
530  do k = 1,inc1%npz
531  do j = inc1%jsc,inc1%jec
532  do i = inc1%isc,inc1%iec
533  zp = zp + inc1%t(i,j,k) * inc2%t(i,j,k)
534  enddo
535  enddo
536  enddo
537 endif
538 
539 !ps
540 if (allocated(inc1%ps)) then
541  do j = inc1%jsc,inc1%jec
542  do i = inc1%isc,inc1%iec
543  zp = zp + inc1%ps(i,j) * inc2%ps(i,j)
544  enddo
545  enddo
546 endif
547 
548 !q
549 if (allocated(inc1%q)) then
550  do k = 1,inc1%npz
551  do j = inc1%jsc,inc1%jec
552  do i = inc1%isc,inc1%iec
553  zp = zp + inc1%q(i,j,k) * inc2%q(i,j,k)
554  enddo
555  enddo
556  enddo
557 endif
558 
559 !qi
560 if (allocated(inc1%qi)) then
561  do k = 1,inc1%npz
562  do j = inc1%jsc,inc1%jec
563  do i = inc1%isc,inc1%iec
564  zp = zp + inc1%qi(i,j,k) * inc2%qi(i,j,k)
565  enddo
566  enddo
567  enddo
568 endif
569 
570 !ql
571 if (allocated(inc1%ql)) then
572  do k = 1,inc1%npz
573  do j = inc1%jsc,inc1%jec
574  do i = inc1%isc,inc1%iec
575  zp = zp + inc1%ql(i,j,k) * inc2%ql(i,j,k)
576  enddo
577  enddo
578  enddo
579 endif
580 
581 !o3
582 if (allocated(inc1%o3)) then
583  do k = 1,inc1%npz
584  do j = inc1%jsc,inc1%jec
585  do i = inc1%isc,inc1%iec
586  zp = zp + inc1%o3(i,j,k) * inc2%o3(i,j,k)
587  enddo
588  enddo
589  enddo
590 endif
591 
592 !psi
593 if (allocated(inc1%psi)) then
594  do k = 1,inc1%npz
595  do j = inc1%jsc,inc1%jec
596  do i = inc1%isc,inc1%iec
597  zp = zp + inc1%psi(i,j,k) * inc2%psi(i,j,k)
598  enddo
599  enddo
600  enddo
601 endif
602 
603 !chi
604 if (allocated(inc1%chi)) then
605  do k = 1,inc1%npz
606  do j = inc1%jsc,inc1%jec
607  do i = inc1%isc,inc1%iec
608  zp = zp + inc1%chi(i,j,k) * inc2%chi(i,j,k)
609  enddo
610  enddo
611  enddo
612 endif
613 
614 !tv
615 if (allocated(inc1%tv)) then
616  do k = 1,inc1%npz
617  do j = inc1%jsc,inc1%jec
618  do i = inc1%isc,inc1%iec
619  zp = zp + inc1%tv(i,j,k) * inc2%tv(i,j,k)
620  enddo
621  enddo
622  enddo
623 endif
624 
625 !qc
626 if (allocated(inc1%qc)) then
627  do k = 1,inc1%npz
628  do j = inc1%jsc,inc1%jec
629  do i = inc1%isc,inc1%iec
630  zp = zp + inc1%qc(i,j,k) * inc2%qc(i,j,k)
631  enddo
632  enddo
633  enddo
634 endif
635 
636 !qic
637 if (allocated(inc1%qic)) then
638  do k = 1,inc1%npz
639  do j = inc1%jsc,inc1%jec
640  do i = inc1%isc,inc1%iec
641  zp = zp + inc1%qic(i,j,k) * inc2%qic(i,j,k)
642  enddo
643  enddo
644  enddo
645 endif
646 
647 !qlc
648 if (allocated(inc1%qlc)) then
649  do k = 1,inc1%npz
650  do j = inc1%jsc,inc1%jec
651  do i = inc1%isc,inc1%iec
652  zp = zp + inc1%qlc(i,j,k) * inc2%qlc(i,j,k)
653  enddo
654  enddo
655  enddo
656 endif
657 
658 !o3c
659 if (allocated(inc1%o3c)) then
660  do k = 1,inc1%npz
661  do j = inc1%jsc,inc1%jec
662  do i = inc1%isc,inc1%iec
663  zp = zp + inc1%o3c(i,j,k) * inc2%o3c(i,j,k)
664  enddo
665  enddo
666  enddo
667 endif
668 
669 !delz
670 if (allocated(inc1%delz)) then
671  do k = 1,inc1%npz
672  do j = inc1%jsc,inc1%jec
673  do i = inc1%isc,inc1%iec
674  zp = zp + inc1%delz(i,j,k) * inc2%delz(i,j,k)
675  enddo
676  enddo
677  enddo
678 endif
679 
680 !delp
681 if (allocated(inc1%delp)) then
682  do k = 1,inc1%npz
683  do j = inc1%jsc,inc1%jec
684  do i = inc1%isc,inc1%iec
685  zp = zp + inc1%delp(i,j,k) * inc2%delp(i,j,k)
686  enddo
687  enddo
688  enddo
689 endif
690 
691 !w
692 if (allocated(inc1%w)) then
693  do k = 1,inc1%npz
694  do j = inc1%jsc,inc1%jec
695  do i = inc1%isc,inc1%iec
696  zp = zp + inc1%w(i,j,k) * inc2%w(i,j,k)
697  enddo
698  enddo
699  enddo
700 endif
701 
702 !Get global dot product
703 call f_comm%allreduce(zp,zprod,fckit_mpi_sum())
704 
705 !For debugging print result:
706 if (f_comm%rank() == 0) print*, "Dot product test result: ", zprod
707 
708 return
709 end subroutine dot_prod
710 
711 ! ------------------------------------------------------------------------------
712 
713 subroutine add_incr(self,rhs)
714 implicit none
715 type(fv3jedi_increment), intent(inout) :: self
716 type(fv3jedi_increment), intent(in) :: rhs
717 
718 integer :: check
719 check = (rhs%iec-rhs%isc+1) - (self%iec-self%isc+1)
720 
721 if (check==0) then
722  if(allocated(rhs%ua )) self%ua = self%ua + rhs%ua
723  if(allocated(rhs%va )) self%va = self%va + rhs%va
724  if(allocated(rhs%t )) self%t = self%t + rhs%t
725  if(allocated(rhs%ps )) self%ps = self%ps + rhs%ps
726  if(allocated(rhs%q )) self%q = self%q + rhs%q
727  if(allocated(rhs%qi )) self%qi = self%qi + rhs%qi
728  if(allocated(rhs%ql )) self%ql = self%ql + rhs%ql
729  if(allocated(rhs%o3 )) self%o3 = self%o3 + rhs%o3
730 
731  if(allocated(rhs%psi )) self%psi = self%psi + rhs%psi
732  if(allocated(rhs%chi )) self%chi = self%chi + rhs%chi
733  if(allocated(rhs%tv )) self%tv = self%tv + rhs%tv
734  if(allocated(rhs%qc )) self%qc = self%qc + rhs%qc
735  if(allocated(rhs%qic )) self%qic = self%qic + rhs%qic
736  if(allocated(rhs%qlc )) self%qlc = self%qlc + rhs%qlc
737  if(allocated(rhs%o3c )) self%o3c = self%o3c + rhs%o3c
738 
739  if(allocated(rhs%w )) self%w = self%w + rhs%w
740  if(allocated(rhs%delz)) self%delz = self%delz + rhs%delz
741  if(allocated(rhs%delp)) self%delp = self%delp + rhs%delp
742 else
743  call abor1_ftn("Increment: add_incr not implemented for low res increment yet")
744 endif
745 
746 return
747 end subroutine add_incr
748 
749 ! ------------------------------------------------------------------------------
750 
751 subroutine diff_incr(lhs,x1,x2)
752 implicit none
753 type(fv3jedi_increment), intent(inout) :: lhs
754 type(fv3jedi_state), intent(in) :: x1
755 type(fv3jedi_state), intent(in) :: x2
756 
757 real(kind=kind_real), allocatable :: x1_ps(:,:), x2_ps(:,:)
758 integer :: check
759 
760 check = (x1%iec-x1%isc+1) - (x2%iec-x2%isc+1)
761 
762 call zeros(lhs)
763 if (check==0) then
764 
765  if(allocated(lhs%ua )) lhs%ua = x1%ua - x2%ua
766  if(allocated(lhs%va )) lhs%va = x1%va - x2%va
767  if(allocated(lhs%t )) lhs%t = x1%t - x2%t
768 
769  if(allocated(lhs%ps)) then
770  allocate(x1_ps(x1%isc:x1%iec,x1%jsc:x1%jec))
771  allocate(x2_ps(x2%isc:x2%iec,x2%jsc:x2%jec))
772  x1_ps = sum(x1%delp,3)
773  x2_ps = sum(x2%delp,3)
774  lhs%ps = x1_ps - x2_ps
775  deallocate(x1_ps,x2_ps)
776  endif
777 
778  if(allocated(lhs%delp)) lhs%delp = x1%delp - x2%delp
779 
780  if(allocated(lhs%q )) lhs%q = x1%q - x2%q
781  if(allocated(lhs%qi )) lhs%qi = x1%qi - x2%qi
782  if(allocated(lhs%ql )) lhs%ql = x1%ql - x2%ql
783  if(allocated(lhs%o3 )) lhs%o3 = x1%o3 - x2%o3
784 
785  if(allocated(lhs%w )) lhs%w = x1%w - x2%w
786  if(allocated(lhs%delz)) lhs%delz = x1%delz - x2%delz
787 
788 else
789 
790  call abor1_ftn("Increment: diff_incr not implemented for low res increment yet")
791 
792 endif
793 
794 return
795 end subroutine diff_incr
796 
797 ! ------------------------------------------------------------------------------
798 
799 subroutine change_resol(inc,rhs)
800 implicit none
801 type(fv3jedi_increment), intent(inout) :: inc
802 type(fv3jedi_increment), intent(in) :: rhs
803 
804 integer :: check
805 check = (rhs%iec-rhs%isc+1) - (inc%iec-inc%isc+1)
806 
807 if (check==0) then
808  call copy(inc, rhs)
809 else
810  call abor1_ftn("Increment: change_resol not implmeneted yet")
811 endif
812 
813 return
814 end subroutine change_resol
815 
816 ! ------------------------------------------------------------------------------
817 
818 subroutine read_file(geom, inc, c_conf, vdate)
819 
820  implicit none
821 
822  type(fv3jedi_geom), intent(inout) :: geom !< Geometry
823  type(fv3jedi_increment), intent(inout) :: inc !< Increment
824  type(c_ptr), intent(in) :: c_conf !< Configuration
825  type(datetime), intent(inout) :: vdate !< DateTime
826 
827  character(len=10) :: restart_type
828 
829  restart_type = config_get_string(c_conf,len(restart_type),"restart_type")
830 
831  if (trim(restart_type) == 'gfs') then
832  call read_fms_restart(geom, inc, c_conf, vdate)
833  elseif (trim(restart_type) == 'geos') then
834  call read_geos_restart(geom, inc, c_conf, vdate)
835  else
836  call abor1_ftn("Increment: read restart type not supported")
837  endif
838 
839  return
840 
841 end subroutine read_file
842 
843 ! ------------------------------------------------------------------------------
844 
845 subroutine write_file(geom, inc, c_conf, vdate)
846 
847  implicit none
848 
849  type(fv3jedi_geom), intent(inout) :: geom !< Geometry
850  type(fv3jedi_increment), intent(in) :: inc !< Increment
851  type(c_ptr), intent(in) :: c_conf !< Configuration
852  type(datetime), intent(inout) :: vdate !< DateTime
853 
854  character(len=10) :: restart_type
855 
856  restart_type = config_get_string(c_conf,len(restart_type),"restart_type")
857 
858  if (trim(restart_type) == 'gfs') then
859  call write_fms_restart(geom, inc, c_conf, vdate)
860  elseif (trim(restart_type) == 'geos') then
861  call write_geos_restart(geom, inc, c_conf, vdate)
862  else
863  call abor1_ftn("Increment: write restart type not supported")
864  endif
865 
866  return
867 
868 end subroutine write_file
869 
870 ! ------------------------------------------------------------------------------
871 
872 subroutine gpnorm(inc, nf, pstat)
873 implicit none
874 type(fv3jedi_increment), intent(in) :: inc
875 integer, intent(in) :: nf
876 real(kind=kind_real), intent(inout) :: pstat(3, nf)
877 
878 integer :: isc, iec, jsc, jec, gs2, gs3
879 
880 !1. Min
881 !2. Max
882 !3. RMS
883 
884 isc = inc%isc
885 iec = inc%iec
886 jsc = inc%jsc
887 jec = inc%jec
888 
889 gs2 = (iec-isc+1)*(jec-jsc+1)
890 gs3 = gs2*inc%npz
891 
892 !ua
893 if (allocated(inc%ua)) then
894  pstat(1,1) = minval(inc%ua(isc:iec,jsc:jec,:))
895  pstat(2,1) = maxval(inc%ua(isc:iec,jsc:jec,:))
896  pstat(3,1) = sqrt((sum(inc%ua(isc:iec,jsc:jec,:))/gs3)**2)
897 endif
898 
899 !va
900 if (allocated(inc%va)) then
901  pstat(1,2) = minval(inc%va(isc:iec,jsc:jec,:))
902  pstat(2,2) = maxval(inc%va(isc:iec,jsc:jec,:))
903  pstat(3,2) = sqrt((sum(inc%va(isc:iec,jsc:jec,:))/gs3)**2)
904 endif
905 
906 !t
907 if (allocated(inc%t)) then
908  pstat(1,3) = minval(inc%t(isc:iec,jsc:jec,:))
909  pstat(2,3) = maxval(inc%t(isc:iec,jsc:jec,:))
910  pstat(3,3) = sqrt((sum(inc%t(isc:iec,jsc:jec,:))/gs3)**2)
911 endif
912 
913 !ps
914 if (allocated(inc%ps)) then
915  pstat(1,4) = minval(inc%ps(isc:iec,jsc:jec))
916  pstat(2,4) = maxval(inc%ps(isc:iec,jsc:jec))
917  pstat(3,4) = sqrt((sum(inc%ps(isc:iec,jsc:jec))/gs2)**2)
918 endif
919 
920 !q
921 if (allocated(inc%q)) then
922  pstat(1,5) = minval(inc%q(isc:iec,jsc:jec,:))
923  pstat(2,5) = maxval(inc%q(isc:iec,jsc:jec,:))
924  pstat(3,5) = sqrt((sum(inc%q(isc:iec,jsc:jec,:))/gs3)**2)
925 endif
926 
927 !qi
928 if (allocated(inc%qi)) then
929  pstat(1,6) = minval(inc%qi(isc:iec,jsc:jec,:))
930  pstat(2,6) = maxval(inc%qi(isc:iec,jsc:jec,:))
931  pstat(3,6) = sqrt((sum(inc%qi(isc:iec,jsc:jec,:))/gs3)**2)
932 endif
933 
934 !ql
935 if (allocated(inc%ql)) then
936  pstat(1,7) = minval(inc%ql(isc:iec,jsc:jec,:))
937  pstat(2,7) = maxval(inc%ql(isc:iec,jsc:jec,:))
938  pstat(3,7) = sqrt((sum(inc%ql(isc:iec,jsc:jec,:))/gs3)**2)
939 endif
940 
941 !o3
942 if (allocated(inc%o3)) then
943  pstat(1,8) = minval(inc%o3(isc:iec,jsc:jec,:))
944  pstat(2,8) = maxval(inc%o3(isc:iec,jsc:jec,:))
945  pstat(3,8) = sqrt((sum(inc%o3(isc:iec,jsc:jec,:))/gs3)**2)
946 endif
947 
948 !psi
949 if (allocated(inc%psi)) then
950  pstat(1,1) = minval(inc%psi(isc:iec,jsc:jec,:))
951  pstat(2,1) = maxval(inc%psi(isc:iec,jsc:jec,:))
952  pstat(3,1) = sqrt((sum(inc%psi(isc:iec,jsc:jec,:))/gs3)**2)
953 endif
954 
955 !chi
956 if (allocated(inc%chi)) then
957  pstat(1,2) = minval(inc%chi(isc:iec,jsc:jec,:))
958  pstat(2,2) = maxval(inc%chi(isc:iec,jsc:jec,:))
959  pstat(3,2) = sqrt((sum(inc%chi(isc:iec,jsc:jec,:))/gs3)**2)
960 endif
961 
962 !tv
963 if (allocated(inc%tv)) then
964  pstat(1,3) = minval(inc%tv(isc:iec,jsc:jec,:))
965  pstat(2,3) = maxval(inc%tv(isc:iec,jsc:jec,:))
966  pstat(3,3) = sqrt((sum(inc%tv(isc:iec,jsc:jec,:))/gs3)**2)
967 endif
968 
969 !qc
970 if (allocated(inc%qc)) then
971  pstat(1,5) = minval(inc%qc(isc:iec,jsc:jec,:))
972  pstat(2,5) = maxval(inc%qc(isc:iec,jsc:jec,:))
973  pstat(3,5) = sqrt((sum(inc%qc(isc:iec,jsc:jec,:))/gs3)**2)
974 endif
975 
976 !qic
977 if (allocated(inc%qic)) then
978  pstat(1,6) = minval(inc%qic(isc:iec,jsc:jec,:))
979  pstat(2,6) = maxval(inc%qic(isc:iec,jsc:jec,:))
980  pstat(3,6) = sqrt((sum(inc%qic(isc:iec,jsc:jec,:))/gs3)**2)
981 endif
982 
983 !qlc
984 if (allocated(inc%qlc)) then
985  pstat(1,7) = minval(inc%qlc(isc:iec,jsc:jec,:))
986  pstat(2,7) = maxval(inc%qlc(isc:iec,jsc:jec,:))
987  pstat(3,7) = sqrt((sum(inc%qlc(isc:iec,jsc:jec,:))/gs3)**2)
988 endif
989 
990 !o3c
991 if (allocated(inc%o3c)) then
992  pstat(1,8) = minval(inc%o3c(isc:iec,jsc:jec,:))
993  pstat(2,8) = maxval(inc%o3c(isc:iec,jsc:jec,:))
994  pstat(3,8) = sqrt((sum(inc%o3c(isc:iec,jsc:jec,:))/gs3)**2)
995 endif
996 
997 !w
998 if (allocated(inc%w)) then
999  pstat(1,9) = minval(inc%w(isc:iec,jsc:jec,:))
1000  pstat(2,9) = maxval(inc%w(isc:iec,jsc:jec,:))
1001  pstat(3,9) = sqrt((sum(inc%w(isc:iec,jsc:jec,:))/gs3)**2)
1002 endif
1003 
1004 !delz
1005 if (allocated(inc%delz)) then
1006  pstat(1,10) = minval(inc%delz(isc:iec,jsc:jec,:))
1007  pstat(2,10) = maxval(inc%delz(isc:iec,jsc:jec,:))
1008  pstat(3,10) = sqrt((sum(inc%delz(isc:iec,jsc:jec,:))/gs3)**2)
1009 endif
1010 
1011 !delp
1012 if (allocated(inc%delp)) then
1013  pstat(1,4) = minval(inc%delp(isc:iec,jsc:jec,:))
1014  pstat(2,4) = maxval(inc%delp(isc:iec,jsc:jec,:))
1015  pstat(3,4) = sqrt((sum(inc%delp(isc:iec,jsc:jec,:))/gs3)**2)
1016 endif
1017 
1018 return
1019 
1020 end subroutine gpnorm
1021 
1022 ! ------------------------------------------------------------------------------
1023 
1024 subroutine incrms(inc, prms)
1025 use fckit_mpi_module, only : fckit_mpi_comm, fckit_mpi_sum
1026 implicit none
1027 type(fv3jedi_increment), intent(in) :: inc
1028 real(kind=kind_real), intent(out) :: prms
1029 
1030 real(kind=kind_real) :: zz
1031 integer i,j,k,ii,nt,ierr,npes,iisum
1032 integer :: isc,iec,jsc,jec,npz
1033 type(fckit_mpi_comm) :: f_comm
1034 
1035 isc = inc%isc
1036 iec = inc%iec
1037 jsc = inc%jsc
1038 jec = inc%jec
1039 npz = inc%npz
1040 
1041 f_comm = fckit_mpi_comm()
1042 
1043 zz = 0.0_kind_real
1044 prms = 0.0_kind_real
1045 ii = 0
1046 
1047 !ua
1048 if (allocated(inc%ua)) then
1049  do k = 1,npz
1050  do j = jsc,jec
1051  do i = isc,iec
1052  zz = zz + inc%ua(i,j,k)**2
1053  ii = ii + 1
1054  enddo
1055  enddo
1056  enddo
1057 endif
1058 
1059 !va
1060 if (allocated(inc%va)) then
1061  do k = 1,npz
1062  do j = jsc,jec
1063  do i = isc,iec
1064  zz = zz + inc%va(i,j,k)**2
1065  ii = ii + 1
1066  enddo
1067  enddo
1068  enddo
1069 endif
1070 
1071 !t
1072 if (allocated(inc%t)) then
1073  do k = 1,npz
1074  do j = jsc,jec
1075  do i = isc,iec
1076  zz = zz + inc%t(i,j,k)**2
1077  ii = ii + 1
1078  enddo
1079  enddo
1080  enddo
1081 endif
1082 
1083 !ps
1084 if (allocated(inc%ps)) then
1085  do j = jsc,jec
1086  do i = isc,iec
1087  zz = zz + inc%ps(i,j)**2
1088  ii = ii + 1
1089  enddo
1090  enddo
1091 endif
1092 
1093 !q
1094 if (allocated(inc%q)) then
1095  do k = 1,npz
1096  do j = jsc,jec
1097  do i = isc,iec
1098  zz = zz + inc%q(i,j,k)**2
1099  ii = ii + 1
1100  enddo
1101  enddo
1102  enddo
1103 endif
1104 
1105 !qi
1106 if (allocated(inc%qi)) then
1107  do k = 1,npz
1108  do j = jsc,jec
1109  do i = isc,iec
1110  zz = zz + inc%qi(i,j,k)**2
1111  ii = ii + 1
1112  enddo
1113  enddo
1114  enddo
1115 endif
1116 
1117 !ql
1118 if (allocated(inc%ql)) then
1119  do k = 1,npz
1120  do j = jsc,jec
1121  do i = isc,iec
1122  zz = zz + inc%ql(i,j,k)**2
1123  ii = ii + 1
1124  enddo
1125  enddo
1126  enddo
1127 endif
1128 
1129 !o3
1130 if (allocated(inc%o3)) then
1131  do k = 1,npz
1132  do j = jsc,jec
1133  do i = isc,iec
1134  zz = zz + inc%o3(i,j,k)**2
1135  ii = ii + 1
1136  enddo
1137  enddo
1138  enddo
1139 endif
1140 
1141 !psi
1142 if (allocated(inc%psi)) then
1143  do k = 1,npz
1144  do j = jsc,jec
1145  do i = isc,iec
1146  zz = zz + inc%psi(i,j,k)**2
1147  ii = ii + 1
1148  enddo
1149  enddo
1150  enddo
1151 endif
1152 
1153 !chi
1154 if (allocated(inc%chi)) then
1155  do k = 1,npz
1156  do j = jsc,jec
1157  do i = isc,iec
1158  zz = zz + inc%chi(i,j,k)**2
1159  ii = ii + 1
1160  enddo
1161  enddo
1162  enddo
1163 endif
1164 
1165 !tv
1166 if (allocated(inc%tv)) then
1167  do k = 1,npz
1168  do j = jsc,jec
1169  do i = isc,iec
1170  zz = zz + inc%tv(i,j,k)**2
1171  ii = ii + 1
1172  enddo
1173  enddo
1174  enddo
1175 endif
1176 
1177 !qc
1178 if (allocated(inc%qc)) then
1179  do k = 1,npz
1180  do j = jsc,jec
1181  do i = isc,iec
1182  zz = zz + inc%qc(i,j,k)**2
1183  ii = ii + 1
1184  enddo
1185  enddo
1186  enddo
1187 endif
1188 
1189 !qic
1190 if (allocated(inc%qic)) then
1191  do k = 1,npz
1192  do j = jsc,jec
1193  do i = isc,iec
1194  zz = zz + inc%qic(i,j,k)**2
1195  ii = ii + 1
1196  enddo
1197  enddo
1198  enddo
1199 endif
1200 
1201 !qlc
1202 if (allocated(inc%qlc)) then
1203  do k = 1,npz
1204  do j = jsc,jec
1205  do i = isc,iec
1206  zz = zz + inc%qlc(i,j,k)**2
1207  ii = ii + 1
1208  enddo
1209  enddo
1210  enddo
1211 endif
1212 
1213 !o3c
1214 if (allocated(inc%o3c)) then
1215  do k = 1,npz
1216  do j = jsc,jec
1217  do i = isc,iec
1218  zz = zz + inc%o3c(i,j,k)**2
1219  ii = ii + 1
1220  enddo
1221  enddo
1222  enddo
1223 endif
1224 
1225 !w
1226 if (allocated(inc%w)) then
1227  do k = 1,npz
1228  do j = jsc,jec
1229  do i = isc,iec
1230  zz = zz + inc%w(i,j,k)**2
1231  ii = ii + 1
1232  enddo
1233  enddo
1234  enddo
1235 endif
1236 
1237 !delz
1238 if (allocated(inc%delz)) then
1239  do k = 1,npz
1240  do j = jsc,jec
1241  do i = isc,iec
1242  zz = zz + inc%delz(i,j,k)**2
1243  ii = ii + 1
1244  enddo
1245  enddo
1246  enddo
1247 endif
1248 
1249 !delp
1250 if (allocated(inc%delp)) then
1251  do k = 1,npz
1252  do j = jsc,jec
1253  do i = isc,iec
1254  zz = zz + inc%delp(i,j,k)**2
1255  ii = ii + 1
1256  enddo
1257  enddo
1258  enddo
1259 endif
1260 
1261 !Get global values
1262 call f_comm%allreduce(zz,prms,fckit_mpi_sum())
1263 call f_comm%allreduce(ii,iisum,fckit_mpi_sum())
1264 
1265 !if (ierr .ne. 0) then
1266 ! print *,'error in incrms/mpi_allreduce, error code=',ierr
1267 !endif
1268 prms = sqrt(prms/real(iisum,kind_real))
1269 
1270 end subroutine incrms
1271 
1272 ! ------------------------------------------------------------------------------
1273 
1274 subroutine dirac(self, c_conf, geom)
1275 
1276 use iso_c_binding
1277 
1278 implicit none
1279 type(fv3jedi_increment), intent(inout) :: self
1280 type(c_ptr), intent(in) :: c_conf !< Configuration
1281 type(fv3jedi_geom), intent(in) :: geom
1282 integer :: ndir,idir,ildir,ifdir,itiledir
1283 integer,allocatable :: ixdir(:),iydir(:)
1284 character(len=3) :: idirchar
1285 
1286 ! Get Diracs positions
1287 ndir = config_get_int(c_conf,"ndir")
1288 allocate(ixdir(ndir))
1289 allocate(iydir(ndir))
1290 
1291 do idir=1,ndir
1292  write(idirchar,'(i3)') idir
1293  ixdir(idir) = config_get_int(c_conf,"ixdir("//trim(adjustl(idirchar))//")")
1294  iydir(idir) = config_get_int(c_conf,"iydir("//trim(adjustl(idirchar))//")")
1295 end do
1296 ildir = config_get_int(c_conf,"ildir")
1297 ifdir = config_get_int(c_conf,"ifdir")
1298 itiledir = config_get_int(c_conf,"itiledir")
1299 
1300 
1301 ! Check
1302 if (ndir<1) call abor1_ftn("Increment: dirac non-positive ndir")
1303 if (any(ixdir<1).or.any(ixdir>self%npx)) then
1304  call abor1_ftn("Increment: dirac invalid ixdir")
1305 endif
1306 if (any(iydir<1).or.any(iydir>geom%size_cubic_grid)) then
1307  call abor1_ftn("Increment: dirac invalid iydir")
1308 endif
1309 if ((ildir<1).or.(ildir>self%npz)) then
1310  call abor1_ftn("Increment: dirac invalid ildir")
1311 endif
1312 if ((ifdir<1).or.(ifdir>5)) then
1313  call abor1_ftn("Increment: dirac invalid ifdir")
1314 endif
1315 if ((itiledir<1).or.(itiledir>6)) then
1316  call abor1_ftn("Increment: dirac invalid itiledir")
1317 endif
1318 
1319 ! Setup Diracs
1320 call zeros(self)
1321 
1322 ! only u, v, T, ps and tracers allowed
1323 do idir=1,ndir
1324 
1325  ! is specified grid point, tile number on this processor
1326  if (geom%ntile == itiledir .and. &
1327  ixdir(idir) >= self%isc .and. ixdir(idir) <= self%iec .and. &
1328  iydir(idir) >= self%jsc .and. iydir(idir) <= self%jec) then
1329  ! If so, perturb desired increment and level
1330  if (ifdir == 1) then
1331  self%ua (ixdir(idir),iydir(idir),ildir) = 1.0
1332  else if (ifdir == 2) then
1333  self%va (ixdir(idir),iydir(idir),ildir) = 1.0
1334  else if (ifdir == 3) then
1335  self%t (ixdir(idir),iydir(idir),ildir) = 1.0
1336  else if (ifdir == 4) then
1337  self%ps (ixdir(idir),iydir(idir) ) = 1.0
1338  else if (ifdir == 5) then
1339  self%q (ixdir(idir),iydir(idir),ildir) = 1.0
1340  else if (ifdir == 6) then
1341  self%qi (ixdir(idir),iydir(idir),ildir) = 1.0
1342  else if (ifdir == 7) then
1343  self%ql (ixdir(idir),iydir(idir),ildir) = 1.0
1344  else if (ifdir == 8) then
1345  self%o3 (ixdir(idir),iydir(idir),ildir) = 1.0
1346  endif
1347  endif
1348 end do
1349 
1350 end subroutine dirac
1351 
1352 ! ------------------------------------------------------------------------------
1353 
1354 subroutine ug_size(self, ug)
1355 use unstructured_grid_mod
1356 implicit none
1357 type(fv3jedi_increment), intent(in) :: self
1358 type(unstructured_grid), intent(inout) :: ug
1359 
1360 integer :: igrid
1361 
1362 ! Set number of grids
1363 if (ug%colocated==1) then
1364  ! Colocatd
1365  ug%ngrid = 1
1366 else
1367  ! Not colocated
1368  ug%ngrid = 1
1369  call abor1_ftn("Increment: Uncolocated grids not coded yet, and not needed")
1370 end if
1371 
1372 ! Allocate grid instances
1373 if (.not.allocated(ug%grid)) allocate(ug%grid(ug%ngrid))
1374 
1375 
1376 ! Set local number of points
1377 ug%grid(1)%nmga = (self%iec - self%isc + 1) * (self%jec - self%jsc + 1)
1378 
1379 ! Set number of levels
1380 ug%grid(1)%nl0 = self%npz
1381 
1382 ! Set number of variables
1383 ug%grid(1)%nv = 8
1384 
1385 if (.not. self%hydrostatic) then
1386  ug%grid(1)%nv = 10
1387 endif
1388 
1389 ! Set number of timeslots
1390 ug%grid(1)%nts = 1
1391 
1392 end subroutine ug_size
1393 
1394 ! ------------------------------------------------------------------------------
1395 
1396 subroutine ug_coord(self, ug, colocated, geom)
1397 use unstructured_grid_mod
1398 implicit none
1399 type(fv3jedi_increment), intent(in) :: self
1400 type(unstructured_grid), intent(inout) :: ug
1401 integer, intent(in) :: colocated
1402 type(fv3jedi_geom), intent(in) :: geom
1403 
1404 integer :: imga,jx,jy,jl
1405 real(kind=kind_real),allocatable :: lon(:),lat(:),area(:),vunit(:,:)
1406 real(kind=kind_real) :: sigmaup,sigmadn
1407 integer :: igrid
1408 
1409 ! Copy colocated
1410 ug%colocated = colocated
1411 
1412 ! Define size
1413 call ug_size(self, ug)
1414 
1415 ! Allocate unstructured grid coordinates
1416 call allocate_unstructured_grid_coord(ug)
1417 
1418 if (ug%colocated==1) then
1419  imga = 0
1420  do jy=self%jsc,self%jec
1421  do jx=self%isc,self%iec
1422  imga = imga+1
1423  ug%grid(1)%lon(imga) = rad2deg*geom%grid_lon(jx,jy)
1424  ug%grid(1)%lat(imga) = rad2deg*geom%grid_lat(jx,jy)
1425  ug%grid(1)%area(imga) = geom%area(jx,jy)
1426  do jl=1,self%npz
1427  sigmaup = geom%ak(jl+1)/101300.0+geom%bk(jl+1) ! si are now sigmas
1428  sigmadn = geom%ak(jl )/101300.0+geom%bk(jl )
1429  ug%grid(1)%vunit(imga,jl) = 0.5*(sigmaup+sigmadn) ! 'fake' sigma coordinates
1430  ug%grid(1)%lmask(imga,jl) = .true.
1431  enddo
1432  enddo
1433  enddo
1434 endif
1435 
1436 end subroutine ug_coord
1437 
1438 ! ------------------------------------------------------------------------------
1439 
1440 subroutine increment_to_ug(self, ug, colocated)
1441 use unstructured_grid_mod
1442 implicit none
1443 type(fv3jedi_increment), intent(in) :: self
1444 type(unstructured_grid), intent(inout) :: ug
1445 integer, intent(in) :: colocated
1446 
1447 integer :: imga,jx,jy,jl
1448 real(kind=kind_real),allocatable :: ptmp(:,:,:)
1449 integer :: igrid
1450 
1451 ! Copy colocated
1452 ug%colocated = colocated
1453 
1454 ! Define size
1455 call ug_size(self, ug)
1456 
1457 ! Allocate unstructured grid increment
1458 call allocate_unstructured_grid_field(ug)
1459 
1460 ! Copy increment
1461 
1462 ug%grid(1)%fld = 0.0_kind_real
1463 
1464 if (ug%colocated==1) then
1465 
1466  if (allocated(self%ua)) then
1467  imga = 0
1468  do jy=self%jsc,self%jec
1469  do jx=self%isc,self%iec
1470  imga = imga+1
1471  do jl=1,self%npz
1472  ug%grid(1)%fld(imga,jl,1,1) = self%ua (jx,jy,jl)
1473  enddo
1474  enddo
1475  enddo
1476  endif
1477 
1478  if (allocated(self%va)) then
1479  imga = 0
1480  do jy=self%jsc,self%jec
1481  do jx=self%isc,self%iec
1482  imga = imga+1
1483  do jl=1,self%npz
1484  ug%grid(1)%fld(imga,jl,2,1) = self%va (jx,jy,jl)
1485  enddo
1486  enddo
1487  enddo
1488  endif
1489 
1490  if (allocated(self%t)) then
1491  imga = 0
1492  do jy=self%jsc,self%jec
1493  do jx=self%isc,self%iec
1494  imga = imga+1
1495  do jl=1,self%npz
1496  ug%grid(1)%fld(imga,jl,3,1) = self%t (jx,jy,jl)
1497  enddo
1498  enddo
1499  enddo
1500  endif
1501 
1502  if (allocated(self%ps)) then
1503  imga = 0
1504  do jy=self%jsc,self%jec
1505  do jx=self%isc,self%iec
1506  imga = imga+1
1507  ug%grid(1)%fld(imga,self%npz,4,1) = self%ps (jx,jy)
1508  enddo
1509  enddo
1510  endif
1511 
1512  if (allocated(self%q)) then
1513  imga = 0
1514  do jy=self%jsc,self%jec
1515  do jx=self%isc,self%iec
1516  imga = imga+1
1517  do jl=1,self%npz
1518  ug%grid(1)%fld(imga,jl,5,1) = self%q (jx,jy,jl)
1519  enddo
1520  enddo
1521  enddo
1522  endif
1523 
1524  if (allocated(self%qi)) then
1525  imga = 0
1526  do jy=self%jsc,self%jec
1527  do jx=self%isc,self%iec
1528  imga = imga+1
1529  do jl=1,self%npz
1530  ug%grid(1)%fld(imga,jl,6,1) = self%qi (jx,jy,jl)
1531  enddo
1532  enddo
1533  enddo
1534  endif
1535 
1536  if (allocated(self%ql)) then
1537  imga = 0
1538  do jy=self%jsc,self%jec
1539  do jx=self%isc,self%iec
1540  imga = imga+1
1541  do jl=1,self%npz
1542  ug%grid(1)%fld(imga,jl,7,1) = self%ql (jx,jy,jl)
1543  enddo
1544  enddo
1545  enddo
1546  endif
1547 
1548  if (allocated(self%o3)) then
1549  imga = 0
1550  do jy=self%jsc,self%jec
1551  do jx=self%isc,self%iec
1552  imga = imga+1
1553  do jl=1,self%npz
1554  ug%grid(1)%fld(imga,jl,8,1) = self%o3 (jx,jy,jl)
1555  enddo
1556  enddo
1557  enddo
1558  endif
1559 
1560  if (allocated(self%psi)) then
1561  imga = 0
1562  do jy=self%jsc,self%jec
1563  do jx=self%isc,self%iec
1564  imga = imga+1
1565  do jl=1,self%npz
1566  ug%grid(1)%fld(imga,jl,1,1) = self%psi (jx,jy,jl)
1567  enddo
1568  enddo
1569  enddo
1570  endif
1571 
1572  if (allocated(self%chi)) then
1573  imga = 0
1574  do jy=self%jsc,self%jec
1575  do jx=self%isc,self%iec
1576  imga = imga+1
1577  do jl=1,self%npz
1578  ug%grid(1)%fld(imga,jl,2,1) = self%chi (jx,jy,jl)
1579  enddo
1580  enddo
1581  enddo
1582  endif
1583 
1584  if (allocated(self%tv)) then
1585  imga = 0
1586  do jy=self%jsc,self%jec
1587  do jx=self%isc,self%iec
1588  imga = imga+1
1589  do jl=1,self%npz
1590  ug%grid(1)%fld(imga,jl,3,1) = self%tv (jx,jy,jl)
1591  enddo
1592  enddo
1593  enddo
1594  endif
1595 
1596  if (allocated(self%qc)) then
1597  imga = 0
1598  do jy=self%jsc,self%jec
1599  do jx=self%isc,self%iec
1600  imga = imga+1
1601  do jl=1,self%npz
1602  ug%grid(1)%fld(imga,jl,5,1) = self%qc (jx,jy,jl)
1603  enddo
1604  enddo
1605  enddo
1606  endif
1607 
1608  if (allocated(self%qic)) then
1609  imga = 0
1610  do jy=self%jsc,self%jec
1611  do jx=self%isc,self%iec
1612  imga = imga+1
1613  do jl=1,self%npz
1614  ug%grid(1)%fld(imga,jl,6,1) = self%qic (jx,jy,jl)
1615  enddo
1616  enddo
1617  enddo
1618  endif
1619 
1620  if (allocated(self%qlc)) then
1621  imga = 0
1622  do jy=self%jsc,self%jec
1623  do jx=self%isc,self%iec
1624  imga = imga+1
1625  do jl=1,self%npz
1626  ug%grid(1)%fld(imga,jl,7,1) = self%qlc (jx,jy,jl)
1627  enddo
1628  enddo
1629  enddo
1630  endif
1631 
1632  if (allocated(self%o3c)) then
1633  imga = 0
1634  do jy=self%jsc,self%jec
1635  do jx=self%isc,self%iec
1636  imga = imga+1
1637  do jl=1,self%npz
1638  ug%grid(1)%fld(imga,jl,8,1) = self%o3c (jx,jy,jl)
1639  enddo
1640  enddo
1641  enddo
1642  endif
1643 
1644  if (allocated(self%w)) then
1645  imga = 0
1646  do jy=self%jsc,self%jec
1647  do jx=self%isc,self%iec
1648  imga = imga+1
1649  do jl=1,self%npz
1650  ug%grid(1)%fld(imga,jl,9,1) = self%w (jx,jy,jl)
1651  enddo
1652  enddo
1653  enddo
1654  endif
1655 
1656  if (allocated(self%delz)) then
1657  imga = 0
1658  do jy=self%jsc,self%jec
1659  do jx=self%isc,self%iec
1660  imga = imga+1
1661  do jl=1,self%npz
1662  ug%grid(1)%fld(imga,jl,10,1) = self%delz(jx,jy,jl)
1663  enddo
1664  enddo
1665  enddo
1666  endif
1667 
1668  if (allocated(self%delp)) then
1669  imga = 0
1670  do jy=self%jsc,self%jec
1671  do jx=self%isc,self%iec
1672  imga = imga+1
1673  do jl=1,self%npz
1674  ug%grid(1)%fld(imga,jl,4,1) = self%delp(jx,jy,jl)
1675  enddo
1676  enddo
1677  enddo
1678  endif
1679 
1680 endif
1681 
1682 end subroutine increment_to_ug
1683 
1684 ! -----------------------------------------------------------------------------
1685 
1686 subroutine increment_from_ug(self, ug)
1687 use unstructured_grid_mod
1688 implicit none
1689 type(fv3jedi_increment), intent(inout) :: self
1690 type(unstructured_grid), intent(in) :: ug
1691 
1692 integer :: imga,jx,jy,jl
1693 real(kind=kind_real),allocatable :: ptmp(:,:,:)
1694 integer :: igrid
1695 
1696 ! Copy increment
1697 
1698 if (allocated(self%ua)) then
1699  imga = 0
1700  do jy=self%jsc,self%jec
1701  do jx=self%isc,self%iec
1702  imga = imga+1
1703  do jl=1,self%npz
1704  self%ua (jx,jy,jl) = ug%grid(1)%fld(imga,jl,1,1)
1705  enddo
1706  enddo
1707  enddo
1708 endif
1709 
1710 if (allocated(self%va)) then
1711  imga = 0
1712  do jy=self%jsc,self%jec
1713  do jx=self%isc,self%iec
1714  imga = imga+1
1715  do jl=1,self%npz
1716  self%va (jx,jy,jl) = ug%grid(1)%fld(imga,jl,2,1)
1717  enddo
1718  enddo
1719  enddo
1720 endif
1721 
1722 if (allocated(self%t)) then
1723  imga = 0
1724  do jy=self%jsc,self%jec
1725  do jx=self%isc,self%iec
1726  imga = imga+1
1727  do jl=1,self%npz
1728  self%t (jx,jy,jl) = ug%grid(1)%fld(imga,jl,3,1)
1729  enddo
1730  enddo
1731  enddo
1732 endif
1733 
1734 if (allocated(self%ps)) then
1735  imga = 0
1736  do jy=self%jsc,self%jec
1737  do jx=self%isc,self%iec
1738  imga = imga+1
1739  self%ps (jx,jy) = ug%grid(1)%fld(imga,self%npz,4,1)
1740  enddo
1741  enddo
1742 endif
1743 
1744 if (allocated(self%q)) then
1745  imga = 0
1746  do jy=self%jsc,self%jec
1747  do jx=self%isc,self%iec
1748  imga = imga+1
1749  do jl=1,self%npz
1750  self%q (jx,jy,jl) = ug%grid(1)%fld(imga,jl,5,1)
1751  enddo
1752  enddo
1753  enddo
1754 endif
1755 
1756 if (allocated(self%qi)) then
1757  imga = 0
1758  do jy=self%jsc,self%jec
1759  do jx=self%isc,self%iec
1760  imga = imga+1
1761  do jl=1,self%npz
1762  self%qi (jx,jy,jl) = ug%grid(1)%fld(imga,jl,6,1)
1763  enddo
1764  enddo
1765  enddo
1766 endif
1767 
1768 if (allocated(self%ql)) then
1769  imga = 0
1770  do jy=self%jsc,self%jec
1771  do jx=self%isc,self%iec
1772  imga = imga+1
1773  do jl=1,self%npz
1774  self%ql (jx,jy,jl) = ug%grid(1)%fld(imga,jl,7,1)
1775  enddo
1776  enddo
1777  enddo
1778 endif
1779 
1780 if (allocated(self%o3)) then
1781  imga = 0
1782  do jy=self%jsc,self%jec
1783  do jx=self%isc,self%iec
1784  imga = imga+1
1785  do jl=1,self%npz
1786  self%o3 (jx,jy,jl) = ug%grid(1)%fld(imga,jl,8,1)
1787  enddo
1788  enddo
1789  enddo
1790 endif
1791 
1792 if (allocated(self%psi)) then
1793  imga = 0
1794  do jy=self%jsc,self%jec
1795  do jx=self%isc,self%iec
1796  imga = imga+1
1797  do jl=1,self%npz
1798  self%psi (jx,jy,jl) = ug%grid(1)%fld(imga,jl,1,1)
1799  enddo
1800  enddo
1801  enddo
1802 endif
1803 
1804 if (allocated(self%chi)) then
1805  imga = 0
1806  do jy=self%jsc,self%jec
1807  do jx=self%isc,self%iec
1808  imga = imga+1
1809  do jl=1,self%npz
1810  self%chi (jx,jy,jl) = ug%grid(1)%fld(imga,jl,2,1)
1811  enddo
1812  enddo
1813  enddo
1814 endif
1815 
1816 if (allocated(self%tv)) then
1817  imga = 0
1818  do jy=self%jsc,self%jec
1819  do jx=self%isc,self%iec
1820  imga = imga+1
1821  do jl=1,self%npz
1822  self%tv (jx,jy,jl) = ug%grid(1)%fld(imga,jl,3,1)
1823  enddo
1824  enddo
1825  enddo
1826 endif
1827 
1828 if (allocated(self%qc)) then
1829  imga = 0
1830  do jy=self%jsc,self%jec
1831  do jx=self%isc,self%iec
1832  imga = imga+1
1833  do jl=1,self%npz
1834  self%qc (jx,jy,jl) = ug%grid(1)%fld(imga,jl,5,1)
1835  enddo
1836  enddo
1837  enddo
1838 endif
1839 
1840 if (allocated(self%qic)) then
1841  imga = 0
1842  do jy=self%jsc,self%jec
1843  do jx=self%isc,self%iec
1844  imga = imga+1
1845  do jl=1,self%npz
1846  self%qic (jx,jy,jl) = ug%grid(1)%fld(imga,jl,6,1)
1847  enddo
1848  enddo
1849  enddo
1850 endif
1851 
1852 if (allocated(self%qlc)) then
1853  imga = 0
1854  do jy=self%jsc,self%jec
1855  do jx=self%isc,self%iec
1856  imga = imga+1
1857  do jl=1,self%npz
1858  self%qlc (jx,jy,jl) = ug%grid(1)%fld(imga,jl,7,1)
1859  enddo
1860  enddo
1861  enddo
1862 endif
1863 
1864 if (allocated(self%o3c)) then
1865  imga = 0
1866  do jy=self%jsc,self%jec
1867  do jx=self%isc,self%iec
1868  imga = imga+1
1869  do jl=1,self%npz
1870  self%o3c (jx,jy,jl) = ug%grid(1)%fld(imga,jl,8,1)
1871  enddo
1872  enddo
1873  enddo
1874 endif
1875 
1876 if (allocated(self%w)) then
1877  imga = 0
1878  do jy=self%jsc,self%jec
1879  do jx=self%isc,self%iec
1880  imga = imga+1
1881  do jl=1,self%npz
1882  self%w (jx,jy,jl) = ug%grid(1)%fld(imga,jl,9,1)
1883  enddo
1884  enddo
1885  enddo
1886 endif
1887 
1888 if (allocated(self%delz)) then
1889  imga = 0
1890  do jy=self%jsc,self%jec
1891  do jx=self%isc,self%iec
1892  imga = imga+1
1893  do jl=1,self%npz
1894  self%delz(jx,jy,jl) = ug%grid(1)%fld(imga,jl,10,1)
1895  enddo
1896  enddo
1897  enddo
1898 endif
1899 
1900 if (allocated(self%delp)) then
1901  imga = 0
1902  do jy=self%jsc,self%jec
1903  do jx=self%isc,self%iec
1904  imga = imga+1
1905  do jl=1,self%npz
1906  self%delp(jx,jy,jl) = ug%grid(1)%fld(imga,jl,4,1)
1907  enddo
1908  enddo
1909  enddo
1910 endif
1911 
1912 end subroutine increment_from_ug
1913 
1914 ! ------------------------------------------------------------------------------
1915 
1916 subroutine check_resolution(x1, x2)
1917 
1918 implicit none
1919 type(fv3jedi_increment), intent(in) :: x1, x2
1920 
1921 if (x1%npx /= x2%npx .or. x1%npz /= x2%npz) then
1922  call abor1_ftn ("Increment: resolution error")
1923 endif
1924 call check(x1)
1925 call check(x2)
1926 
1927 end subroutine check_resolution
1928 
1929 ! ------------------------------------------------------------------------------
1930 
1931 subroutine check(self)
1932 implicit none
1933 type(fv3jedi_increment), intent(in) :: self
1934 logical :: bad
1935 
1936 bad = .false.
1937 
1938 if (bad) then
1939  call abor1_ftn ("Increment: increment not consistent")
1940 endif
1941 
1942 end subroutine check
1943 
1944 ! ------------------------------------------------------------------------------
1945 
1946 subroutine jnormgrad(self,geom,ref,c_conf)
1947 
1948 use mpp_domains_mod, only: mpp_global_sum, bitwise_efp_sum
1949 
1950 implicit none
1951 type(fv3jedi_increment) :: self
1952 type(fv3jedi_geom) :: geom
1953 type(fv3jedi_state) :: ref !To linearize around if nl
1954 type(c_ptr) :: c_conf
1955 
1956 integer :: i,j,k
1957 integer :: isc,iec,jsc,jec,npz
1958 real(kind=kind_real), allocatable :: cellweight(:,:,:), ref_ps(:,:)
1959 
1960 real(kind=kind_real) :: global_area
1961 
1962 real(kind=kind_real) :: ufac
1963 real(kind=kind_real) :: tfac, tref
1964 real(kind=kind_real) :: qfac, qeps
1965 real(kind=kind_real) :: pfac, pref
1966 
1967 
1968 !Code to compute a vector norm for an increment, e.g. the energy norm for FSOI
1969 
1970 isc = self%isc
1971 iec = self%iec
1972 jsc = self%jsc
1973 jec = self%jec
1974 npz = self%npz
1975 
1976 ! Constants
1977 ! ---------
1978 tref = config_get_real(c_conf,"Tref")
1979 qeps = config_get_real(c_conf,"qepsilon")
1980 pref = config_get_real(c_conf,"pref")
1981 
1982 ufac = 0.5_kind_real
1983 tfac = 0.5_kind_real*cp/tref
1984 qfac = 0.5_kind_real*qeps*alhl*alhl/(cp*tref)
1985 pfac = 0.5_kind_real*rgas*tref/pref**2
1986 
1987 ! Compute grid weighting based on volume
1988 ! --------------------------------------
1989 
1990 global_area = mpp_global_sum(geom%domain, geom%area, flags=bitwise_efp_sum)
1991 
1992 allocate(ref_ps(isc:iec,jsc:jec))
1993 ref_ps = sum(ref%delp,3)
1994 
1995 allocate(cellweight(isc:iec,jsc:jec,1:npz))
1996 do k = 1, npz
1997  do j = jsc,jec
1998  do i = isc,iec
1999  cellweight(i,j,k) = (ref%delp(i,j,k)/ref_ps(i,j)) * geom%area(i,j)/global_area
2000  enddo
2001  enddo
2002 enddo
2003 
2004 !ua
2005 do k = 1, npz
2006  do j = jsc,jec
2007  do i = isc,iec
2008  self%ua(i,j,k) = ufac * 2.0_kind_real * ref%ua(i,j,k) * cellweight(i,j,k)
2009  enddo
2010  enddo
2011 enddo
2012 
2013 !va
2014 do k = 1, npz
2015  do j = jsc,jec
2016  do i = isc,iec
2017  self%va(i,j,k) = ufac * 2.0_kind_real * ref%va(i,j,k) * cellweight(i,j,k)
2018  enddo
2019  enddo
2020 enddo
2021 
2022 !T
2023 do k = 1, npz
2024  do j = jsc,jec
2025  do i = isc,iec
2026  self%t(i,j,k) = tfac * 2.0_kind_real * ref%T(i,j,k) * cellweight(i,j,k)
2027  enddo
2028  enddo
2029 enddo
2030 
2031 !q
2032 do k = 1, npz
2033  do j = jsc,jec
2034  do i = isc,iec
2035  self%q(i,j,k) = qfac * 2.0_kind_real * ref%q (i,j,k) * cellweight(i,j,k)
2036  enddo
2037  enddo
2038 enddo
2039 
2040 !ps
2041 if (allocated(self%ps)) then
2042  do j = jsc,jec
2043  do i = isc,iec
2044  self%ps(i,j) = pfac * 2.0_kind_real * ref_ps(i,j) * cellweight(i,j,npz) / (ref%delp(i,j,npz)/ref_ps(i,j))
2045  enddo
2046  enddo
2047 else
2048  call abor1_ftn("Increment: JGradNorm does not support not using Ps in the increment yet")
2049 endif
2050 
2051 deallocate(cellweight)
2052 deallocate(ref_ps)
2053 
2054 end subroutine jnormgrad
2055 
2056 ! ------------------------------------------------------------------------------
2057 
2058 end module fv3jedi_increment_mod
fv3jedi_increment_mod::gpnorm
subroutine, public gpnorm(inc, nf, pstat)
Definition: fv3jedi_increment_mod.f90:873
unstructured_grid_mod::allocate_unstructured_grid_coord
subroutine, public allocate_unstructured_grid_coord(self)
Definition: unstructured_grid_mod.F90:90
random_vectors_mod::random_vector
Fortran generic for generating random 1d, 2d and 3d arrays.
Definition: fv3-jedi/src/Utilities/random_vectors_mod.f90:20
fv3jedi_increment_utils_mod::fv3jedi_increment
Fortran derived type to hold FV3JEDI increment.
Definition: fv3jedi_increment_utils_mod.F90:18
fv3jedi_increment_mod::dirac
subroutine, public dirac(self, c_conf, geom)
Definition: fv3jedi_increment_mod.f90:1275
fv3jedi_increment_mod::change_resol
subroutine, public change_resol(inc, rhs)
Definition: fv3jedi_increment_mod.f90:800
fv3jedi_increment_mod::self_schur
subroutine, public self_schur(self, rhs)
Definition: fv3jedi_increment_mod.f90:332
fv3jedi_constants_mod::rad2deg
real(kind=kind_real), parameter, public rad2deg
Definition: fv3jedi_constants_mod.f90:13
fv3jedi_increment_mod::axpy_state
subroutine, public axpy_state(self, zz, rhs)
Definition: fv3jedi_increment_mod.f90:457
fv3jedi_increment_mod::write_file
subroutine, public write_file(geom, inc, c_conf, vdate)
Definition: fv3jedi_increment_mod.f90:846
fv3jedi_increment_mod::dot_prod
subroutine, public dot_prod(inc1, inc2, zprod)
Definition: fv3jedi_increment_mod.f90:491
fv3jedi_increment_mod::self_add
subroutine, public self_add(self, rhs)
Definition: fv3jedi_increment_mod.f90:301
fv3jedi_increment_mod::copy
subroutine, public copy(self, rhs)
Definition: fv3jedi_increment_mod.f90:254
fv3jedi_constants_mod::rgas
real(kind=kind_real), parameter, public rgas
Definition: fv3jedi_constants_mod.f90:39
fv3jedi_state_utils_mod::fv3jedi_state
Fortran derived type to hold FV3JEDI state.
Definition: fv3jedi_state_utils_mod.F90:18
fv3jedi_increment_mod::jnormgrad
subroutine, public jnormgrad(self, geom, ref, c_conf)
Definition: fv3jedi_increment_mod.f90:1947
fv3jedi_increment_mod::incrms
subroutine, public incrms(inc, prms)
Definition: fv3jedi_increment_mod.f90:1025
fv3jedi_increment_mod::diff_incr
subroutine, public diff_incr(lhs, x1, x2)
Definition: fv3jedi_increment_mod.f90:752
fv3jedi_increment_mod::self_sub
subroutine, public self_sub(self, rhs)
Definition: fv3jedi_increment_mod.f90:363
fv3jedi_geom_mod::fv3jedi_geom
Fortran derived type to hold geometry data for the FV3JEDI model.
Definition: fv3jedi_geom_mod.f90:34
fv3jedi_getvalues_mod::getvalues_ad
subroutine, public getvalues_ad(geom, inc, locs, vars, gom, traj)
Definition: fv3jedi_getvalues_mod.f90:763
fv3jedi_increment_mod::ug_coord
subroutine, public ug_coord(self, ug, colocated, geom)
Definition: fv3jedi_increment_mod.f90:1397
fv3jedi_increment_io_mod::write_fms_restart
subroutine, public write_fms_restart(geom, incr, c_conf, vdate)
Definition: fv3jedi_increment_io_mod.f90:48
mpp_domains_mod
Definition: mpp_domains.F90:128
fv3jedi_increment_mod::random
subroutine, public random(self)
Definition: fv3jedi_increment_mod.f90:223
mpp_domains_mod::mpp_global_sum
Definition: mpp_domains.F90:2390
fv3jedi_vars_mod
Fortran module to handle variables for the FV3JEDI model.
Definition: fv3jedi_vars_mod.f90:10
fv3jedi_increment_mod::delete
subroutine, public delete(self)
Definition: fv3jedi_increment_mod.f90:132
fv3jedi_constants_mod
Definition: fv3jedi_constants_mod.f90:6
fv3jedi_increment_mod::increment_to_ug
subroutine, public increment_to_ug(self, ug, colocated)
Definition: fv3jedi_increment_mod.f90:1441
fv3jedi_getvalues_mod::getvalues_tl
subroutine, public getvalues_tl(geom, inc, locs, vars, gom, traj)
Definition: fv3jedi_getvalues_mod.f90:601
unstructured_grid_mod::unstructured_grid
Definition: unstructured_grid_mod.F90:37
fv3jedi_increment_mod::add_incr
subroutine, public add_incr(self, rhs)
Definition: fv3jedi_increment_mod.f90:714
unstructured_grid_mod::allocate_unstructured_grid_field
subroutine, public allocate_unstructured_grid_field(self)
Definition: unstructured_grid_mod.F90:108
fv3jedi_increment_mod::check
subroutine check(self)
Definition: fv3jedi_increment_mod.f90:1932
fv3jedi_increment_utils_mod
Utilities for increment for the FV3JEDI model.
Definition: fv3jedi_increment_utils_mod.F90:8
fv3jedi_increment_io_mod::write_geos_restart
subroutine, public write_geos_restart(geom, incr, c_conf, vdate)
Definition: fv3jedi_increment_io_mod.f90:211
fv3jedi_increment_mod
Handle increment for the FV3JEDI model.
Definition: fv3jedi_increment_mod.f90:8
fv3jedi_increment_mod::axpy_inc
subroutine, public axpy_inc(self, zz, rhs)
Definition: fv3jedi_increment_mod.f90:425
fv3jedi_getvalues_mod
Definition: fv3jedi_getvalues_mod.f90:1
random_vectors_mod
Fortran module for generating random vectors.
Definition: fv3-jedi/src/Utilities/random_vectors_mod.f90:9
fv3jedi_state_utils_mod
Utilities for state for the FV3JEDI model.
Definition: fv3jedi_state_utils_mod.F90:8
fv3jedi_increment_mod::read_file
subroutine, public read_file(geom, inc, c_conf, vdate)
Definition: fv3jedi_increment_mod.f90:819
fv3jedi_geom_mod
Fortran module handling geometry for the FV3 model.
Definition: fv3jedi_geom_mod.f90:8
fv3jedi_kinds_mod
Definition: fv3jedi_kinds_mod.f90:6
fv3jedi_increment_mod::ug_size
subroutine ug_size(self, ug)
Definition: fv3jedi_increment_mod.f90:1355
fv3jedi_increment_mod::zeros
subroutine, public zeros(self)
Definition: fv3jedi_increment_mod.f90:161
fv3jedi_increment_io_mod
Definition: fv3jedi_increment_io_mod.f90:1
fv3jedi_constants_mod::constoz
real(kind=kind_real), parameter, public constoz
Definition: fv3jedi_constants_mod.f90:57
fv3jedi_increment_mod::increment_from_ug
subroutine, public increment_from_ug(self, ug)
Definition: fv3jedi_increment_mod.f90:1687
fv3jedi_vars_mod::fv3jedi_vars
Fortran derived type to represent fv3jedi model variables.
Definition: fv3jedi_vars_mod.f90:22
fv3jedi_increment_mod::check_resolution
subroutine check_resolution(x1, x2)
Definition: fv3jedi_increment_mod.f90:1917
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:13
unstructured_grid_mod
Fortran module for handling generic unstructured grid.
Definition: unstructured_grid_mod.F90:8
fv3jedi_increment_mod::ones
subroutine ones(self)
Definition: fv3jedi_increment_mod.f90:192
fv3jedi_increment_mod::create
subroutine, public create(self, geom, vars)
Definition: fv3jedi_increment_mod.f90:42
fv3jedi_increment_mod::self_mul
subroutine, public self_mul(self, zz)
Definition: fv3jedi_increment_mod.f90:394
fv3jedi_increment_io_mod::read_geos_restart
subroutine, public read_geos_restart(geom, incr, c_conf, vdate)
Definition: fv3jedi_increment_io_mod.f90:196
fv3jedi_constants_mod::alhl
real(kind=kind_real), parameter, public alhl
Definition: fv3jedi_constants_mod.f90:25
fv3jedi_constants_mod::cp
real(kind=kind_real), parameter, public cp
Definition: fv3jedi_constants_mod.f90:40
fv3jedi_increment_io_mod::read_fms_restart
subroutine, public read_fms_restart(geom, incr, c_conf, vdate)
Definition: fv3jedi_increment_io_mod.f90:33