FV3DYCORE  Version 2.0.0
fv_nggps_diag.F90
Go to the documentation of this file.
1 
2 !***********************************************************************
3 !* GNU Lesser General Public License
4 !*
5 !* This file is part of the FV3 dynamical core.
6 !*
7 !* The FV3 dynamical core is free software: you can redistribute it
8 !* and/or modify it under the terms of the
9 !* GNU Lesser General Public License as published by the
10 !* Free Software Foundation, either version 3 of the License, or
11 !* (at your option) any later version.
12 !*
13 !* The FV3 dynamical core is distributed in the hope that it will be
14 !* useful, but WITHOUT ANYWARRANTY; without even the implied warranty
15 !* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16 !* See the GNU General Public License for more details.
17 !*
18 !* You should have received a copy of the GNU Lesser General Public
19 !* License along with the FV3 dynamical core.
20 !* If not, see <http://www.gnu.org/licenses/>.
21 !***********************************************************************
22 
27 
29 
30 ! <table>
31 ! <tr>
32 ! <th>Module Name</th>
33 ! <th>Functions Included</th>
34 ! </tr>
35 ! <tr>
36 ! <td>constants_mod</td>
37 ! <td>kappa, grav, rdgas</td>
38 ! </tr>
39 ! <tr>
40 ! <td>diag_manager_mod</td>
41 ! <td>register_diag_field, send_data</td>
42 ! </tr>
43 ! <tr>
44 ! <td>field_manager_mod</td>
45 ! <td>MODEL_ATMOS</td>
46 ! </tr>
47 ! <tr>
48 ! <td>fms_io_mod</td>
49 ! <td>set_domain, nullify_domain</td>
50 ! </tr>
51 ! <tr>
52 ! <td>fv_arrays_mod</td>
53 ! <td>fv_atmos_type</td>
54 ! </tr>
55 ! <tr>
56 ! <td>fv_diagnostics_mod</td>
57 ! <td>range_check</td>
58 ! </tr>
59 ! <tr>
60 ! <td>mpp_mod</td>
61 ! <td>mpp_pe, mpp_root_pe</td>
62 ! </tr>
63 ! <tr>
64 ! <td>tracer_manager_mod</td>
65 ! <td>get_tracer_names, get_number_tracers, get_tracer_index</td>
66 ! </tr>
67 ! </table>
68 
69  use mpp_mod, only: mpp_pe, mpp_root_pe,fatal,mpp_error
70  use constants_mod, only: grav, rdgas
71  use fms_io_mod, only: set_domain, nullify_domain
72  use time_manager_mod, only: time_type, get_time
73  use diag_manager_mod, only: register_diag_field, send_data
74  use diag_axis_mod, only: get_axis_global_length, get_diag_axis, get_diag_axis_name
75  use diag_data_mod, only: output_fields, max_output_fields
76  use diag_util_mod, only: find_input_field
77  use tracer_manager_mod, only: get_tracer_names, get_number_tracers, get_tracer_index
78  use field_manager_mod, only: model_atmos
82  use fv_arrays_mod, only: fv_atmos_type
83  use mpp_domains_mod, only: domain1d, domainug
84 #ifdef MULTI_GASES
85  use multi_gases_mod, only: virq
86 #endif
87 
88  implicit none
89  private
90 
91  real, parameter:: missing_value = -1.e10
92  real, parameter:: stndrd_atmos_ps = 101325.
93  real, parameter:: stndrd_atmos_lapse = 0.0065
94 
95  logical master
113  integer :: isco, ieco, jsco, jeco, npzo, ncnsto
114  integer :: isdo, iedo, jsdo, jedo
115  integer :: nlevs
116  logical :: hydrostatico
117  integer, allocatable :: id_tracer(:), all_axes(:)
118  integer, allocatable :: kstt_tracer(:), kend_tracer(:)
119  real, allocatable :: ak(:), bk(:)
120  character(20),allocatable :: axis_name(:),axis_name_vert(:)
121 
122  logical :: module_is_initialized=.false.
123  logical :: use_wrtgridcomp_output=.false.
124  integer :: sphum, liq_wat, ice_wat
125  integer :: rainwat, snowwat, graupel
126  real :: vrange(2) = (/ -330., 330. /)
127  real :: wrange(2) = (/ -100., 100. /)
128  real :: trange(2) = (/ 100., 350. /)
129  real :: skrange(2) = (/ -10000000.0, 10000000.0 /)
130 
131 ! file name
132  character(len=64) :: file_name = 'gfs_dyn'
133 
134 ! tracers
135  character(len=128) :: tname
136  character(len=256) :: tlongname, tunits
137 
138 ! wrtout buffer
139  real(4), dimension(:,:,:), allocatable, target :: buffer_dyn
140  real(4), dimension(:,:,:,:), allocatable, target :: windvect
141  real(4), dimension(:,:), allocatable, target :: psurf
142  real, dimension(:,:), allocatable :: lon, lat
143  real, dimension(:,:),allocatable :: up2,dn2,uhmax03,uhmin03
144  real, dimension(:,:),allocatable :: uhmax25,uhmin25,maxvort01
145  real, dimension(:,:),allocatable :: maxvorthy1,maxvort02
147 #ifdef use_WRTCOMP
148  public :: fv_dyn_bundle_setup
149 #endif
150 
151 contains
152 
153  subroutine fv_nggps_diag_init(Atm, axes, Time)
154  type(fv_atmos_type), intent(inout), target :: Atm(:)
155  integer, intent(in) :: axes(4)
156  type(time_type), intent(in) :: Time
157  integer :: n, i, j, nz
158 
159  n = 1
160  ncnsto = atm(1)%ncnst
161  npzo = atm(1)%npz
162  isco = atm(n)%bd%isc; ieco = atm(n)%bd%iec
163  jsco = atm(n)%bd%jsc; jeco = atm(n)%bd%jec
164  isdo = atm(n)%bd%isd; iedo = atm(n)%bd%ied
165  jsdo = atm(n)%bd%jsd; jedo = atm(n)%bd%jed
166  hydrostatico = atm(n)%flagstruct%hydrostatic
167 
168  call set_domain(atm(1)%domain) ! Set domain so that diag_manager can access tile information
169 
170  sphum = get_tracer_index(model_atmos, 'sphum')
171  liq_wat = get_tracer_index(model_atmos, 'liq_wat')
172  ice_wat = get_tracer_index(model_atmos, 'ice_wat')
173 
174  rainwat = get_tracer_index(model_atmos, 'rainwat')
175  snowwat = get_tracer_index(model_atmos, 'snowwat')
176  graupel = get_tracer_index(model_atmos, 'graupel')
177 
178 !--------------------------------------------------------------
179 ! Register main prognostic fields: ps, (u,v), t, omega (dp/dt)
180 !--------------------------------------------------------------
181  allocate(id_tracer(ncnsto))
182  allocate(kstt_tracer(ncnsto), kend_tracer(ncnsto))
183  id_tracer(:) = 0
184  kstt_tracer(:) = 0
185  kend_tracer(:) = 0
186 
187  if (atm(n)%flagstruct%write_3d_diags) then
188 !-------------------
189 ! A grid winds (lat-lon)
190 !-------------------
191  id_ua = register_diag_field( trim(file_name), 'ucomp', axes(1:3), time, &
192  'zonal wind', 'm/sec', missing_value=missing_value, range=vrange )
193  if (id_ua>0) then
194  kstt_ua = 1; kend_ua = npzo
195  nlevs = nlevs + npzo
196  endif
197 
198  id_va = register_diag_field( trim(file_name), 'vcomp', axes(1:3), time, &
199  'meridional wind', 'm/sec', missing_value=missing_value, range=vrange)
200  if (id_va>0) then
202  nlevs = nlevs + npzo
203  endif
204 
205  if(id_ua>0 .and. id_va>0) then
207  allocate(windvect(3,isco:ieco,jsco:jeco,npzo))
208  windvect = 0.
209  endif
210 
211  if( atm(n)%flagstruct%hydrostatic ) then
212  id_pfhy = register_diag_field( trim(file_name), 'pfhy', axes(1:3), time, &
213  'hydrostatic pressure', 'pa', missing_value=missing_value )
214  if (id_pfhy>0) then
216  nlevs = nlevs + npzo
217  endif
218  else
219  id_pfnh = register_diag_field( trim(file_name), 'pfnh', axes(1:3), time, &
220  'non-hydrostatic pressure', 'pa', missing_value=missing_value )
221  if (id_pfnh>0) then
223  nlevs = nlevs + npzo
224  endif
225  id_w = register_diag_field( trim(file_name), 'w', axes(1:3), time, &
226  'vertical wind', 'm/sec', missing_value=missing_value, range=wrange )
227  if (id_w>0) then
228  kstt_w = nlevs+1; kend_w = nlevs+npzo
229  nlevs = nlevs + npzo
230  endif
231  id_delz = register_diag_field( trim(file_name), 'delz', axes(1:3), time, &
232  'height thickness', 'm', missing_value=missing_value )
233  if (id_delz>0) then
235  nlevs = nlevs + npzo
236  endif
237  endif
238 
239  id_omga = register_diag_field( trim(file_name), 'omga', axes(1:3), time, &
240  'Vertical pressure velocity', 'pa/sec', missing_value=missing_value )
241  if (id_omga>0) then
243  nlevs = nlevs + npzo
244  endif
245 
246  id_pt = register_diag_field( trim(file_name), 'temp', axes(1:3), time, &
247  'temperature', 'K', missing_value=missing_value, range=trange )
248  if (id_pt>0) then
250  nlevs = nlevs + npzo
251  endif
252 
253  id_delp = register_diag_field( trim(file_name), 'delp', axes(1:3), time, &
254  'pressure thickness', 'pa', missing_value=missing_value )
255  if (id_delp>0) then
257  nlevs = nlevs + npzo
258  endif
259 
260  !--- diagnostic output for skeb: dissipation estimate
261  id_diss = register_diag_field( trim(file_name), 'diss_est', axes(1:3), time, &
262  'dissipation estimate', 'none', missing_value=missing_value, range=skrange )
263  if (id_delp>0) then
265  nlevs = nlevs + npzo
266  endif
267 
268 !--------------------
269 ! Tracer diagnostics:
270 !--------------------
271  do i=1, ncnsto
272  call get_tracer_names ( model_atmos, i, tname, tlongname, tunits )
273  id_tracer(i) = register_diag_field( file_name, trim(tname), &
274  axes(1:3), time, trim(tlongname), &
276  if (id_tracer(i)>0) then
277  kstt_tracer(i) = nlevs+1; kend_tracer(i) = nlevs+npzo
278  nlevs = nlevs + npzo
279  endif
280  enddo
281 !
282  id_ps = register_diag_field( trim(file_name), 'ps', axes(1:2), time, &
283  'surface pressure', 'pa', missing_value=missing_value )
284  if( id_ps > 0) then
285  kstt_ps = nlevs+1; kend_ps = nlevs+1
286  nlevs = nlevs + 1
287  allocate(psurf(isco:ieco,jsco:jeco))
288  endif
289 !
290  id_hs = register_diag_field( trim(file_name), 'hs', axes(1:2), time, &
291  'surface geopotential height', 'gpm', missing_value=missing_value )
292  if( id_hs > 0) then
293  kstt_hs = nlevs+1; kend_hs = nlevs+1
294  nlevs = nlevs + 1
295  endif
296 !
297  id_dbz = register_diag_field( trim(file_name), 'reflectivity', axes(1:3), time, &
298  'Stoelinga simulated reflectivity', 'dBz', missing_value=missing_value)
299  if( rainwat > 0 .and. id_dbz > 0) then
301  nlevs = nlevs + npzo
302  endif
303  id_ustm = register_diag_field( trim(file_name), 'ustm',axes(1:2), time, &
304  'u comp of storm motion', 'm/s', missing_value=missing_value )
305  if( id_ustm > 0) then
306  kstt_ustm = nlevs+1; kend_ustm = nlevs+1
307  nlevs = nlevs + 1
308  endif
309  id_vstm = register_diag_field( trim(file_name), 'vstm',axes(1:2), time, &
310  'v comp of storm motion', 'm/s', missing_value=missing_value )
311  if( id_vstm > 0) then
312  kstt_vstm = nlevs+1; kend_vstm = nlevs+1
313  nlevs = nlevs + 1
314  endif
315 
316  id_srh01 = register_diag_field( trim(file_name), 'srh01',axes(1:2), time, &
317  '0-1km storm rel. helicity', 'm/s**2', missing_value=missing_value )
318  if( id_srh01 > 0) then
320  nlevs = nlevs + 1
321  endif
322  id_srh03 = register_diag_field( trim(file_name), 'srh03',axes(1:2), time, &
323  '0-3km storm rel. helicity', 'm/s**2', missing_value=missing_value )
324  if( id_srh03 > 0) then
326  nlevs = nlevs + 1
327  endif
328  id_maxvort01 = register_diag_field( trim(file_name), 'maxvort01',axes(1:2), time, &
329  'Max hourly 0-1km vert vorticity', '1/s', missing_value=missing_value )
330  if( id_maxvort01 > 0) then
331  allocate ( maxvort01(isco:ieco,jsco:jeco) )
333  nlevs = nlevs + 1
334  endif
335  id_maxvort02 = register_diag_field( trim(file_name), 'maxvort02',axes(1:2), time, &
336  'Max hourly 0-2km vert vorticity', '1/s', missing_value=missing_value )
337  if( id_maxvort02 > 0) then
338  allocate ( maxvort02(isco:ieco,jsco:jeco) )
340  nlevs = nlevs + 1
341  endif
342  id_maxvorthy1 = register_diag_field( trim(file_name), 'maxvorthy1',axes(1:2), time, &
343  'Max hourly hybrid lev1 vert. vorticity', '1/s', missing_value=missing_value )
344  if( id_maxvorthy1 > 0) then
345  allocate ( maxvorthy1(isco:ieco,jsco:jeco) )
347  nlevs = nlevs + 1
348  endif
349  id_wmaxup = register_diag_field( trim(file_name), 'wmaxup',axes(1:2), time, &
350  'Max hourly updraft velocity', 'm/s', missing_value=missing_value )
351  if( id_wmaxup > 0) then
352  allocate ( up2(isco:ieco,jsco:jeco) )
353  kstt_wup = nlevs+1; kend_wup = nlevs+1
354  nlevs = nlevs + 1
355  endif
356  id_wmaxdn = register_diag_field( trim(file_name), 'wmaxdn',axes(1:2), time, &
357  'Max hourly downdraft velocity', 'm/s', missing_value=missing_value )
358 ! write (0,*)'id_wmaxdn in fv_nggps=',id_wmaxdn
359  if( id_wmaxdn > 0) then
360  allocate ( dn2(isco:ieco,jsco:jeco) )
361  kstt_wdn = nlevs+1; kend_wdn = nlevs+1
362  nlevs = nlevs + 1
363  endif
364  id_uhmax03 = register_diag_field( trim(file_name), 'uhmax03',axes(1:2), time, &
365  'Max hourly max 0-3km updraft helicity', 'm/s**2', missing_value=missing_value )
366 ! write (0,*)'id_uhmax03 in fv_nggps=',id_uhmax03
367  if( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmax03 > 0 ) then
368  allocate ( uhmax03(isco:ieco,jsco:jeco) )
370  nlevs = nlevs + 1
371  endif
372 !
373  id_uhmin03 = register_diag_field( trim(file_name), 'uhmin03',axes(1:2), time, &
374  'Max hourly min 0-3km updraft helicity', 'm/s**2', missing_value=missing_value )
375  if( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmin03 > 0 ) then
376  allocate ( uhmin03(isco:ieco,jsco:jeco) )
378  nlevs = nlevs + 1
379  endif
380 !
381  id_uhmax25 = register_diag_field( trim(file_name), 'uhmax25',axes(1:2), time, &
382  'Max hourly max 2-5km updraft helicity', 'm/s**2', missing_value=missing_value )
383  if( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmax25 > 0 ) then
384  allocate ( uhmax25(isco:ieco,jsco:jeco) )
386  nlevs = nlevs + 1
387  endif
388 !
389  id_uhmin25 = register_diag_field( trim(file_name), 'uhmin25',axes(1:2), time, &
390  'Max hourly min 2-5km updraft helicity', 'm/s**2', missing_value=missing_value )
391  if( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmin25 > 0 ) then
392  allocate ( uhmin25(isco:ieco,jsco:jeco) )
394  nlevs = nlevs + 1
395  endif
396 !
397  nz = size(atm(1)%ak)
398  allocate(ak(nz))
399  allocate(bk(nz))
400  do i=1,nz
401  ak(i) = atm(1)%ak(i)
402  bk(i) = atm(1)%bk(i)
403  enddo
404 ! print *,'in ngpps diag init, ak=',ak(1:5),' bk=',bk(1:5)
405 
406 ! get lon,lat information
407  if(.not.allocated(lon)) then
408  allocate(lon(isco:ieco,jsco:jeco))
409  do j=jsco,jeco
410  do i=isco,ieco
411  lon(i,j) = atm(n)%gridstruct%agrid(i,j,1)
412  enddo
413  enddo
414 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,lon=',lon(isco,jsco),lon(ieco-2:ieco,jeco-2:jeco)*180./3.14157
415  endif
416  if(.not.allocated(lat)) then
417  allocate(lat(isco:ieco,jsco:jeco))
418  do j=jsco,jeco
419  do i=isco,ieco
420  lat(i,j) = atm(n)%gridstruct%agrid(i,j,2)
421  enddo
422  enddo
423 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,lat=',lat(isco,jsco),lat(ieco-2:ieco,jeco-2:jeco)*180./3.14157
424  endif
425  endif
426 !
427 !------------------------------------
428 ! use wrte grid component for output
429 !------------------------------------
430  use_wrtgridcomp_output = .false.
431 
432  end subroutine fv_nggps_diag_init
433 
434 
435  subroutine fv_nggps_diag(Atm, zvir, Time)
437  type(fv_atmos_type), intent(inout) :: Atm(:)
438  real, intent(in):: zvir
439  type(time_type), intent(in) :: Time
440 
441  integer :: i, j, k, n, ngc, nq, itrac
442  logical :: bad_range
443  real :: ptop, allmax
444  real, allocatable :: wk(:,:,:), wk2(:,:,:)
445  real, dimension(:,:),allocatable :: ustm,vstm,srh01,srh03
446 
447  n = 1
448  ngc = atm(n)%ng
449  ptop = atm(n)%ak(1)
450  allmax = -20.
451  nq = size (atm(n)%q,4)
452  allocate ( wk(isco:ieco,jsco:jeco,npzo) )
453  allocate ( wk2(isco:ieco,jsco:jeco,npzo) )
454  allocate ( ustm(isco:ieco,jsco:jeco) )
455  allocate ( vstm(isco:ieco,jsco:jeco) )
456  allocate ( srh01(isco:ieco,jsco:jeco) )
457  allocate ( srh03(isco:ieco,jsco:jeco) )
458  if ( atm(n)%flagstruct%range_warn ) then
459  call range_check('DELP', atm(n)%delp, isco, ieco, jsco, jeco, ngc, npzo, atm(n)%gridstruct%agrid, &
460  0.01*ptop, 200.e2, bad_range)
461  call range_check('UA', atm(n)%ua, isco, ieco, jsco, jeco, ngc, npzo, atm(n)%gridstruct%agrid, &
462  -250., 250., bad_range)
463  call range_check('VA', atm(n)%va, isco, ieco, jsco, jeco, ngc, npzo, atm(n)%gridstruct%agrid, &
464  -250., 250., bad_range)
465  call range_check('TA', atm(n)%pt, isco, ieco, jsco, jeco, ngc, npzo, atm(n)%gridstruct%agrid, &
466  150., 350., bad_range) !DCMIP ICs have very low temperatures
467  endif
468  !--- A-GRID WINDS
469  if ( .not. allocated(buffer_dyn)) allocate(buffer_dyn(isco:ieco,jsco:jeco,nlevs))
470  if(id_ua > 0) call store_data(id_ua, atm(n)%ua(isco:ieco,jsco:jeco,:), time, kstt_ua, kend_ua)
471 
472  if(id_va > 0) call store_data(id_va, atm(n)%va(isco:ieco,jsco:jeco,:), time, kstt_va, kend_va)
473 
474  !--- set up 3D wind vector
475  if(id_ua>0 .and. id_va>0) then
476  do k=1,npzo
477  do j=jsco,jeco
478  do i=isco,ieco
479  windvect(1,i,j,k) = atm(n)%ua(i,j,k)*cos(lon(i,j)) - atm(n)%va(i,j,k)*sin(lat(i,j))*sin(lon(i,j))
480  windvect(2,i,j,k) = atm(n)%ua(i,j,k)*sin(lon(i,j)) + atm(n)%va(i,j,k)*sin(lat(i,j))*cos(lon(i,j))
481  windvect(3,i,j,k) = atm(n)%va(i,j,k)*cos(lat(i,j))
482  enddo
483  enddo
484  enddo
485  endif
486 
487  !--- W (non-hydrostatic)
488  if ( .not.atm(n)%flagstruct%hydrostatic .and. id_w>0 ) then
489  call store_data(id_w, atm(n)%w(isco:ieco,jsco:jeco,:), time, kstt_w, kend_w)
490  endif
491 
492  !--- OMGA (non-hydrostatic)
493  if ( id_omga>0 ) then
494  call store_data(id_omga, atm(n)%omga(isco:ieco,jsco:jeco,:), time, kstt_omga, kend_omga)
495  endif
496 
497  !--- TEMPERATURE
498  if(id_pt > 0) call store_data(id_pt, atm(n)%pt(isco:ieco,jsco:jeco,:), time, kstt_pt, kend_pt)
499 
500  !--- TRACERS
501  do itrac=1, ncnsto
502  call get_tracer_names (model_atmos, itrac, tname)
503  if (id_tracer(itrac) > 0 .and. itrac.gt.nq) then
504  call store_data (id_tracer(itrac), atm(n)%qdiag(isco:ieco,jsco:jeco,:,itrac), time, &
505  kstt_tracer(itrac),kend_tracer(itrac) )
506  else
507  call store_data (id_tracer(itrac), atm(n)%q(isco:ieco,jsco:jeco,:,itrac), time, &
508  kstt_tracer(itrac),kend_tracer(itrac) )
509  endif
510  enddo
511 
512  !--- DELZ (non-hydrostatic)
513  if((.not. atm(n)%flagstruct%hydrostatic) .and. id_delz > 0) then
514  do k=1,npzo
515  do j=jsco,jeco
516  do i=isco,ieco
517  wk(i,j,k) = atm(n)%delz(i,j,k)
518  enddo
519  enddo
520  enddo
521  call store_data(id_delz, wk, time, kstt_delz, kend_delz)
522  endif
523 
524  !--- PRESSURE (hydrostatic)
525  if( atm(n)%flagstruct%hydrostatic .and. id_pfhy > 0 ) then
526  do k=1,npzo
527  do j=jsco,jeco
528  do i=isco,ieco
529  wk(i,j,k) = 0.5 *(atm(n)%pe(i,k,j)+atm(n)%pe(i,k+1,j))
530  enddo
531  enddo
532  enddo
533  call store_data(id_pfhy, wk, time, kstt_pfhy, kend_pfhy)
534  endif
535 
536 #ifdef GFS_PHYS
537  !--- DELP
538  if(id_delp > 0 .or. ((.not. atm(n)%flagstruct%hydrostatic) .and. id_pfnh > 0)) then
539  do k=1,npzo
540  do j=jsco,jeco
541  do i=isco,ieco
542  wk(i,j,k) = atm(n)%delp(i,j,k)*(1.-sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
543  enddo
544  enddo
545  enddo
546  call store_data(id_delp, wk, time, kstt_delp, kend_delp)
547  endif
548 
549  !--- Surface Pressure (PS)
550  ! Re-compute pressure (dry_mass + water_vapor) surface pressure
551  if(id_ps > 0) then
552  do k=1,npzo
553  do j=jsco,jeco
554  do i=isco,ieco
555  wk(i,j,k) = atm(n)%delp(i,j,k)*(1.-sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
556  enddo
557  enddo
558  enddo
559  do j=jsco,jeco
560  do i=isco,ieco
561  psurf(i,j) = ptop
562  do k=npzo,1,-1
563  psurf(i,j) = psurf(i,j) + wk(i,j,k)
564  enddo
565  enddo
566  enddo
567  endif
568 
569  !--- PRESSURE (non-hydrostatic)
570  if( (.not. atm(n)%flagstruct%hydrostatic) .and. id_pfnh > 0) then
571  do k=1,npzo
572  do j=jsco,jeco
573  do i=isco,ieco
574  wk(i,j,k) = -wk(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*atm(n)%pt(i,j,k)
575 #ifdef MULTI_GASES
576  wk(i,j,k) = wk(i,j,k) * virq(atm(n)%q(i,j,k,:))
577 #else
578  wk(i,j,k) = wk(i,j,k) * (1.+zvir*atm(n)%q(i,j,k,sphum))
579 #endif
580  enddo
581  enddo
582  enddo
583  call store_data(id_pfnh, wk, time, kstt_pfnh, kend_pfnh)
584  endif
585 #else
586  !--- DELP
587  if(id_delp > 0) call store_data(id_delp, atm(n)%delp(isco:ieco,jsco:jeco,:), time, kstt_delp)
588 
589  !--- Surface Pressure (PS)
590  if( id_ps > 0) then
591  do j=jsco,jeco
592  do i=isco,ieco
593  psurf(i,j) = atm(n)%ps(i,j)
594  enddo
595  enddo
596  endif
597 
598  !--- PRESSURE (non-hydrostatic)
599  if( (.not. atm(n)%flagstruct%hydrostatic) .and. id_pfnh > 0) then
600  do k=1,npzo
601  do j=jsco,jeco
602  do i=isco,ieco
603  wk(i,j,k) = -atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas*atm(n)%pt(i,j,k)
604 #ifdef MULTI_GASES
605  wk(i,j,k) = wk(i,j,k)*virq(atm(n)%q(i,j,k,:)
606 #else
607  wk(i,j,k) = wk(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,sphum))
608 #endif
609  enddo
610  enddo
611  enddo
612  call store_data(id_pfnh, wk, time, kstt_pfnh, kend_pfnh)
613  endif
614 #endif
615 
616  !--- DISS_EST (skeb: dissipation estimate)
617  if(id_diss > 0) call store_data(id_diss, atm(n)%diss_est(isco:ieco,jsco:jeco,:), time, kstt_diss, kend_diss)
618 !
619  if(id_ps > 0) then
620  if( use_wrtgridcomp_output ) then
621  do j=jsco,jeco
622  do i=isco,ieco
623  wk(i,j,1) = (psurf(i,j)/stndrd_atmos_ps)**(rdgas/grav*stndrd_atmos_lapse)
624  enddo
625  enddo
626  else
627  do j=jsco,jeco
628  do i=isco,ieco
629  wk(i,j,1) = psurf(i,j)
630  enddo
631  enddo
632  endif
633 ! print *,'in comput ps, i=',isco,'j=',jsco,'psurf=',psurf(isco,jsco),'stndrd_atmos_ps=',stndrd_atmos_ps, &
634 ! 'rdgas=',rdgas,'grav=',grav,'stndrd_atmos_lapse=',stndrd_atmos_lapse,rdgas/grav*stndrd_atmos_lapse
635  call store_data(id_ps, wk, time, kstt_ps, kend_ps)
636  endif
637 
638  if( id_hs > 0) then
639  do j=jsco,jeco
640  do i=isco,ieco
641  wk(i,j,1) = atm(n)%phis(i,j)/grav
642  enddo
643  enddo
644  call store_data(id_hs, wk, time, kstt_hs, kend_hs)
645  endif
646 
647  !--- 3-D Reflectivity field
648  if ( rainwat > 0 .and. id_dbz>0) then
649  call dbzcalc(atm(n)%q, atm(n)%pt, atm(n)%delp, atm(n)%peln, atm(n)%delz, &
650  wk, wk2, allmax, atm(n)%bd, npzo, atm(n)%ncnst, atm(n)%flagstruct%hydrostatic, &
651  zvir, .false., .false., .false., .true. ) ! GFDL MP has constant N_0 intercept
652  call store_data(id_dbz, wk, time, kstt_dbz, kend_dbz)
653  endif
654 
655  deallocate ( wk )
656  !---u and v comp of storm motion, 0-1, 0-3km SRH
657  if ( id_ustm > 0 .or. id_vstm > 0 .or. id_srh01 > 0 .or. id_srh03 > 0) then
658  if ( id_ustm > 0 .and. id_vstm > 0 .and. id_srh01 > 0 .and. id_srh03 > 0) then
659  call bunkers_vector(isco,ieco,jsco,jeco,ngc,npzo,zvir,sphum,ustm,vstm, &
660  atm(n)%ua,atm(n)%va, atm(n)%delz, atm(n)%q, &
661  atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln,&
662  atm(n)%phis, grav)
663 
664  call helicity_relative_caps(isco,ieco,jsco,jeco,ngc,npzo,zvir,sphum,srh01, &
665  ustm, vstm,atm(n)%ua, atm(n)%va, atm(n)%delz, &
666  atm(n)%q,atm(n)%flagstruct%hydrostatic, &
667  atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 1.e3)
668 
669  call helicity_relative_caps(isco,ieco,jsco,jeco,ngc,npzo,zvir,sphum,srh03, &
670  ustm, vstm,atm(n)%ua, atm(n)%va, atm(n)%delz, &
671  atm(n)%q,atm(n)%flagstruct%hydrostatic, &
672  atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 3.e3)
673 
674  call store_data(id_ustm, ustm, time, kstt_ustm, kend_ustm)
675  call store_data(id_vstm, vstm, time, kstt_vstm, kend_vstm)
676  call store_data(id_srh01, srh01, time, kstt_srh01, kend_srh01)
677  call store_data(id_srh03, srh03, time, kstt_srh03, kend_srh03)
678  else
679  print *,'Missing fields in diag_table'
680  print *,'Make sure the following are listed in the diag_table under gfs_dyn:'
681  print *,'ustm,vstm,srh01,shr03'
682  call mpp_error(fatal, 'Missing fields in diag_table')
683  stop
684  endif
685  endif
686  deallocate ( ustm )
687  deallocate ( vstm )
688  deallocate ( srh01 )
689  deallocate ( srh03 )
690 
691  !--- max hourly 0-1km vert. vorticity
692  if ( id_maxvort01 > 0) then
694  endif
695  !--- max hourly 0-2km vert. vorticity
696  if ( id_maxvort02 > 0) then
698  endif
699  !--- max hourly hybrid lev 1 vert. vorticity
700  if ( id_maxvorthy1 > 0) then
702  endif
703 !
704  !--- max hourly updraft velocity
705  if ( id_wmaxup > 0) then
706  call store_data(id_wmaxup, up2, time, kstt_wup, kend_wup)
707  endif
708  !--- max hourly downdraft velocity
709  if ( id_wmaxdn > 0) then
710  call store_data(id_wmaxdn, dn2, time, kstt_wdn, kend_wdn)
711  endif
712  !--- max hourly max 0-3km updraft helicity
713  if ( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmax03 > 0) then
715  endif
716 !
717  !--- max hourly min 0-3km updraft helicity
718  if ( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmin03 > 0) then
720  endif
721 !
722  !--- max hourly max 2-5km updraft helicity
723  if ( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmax25 > 0) then
725  endif
726 !
727  !--- max hourly min 2-5km updraft helicity
728  if ( .not.atm(n)%flagstruct%hydrostatic .and. id_uhmin25 > 0) then
730  endif
731 
732  call nullify_domain()
733 
734  end subroutine fv_nggps_diag
735 
736  subroutine fv_nggps_tavg(Atm, Time_step_atmos,avg_max_length,zvir)
737  type(fv_atmos_type), intent(inout) :: Atm(:)
738  type(time_type), intent(in) :: Time_step_atmos
739  real, intent(in):: zvir
740  integer :: i, j, k, n, ngc, nq, itrac
741  integer seconds, days, nsteps_per_reset
742  logical, save :: first_call=.true.
743  real, save :: first_time = 0.
744  integer, save :: kdtt = 0
745  real :: avg_max_length
746  real,dimension(:,:,:),allocatable :: vort
747  n = 1
748  ngc = atm(n)%ng
749  nq = size (atm(n)%q,4)
750 !
751 !Check if any of the max hourly fields are being requested otherwise skip
752 !
753  if(id_wmaxup > 0 .or. id_wmaxdn > 0 .or. id_uhmax03 > 0 .or. id_uhmin03 > 0 &
754  .or. id_uhmax25 > 0 .or. id_uhmin25 > 0 .or. id_maxvort01 > 0 &
755  .or. id_maxvorthy1 > 0 .or. id_maxvort02 > 0) then
756 !Make sure the group of max hrly fields listed below are ALL present otherwise
757 !abort
758 !
759  if(id_wmaxup > 0 .and. id_wmaxdn > 0 .and. id_uhmax03 > 0 .and. id_uhmin03 > 0 &
760  .and. id_uhmax25 > 0 .and. id_uhmin25 > 0 .and. id_maxvort01 > 0 &
761  .and. id_maxvorthy1 > 0 .and. id_maxvort02 > 0) then
762  allocate ( vort(isco:ieco,jsco:jeco,npzo) )
763  if (first_call) then
764  call get_time (time_step_atmos, seconds, days)
765  first_time=seconds
766  first_call=.false.
767  kdtt=0
768  endif
769  nsteps_per_reset = nint(avg_max_length/first_time)
770  do j=jsco,jeco
771  do i=isco,ieco
772  if(mod(kdtt,nsteps_per_reset)==0)then
773  up2(i,j) = -999.
774  dn2(i,j) = 999.
775  maxvorthy1(i,j)= 0.
776  maxvort01(i,j)= 0.
777  maxvort02(i,j)= 0.
778  endif
779  enddo
780  enddo
782  npzo,atm(n)%u,atm(n)%v,vort, &
783  atm(n)%gridstruct%dx,atm(n)%gridstruct%dy,&
784  atm(n)%gridstruct%rarea)
786  call max_vorticity(isco,ieco,jsco,jeco,ngc,npzo,zvir, &
787  sphum,atm(n)%delz,atm(n)%q, &
788  atm(n)%flagstruct%hydrostatic, &
789  atm(n)%pt,atm(n)%peln,atm(n)%phis,grav, &
790  vort,maxvort01,0., 1.e3)
791  call max_vorticity(isco,ieco,jsco,jeco,ngc,npzo,zvir, &
792  sphum,atm(n)%delz,atm(n)%q, &
793  atm(n)%flagstruct%hydrostatic, &
794  atm(n)%pt,atm(n)%peln,atm(n)%phis,grav, &
795  vort,maxvort02,0., 2.e3)
796  if( .not.atm(n)%flagstruct%hydrostatic ) then
797  call max_vv(isco,ieco,jsco,jeco,npzo,ngc,up2,dn2,atm(n)%pe,atm(n)%w)
798  do j=jsco,jeco
799  do i=isco,ieco
800  if(mod(kdtt,nsteps_per_reset)==0)then
801  uhmax03(i,j)= 0.
802  uhmin03(i,j)= 0.
803  uhmax25(i,j)= 0.
804  uhmin25(i,j)= 0.
805  endif
806  enddo
807  enddo
808 
809  call max_uh(isco,ieco,jsco,jeco,ngc,npzo,zvir, &
810  sphum,uhmax03,uhmin03,atm(n)%w,vort,atm(n)%delz, &
811  atm(n)%q,atm(n)%flagstruct%hydrostatic, &
812  atm(n)%pt,atm(n)%peln,atm(n)%phis,grav, &
813  0., 3.e3)
814  call max_uh(isco,ieco,jsco,jeco,ngc,npzo,zvir, &
815  sphum,uhmax25,uhmin25,atm(n)%w,vort,atm(n)%delz, &
816  atm(n)%q,atm(n)%flagstruct%hydrostatic, &
817  atm(n)%pt,atm(n)%peln,atm(n)%phis,grav, &
818  2.e3, 5.e3)
819  endif
820  kdtt=kdtt+1
821  deallocate (vort)
822  else
823  print *,'Missing max/min hourly field in diag_table'
824  print *,'Make sure the following are listed in the diag_table under gfs_dyn:'
825  print *,'wmaxup,wmaxdn,uhmax03,uhmin03,uhmax25,uhmin25,maxvort01,maxvort02 and maxvorthy1'
826  call mpp_error(fatal, 'Missing max hourly fields in diag_table')
827  stop
828  endif
829  endif
830  end subroutine fv_nggps_tavg
831 !
832  subroutine store_data(id, work, Time, nstt, nend)
833  integer, intent(in) :: id
834  integer, intent(in) :: nstt, nend
835  real, intent(in) :: work(isco:ieco,jsco:jeco,nend-nstt+1)
836  type(time_type), intent(in) :: Time
837 !
838  integer k,j,i,kb
839  logical :: used
840 !
841  if( id > 0 ) then
842  if( use_wrtgridcomp_output ) then
843  do k=1,nend-nstt+1
844  do j= jsco,jeco
845  do i=isco,ieco
846  kb = k + nstt - 1
847  buffer_dyn(i,j,kb) = work(i,j,k)
848  enddo
849  enddo
850  enddo
851  else
852  used = send_data(id, work, time)
853  endif
854  endif
855 
856  end subroutine store_data
857 
858 #ifdef use_WRTCOMP
859 
860  subroutine fv_dyn_bundle_setup(axes, dyn_bundle, fcst_grid, quilting, rc)
861 !
862 !-------------------------------------------------------------
863 !*** set esmf bundle for dyn output fields
864 !------------------------------------------------------------
865 !
866  use esmf
867  use diag_data_mod, ONLY: diag_atttype
868 !
869  integer, intent(in) :: axes(:)
870  type(esmf_fieldbundle),intent(inout) :: dyn_bundle
871  type(esmf_grid),intent(inout) :: fcst_grid
872  logical,intent(in) :: quilting
873  integer,intent(out) :: rc
874 
875 
876 !*** local variables
877  integer i, j, k, n
878  integer num_axes, id, axis_length, direction, edges
879  integer num_attributes, num_field_dyn, axis_typ
880  character(255) :: units, long_name, cart_name,axis_direct,edgesS
881  character(128) :: output_name, output_file, output_file1, dynbdl_name, shydrostatic
882  integer currdate(6), idx1
883  logical l3Dvector
884  type(domain1d) :: Domain
885  type(domainug) :: DomainU
886  real,dimension(:),allocatable :: axis_data
887  type(diag_atttype),dimension(:),allocatable :: attributes
888  character(2) axis_id
889 
890  type(esmf_field) :: field
891 !
892 !jwtest
893 ! integer :: fieldcount
894 ! character(128) :: fld_outfilename
895 ! character(128),dimension(:),allocatable :: fieldnamelist
896 ! type(ESMF_Field),dimension(:),allocatable :: fieldlist
897 !
898 !------------------------------------------------------------
899 
900 ! initialize RC
901  rc = esmf_success
902 
903 !--- use wrte grid component for output
904  use_wrtgridcomp_output = quilting
905 
906 ! data type
907  if(.not. allocated(buffer_dyn))allocate(buffer_dyn(isco:ieco,jsco:jeco,nlevs))
908  buffer_dyn=0.
909  num_field_dyn = 0.
910 !
911 ! set output files
912  call esmf_fieldbundleget(dyn_bundle, name=dynbdl_name,rc=rc)
913  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
914  line=__line__, &
915  file=__file__)) &
916  return ! bail out
917  idx1 = index(dynbdl_name,'_bilinear')
918  if(idx1 > 0) then
919  output_file = dynbdl_name(1:idx1-1)
920  else
921  output_file = 'dyn'
922  endif
923 !
924 !------------------------------------------------------------
925 !*** add attributes to the bundle such as subdomain limtis,
926 !*** axes, output time, etc
927 !------------------------------------------------------------
928 !
929 !*** add attributes
930  num_axes = size(axes)
931  allocate(all_axes(num_axes))
932  all_axes(1:num_axes) = axes(1:num_axes)
933 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,num_axes=',num_axes, 'axes=',axes
934 !
935 !*** add global attributes in the field bundle:
936  call esmf_attributeadd(dyn_bundle, convention="NetCDF", purpose="FV3", &
937  attrlist=(/"hydrostatic", &
938  "ncnsto ", &
939  "ak ", &
940  "bk "/), rc=rc)
941  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
942  line=__line__, &
943  file=__file__)) &
944  return ! bail out
945  if (hydrostatico ) then
946  shydrostatic = 'hydrostatic'
947  else
948  shydrostatic = 'non-hydrostatic'
949  endif
950  call esmf_attributeset(dyn_bundle, convention="NetCDF", purpose="FV3", &
951  name="hydrostatic", value=trim(shydrostatic), rc=rc)
952  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
953  line=__line__, &
954  file=__file__)) &
955  return ! bail out
956 !
957  call esmf_attributeset(dyn_bundle, convention="NetCDF", purpose="FV3", &
958  name="ncnsto", value=ncnsto, rc=rc)
959  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
960  line=__line__, &
961  file=__file__)) &
962  return ! bail out
963 !
964  call esmf_attributeset(dyn_bundle, convention="NetCDF", purpose="FV3", &
965  name="ak", valuelist=ak, rc=rc)
966 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,after add ak, rc=',rc
967  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
968  line=__line__, &
969  file=__file__)) &
970  return ! bail out
971 !
972  call esmf_attributeset(dyn_bundle, convention="NetCDF", purpose="FV3", &
973  name="bk", valuelist=bk, rc=rc)
974 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,after add bk, rc=',rc
975  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
976  line=__line__, &
977  file=__file__)) &
978  return ! bail out
979 !
980 !*** get axis names
981  allocate(axis_name(num_axes))
982  do id = 1,num_axes
983  call get_diag_axis_name( axes(id), axis_name(id))
984  enddo
985  if( num_axes>2 ) then
986  allocate(axis_name_vert(num_axes-2))
987  do id=3,num_axes
988  axis_name_vert(id-2) = axis_name(id)
989  enddo
990 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,axis_name_vert=',axis_name_vert
991  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
992  attrlist=(/"vertical_dim_labels"/), rc=rc)
993  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
994  line=__line__, &
995  file=__file__)) &
996  return ! bail out
997  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
998  name="vertical_dim_labels", valuelist=axis_name_vert, rc=rc)
999  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1000  line=__line__, &
1001  file=__file__)) &
1002  return ! bail out
1003  endif
1004 
1005  do id = 1,num_axes
1006  axis_length = get_axis_global_length(axes(id))
1007  allocate(axis_data(axis_length))
1008  call get_diag_axis( axes(id), axis_name(id), units, long_name, cart_name, &
1009  direction, edges, domain, domainu, axis_data, &
1010  num_attributes=num_attributes, &
1011  attributes=attributes)
1012 !
1013  edgess=''
1014  do i = 1,num_axes
1015  if(axes(i) == edges) edgess=axis_name(i)
1016  enddo
1017 
1018 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,id=',id,'edges=',edges,rc, &
1019 ! 'num_attributes=',num_attributes,'edgesS=',trim(edgesS)
1020 !
1021 ! Add vertical dimension Attributes to Grid
1022  if( id>2 ) then
1023 ! if(mpp_pe()==mpp_root_pe())print *,' in dyn add grid, axis_name=', &
1024 ! trim(axis_name(id)),'axis_data=',axis_data
1025 !
1026 ! Previous definition using variable-length character arrays violates the Fortran standards.
1027 ! While this worked with Intel compilers, it caused the model to crash in different places
1028 ! with both GNU and PGI. Compilers should throw an error at compile time, but it seems that
1029 ! they can't handle the "trim(...) // ..." expressions.
1030 ! The Standard (Fortran 2003) way to do this correctly is to tell the array constructor
1031 ! how long to make the fixed array of characters:
1032 !
1033 ! call ESMF_AttributeAdd(fcst_grid, convention="NetCDF", purpose="FV3", &
1034 ! attrList=(/ character(128) :: trim(axis_name(id)),trim(axis_name(id))//":long_name", &
1035 ! trim(axis_name(id))//":units", trim(axis_name(id))//":cartesian_axis", &
1036 ! trim(axis_name(id))//":positive", trim(axis_name(id))//":edges"/), rc=rc)
1037 !
1038 ! However this fails for GNU and PGI, see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=85547
1039 ! The easiest and safest way forward is to define the attributes one by one as it is done
1040 ! as it is done below in add_field_to_bundle.
1041 !
1042  ! Add attributes one by one
1043  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1044  attrlist=(/trim(axis_name(id))/), rc=rc)
1045  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1046  line=__line__, &
1047  file=__file__)) &
1048  return ! bail out
1049  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1050  attrlist=(/trim(axis_name(id))//":long_name"/), rc=rc)
1051  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1052  line=__line__, &
1053  file=__file__)) &
1054  return ! bail out
1055  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1056  attrlist=(/trim(axis_name(id))//":units"/), rc=rc)
1057  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1058  line=__line__, &
1059  file=__file__)) &
1060  return ! bail out
1061  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1062  attrlist=(/trim(axis_name(id))//":cartesian_axis"/), rc=rc)
1063  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1064  line=__line__, &
1065  file=__file__)) &
1066  return ! bail out
1067  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1068  attrlist=(/trim(axis_name(id))//":positive"/), rc=rc)
1069  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1070  line=__line__, &
1071  file=__file__)) &
1072  return ! bail out
1073  if(trim(edgess)/='') then
1074  call esmf_attributeadd(fcst_grid, convention="NetCDF", purpose="FV3", &
1075  attrlist=(/trim(axis_name(id))//":edges"/), rc=rc)
1076  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1077  line=__line__, &
1078  file=__file__)) &
1079  return ! bail out
1080  endif
1081  ! Set attributes
1082  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1083  name=trim(axis_name(id)), valuelist=axis_data, rc=rc)
1084  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1085  line=__line__, &
1086  file=__file__)) &
1087  return ! bail out
1088  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1089  name=trim(axis_name(id))//":long_name", value=trim(long_name), rc=rc)
1090  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1091  line=__line__, &
1092  file=__file__)) &
1093  return ! bail out
1094  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1095  name=trim(axis_name(id))//":units", value=trim(units), rc=rc)
1096  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1097  line=__line__, &
1098  file=__file__)) &
1099  return ! bail out
1100  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1101  name=trim(axis_name(id))//":cartesian_axis", value=trim(cart_name), rc=rc)
1102  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1103  line=__line__, &
1104  file=__file__)) &
1105  return ! bail out
1106  if(direction>0) then
1107  axis_direct="up"
1108  else
1109  axis_direct="down"
1110  endif
1111  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1112  name=trim(axis_name(id))//":positive", value=trim(axis_direct), rc=rc)
1113  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1114  line=__line__, &
1115  file=__file__)) &
1116  return ! bail out
1117  if(trim(edgess)/='') then
1118  call esmf_attributeset(fcst_grid, convention="NetCDF", purpose="FV3", &
1119  name=trim(axis_name(id))//":edges", value=trim(edgess), rc=rc)
1120  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1121  line=__line__, &
1122  file=__file__)) &
1123  return ! bail out
1124  endif
1125 
1126  endif
1127 
1128  deallocate(axis_data)
1129  enddo
1130 !
1131 !*** add esmf fields
1132  if(id_ua > 0) then
1133  call find_outputname(trim(file_name),'ucomp',output_name)
1134 ! if(mpp_pe()==mpp_root_pe()) print *,'ucomp output name is ',trim(output_name)
1135  call add_field_to_bundle(trim(output_name),'zonal wind', 'm/sec', "time: point", &
1136  axes(1:3), fcst_grid, kstt_ua,kend_ua, dyn_bundle, output_file, &
1137  range=vrange, rcd=rc)
1138  if(rc==0) num_field_dyn=num_field_dyn+1
1139  endif
1140 !
1141  if(id_va > 0) then
1142  call find_outputname(trim(file_name),'vcomp',output_name)
1143  call add_field_to_bundle(trim(output_name),'meridional wind', 'm/sec', "time: point", &
1144  axes(1:3), fcst_grid, kstt_va,kend_va, dyn_bundle, output_file, &
1145  range=vrange,rcd=rc)
1146  if(rc==0) num_field_dyn=num_field_dyn+1
1147  endif
1148 !
1149 !*** create 3D vector from local u/v winds
1150  if(id_ua > 0 .and. id_va > 0) then
1151  output_name = "windvector"
1152  output_file1 = 'none'
1153  l3dvector = .true.
1154  call add_field_to_bundle(trim(output_name),'3D cartisian wind vector', 'm/sec', "time: point", &
1155  axes(1:3), fcst_grid, kstt_windvect,kend_windvect, dyn_bundle, output_file1, range=vrange, &
1156  l3dvector=l3dvector,rcd=rc)
1157  endif
1158 !
1159  if ( .not.hydrostatico ) then
1160  if( id_w>0 ) then
1161  call find_outputname(trim(file_name),'w',output_name)
1162  call add_field_to_bundle(trim(output_name),'vertical wind', 'm/sec', "time: point", &
1163  axes(1:3), fcst_grid, kstt_w,kend_w, dyn_bundle, output_file, &
1164  range=wrange, rcd=rc)
1165  if(rc==0) num_field_dyn=num_field_dyn+1
1166  endif
1167  if( id_pfnh>0 ) then
1168  call find_outputname(trim(file_name),'pfnh',output_name)
1169  call add_field_to_bundle(trim(output_name),'non-hydrostatic pressure', 'pa', "time: point", &
1170  axes(1:3), fcst_grid, kstt_pfnh,kend_pfnh, dyn_bundle, output_file, rcd=rc)
1171  if(rc==0) num_field_dyn=num_field_dyn+1
1172  endif
1173  if( id_delz>0 ) then
1174  call find_outputname(trim(file_name),'delz',output_name)
1175  call add_field_to_bundle(trim(output_name),'height thickness', 'm', "time: point", &
1176  axes(1:3), fcst_grid, kstt_delz,kend_delz, dyn_bundle, output_file, rcd=rc)
1177  if(rc==0) num_field_dyn=num_field_dyn+1
1178  endif
1179  else
1180  if( id_pfhy>0 ) then
1181  call find_outputname(trim(file_name),'pfhy',output_name)
1182  call add_field_to_bundle(trim(output_name),'hydrostatic pressure', 'pa', "time: point", &
1183  axes(1:3), fcst_grid, kstt_pfhy,kend_pfhy, dyn_bundle, output_file, rcd=rc)
1184  if(rc==0) num_field_dyn=num_field_dyn+1
1185  endif
1186  endif
1187 !
1188  if( id_omga>0 ) then
1189  call find_outputname(trim(file_name),'omga',output_name)
1190  call add_field_to_bundle(trim(output_name),'Vertical pressure velocity', 'pa/sec', "time: point", &
1191  axes(1:3), fcst_grid, kstt_omga,kend_omga, dyn_bundle, output_file, rcd=rc)
1192  if(rc==0) num_field_dyn=num_field_dyn+1
1193  endif
1194 !
1195  if(id_pt > 0) then
1196  call find_outputname(trim(file_name),'temp',output_name)
1197  call add_field_to_bundle(trim(output_name),'temperature', 'K', "time: point", &
1198  axes(1:3), fcst_grid, kstt_pt,kend_pt, dyn_bundle, output_file, &
1199  range=trange,rcd=rc)
1200  if(rc==0) num_field_dyn=num_field_dyn+1
1201  endif
1202 !
1203  if( id_delp > 0) then
1204  call find_outputname(trim(file_name),'delp',output_name)
1205  call add_field_to_bundle(trim(output_name),'pressure thickness', 'pa', "time: point", &
1206  axes(1:3), fcst_grid, kstt_delp,kend_delp, dyn_bundle, output_file, rcd=rc)
1207  if(rc==0) num_field_dyn=num_field_dyn+1
1208  endif
1209 !
1210 ! tracers
1211  do i=1, ncnsto
1212  call get_tracer_names ( model_atmos, i, tname, tlongname, tunits )
1213  if (id_tracer(i)>0) then
1214  call find_outputname(trim(file_name),trim(tname),output_name)
1215  call add_field_to_bundle(trim(output_name),trim(tlongname), trim(tunits), "time: point", &
1216  axes(1:3), fcst_grid, kstt_tracer(i),kend_tracer(i), dyn_bundle, output_file, rcd=rc)
1217  if(rc==0) num_field_dyn=num_field_dyn+1
1218  endif
1219 ! if(mpp_pe()==mpp_root_pe())print *,'in fv_dyn bundle,add trac,i=',i,'output_name=',trim(output_name),' rc=',rc
1220  enddo
1221 !
1222 !
1223  if( id_ps > 0) then
1224  call find_outputname(trim(file_name),'ps',output_name)
1225  call add_field_to_bundle(trim(output_name),'surface pressure', 'pa', "time: point", &
1226  axes(1:2), fcst_grid, kstt_ps,kend_ps, dyn_bundle, output_file, rcd=rc)
1227  if(rc==0) num_field_dyn=num_field_dyn+1
1228  endif
1229 !
1230  if( id_hs > 0) then
1231  call find_outputname(trim(file_name),'hs',output_name)
1232  call add_field_to_bundle(trim(output_name),'surface geopotential height', 'gpm', "time: point", &
1233  axes(1:2), fcst_grid, kstt_hs,kend_hs, dyn_bundle, output_file, rcd=rc)
1234  if(rc==0) num_field_dyn=num_field_dyn+1
1235  endif
1236 !
1237  if(id_dbz > 0) then
1238  call find_outputname(trim(file_name),'reflectivity',output_name)
1239 ! if(mpp_pe()==mpp_root_pe())print *,'reflectivity, output name=',trim(output_name)
1240  call add_field_to_bundle(trim(output_name),'Stoelinga simulated reflectivity', 'dBz', "time: point", &
1241  axes(1:3), fcst_grid, kstt_dbz,kend_dbz, dyn_bundle, output_file, rcd=rc)
1242  if(rc==0) num_field_dyn=num_field_dyn+1
1243  endif
1244  if(id_ustm > 0 .and. id_vstm > 0 .and. id_srh01 > 0 .and. id_srh03 > 0) then
1245  call find_outputname(trim(file_name),'ustm',output_name)
1246  if(mpp_pe()==mpp_root_pe())print *,'u comp. of storm motion, output name=',trim(output_name)
1247  call add_field_to_bundle(trim(output_name),'u comp of storm motion', 'm/s', "time: point", &
1248  axes(1:2), fcst_grid, kstt_ustm,kend_ustm, dyn_bundle, output_file, rcd=rc)
1249  if(rc==0) num_field_dyn=num_field_dyn+1
1250 
1251  call find_outputname(trim(file_name),'vstm',output_name)
1252  if(mpp_pe()==mpp_root_pe())print *,'v comp. of storm motion, output name=',trim(output_name)
1253  call add_field_to_bundle(trim(output_name),'v comp of storm motion', 'm/s', "time: point", &
1254  axes(1:2), fcst_grid, kstt_vstm,kend_vstm, dyn_bundle, output_file, rcd=rc)
1255  if(rc==0) num_field_dyn=num_field_dyn+1
1256 
1257  call find_outputname(trim(file_name),'srh01',output_name)
1258  if(mpp_pe()==mpp_root_pe())print *,'0-1km srh, output name=',trim(output_name)
1259  call add_field_to_bundle(trim(output_name),'0-1km srh', 'm/s**2', "time: point", &
1260  axes(1:2), fcst_grid, kstt_srh01,kend_srh01, dyn_bundle, output_file, rcd=rc)
1261  if(rc==0) num_field_dyn=num_field_dyn+1
1262 
1263  call find_outputname(trim(file_name),'srh03',output_name)
1264  if(mpp_pe()==mpp_root_pe())print *,'0-3km srh, output name=',trim(output_name)
1265  call add_field_to_bundle(trim(output_name),'0-3km srh', 'm/s**2', "time: point", &
1266  axes(1:2), fcst_grid, kstt_srh03,kend_srh03, dyn_bundle, output_file, rcd=rc)
1267  if(rc==0) num_field_dyn=num_field_dyn+1
1268  endif
1269 
1270 
1271  if(id_maxvort01 > 0) then
1272  call find_outputname(trim(file_name),'maxvort01',output_name)
1273  if(mpp_pe()==mpp_root_pe())print *,'max hourly 0-1km vert. vorticity, output name=',trim(output_name)
1274  call add_field_to_bundle(trim(output_name),'Max hourly 0-1km vert. vorticity', '1/s', "time: point", &
1275  axes(1:2), fcst_grid, kstt_maxvort01,kend_maxvort01, dyn_bundle, output_file, rcd=rc)
1276  if(rc==0) num_field_dyn=num_field_dyn+1
1277  endif
1278  if(id_maxvort02 > 0) then
1279  call find_outputname(trim(file_name),'maxvort02',output_name)
1280  if(mpp_pe()==mpp_root_pe())print *,'max hourly 0-2km vert. vorticity, output name=',trim(output_name)
1281  call add_field_to_bundle(trim(output_name),'Max hourly 0-2km vert. vorticity', '1/s', "time: point", &
1282  axes(1:2), fcst_grid, kstt_maxvort02,kend_maxvort02, dyn_bundle, output_file, rcd=rc)
1283  if(rc==0) num_field_dyn=num_field_dyn+1
1284  endif
1285  if(id_maxvorthy1 > 0) then
1286  call find_outputname(trim(file_name),'maxvorthy1',output_name)
1287  if(mpp_pe()==mpp_root_pe())print *,'max hourly lev 1 vert. vorticity output name=',trim(output_name)
1288  call add_field_to_bundle(trim(output_name),'Max hourly lev 1 vert vort.', '1/s', "time: point", &
1289  axes(1:2), fcst_grid, kstt_maxvorthy1,kend_maxvorthy1, dyn_bundle, output_file, rcd=rc)
1290  if(rc==0) num_field_dyn=num_field_dyn+1
1291  endif
1292  if(id_wmaxup > 0) then
1293  call find_outputname(trim(file_name),'wmaxup',output_name)
1294  if(mpp_pe()==mpp_root_pe())print *,'max hourly updraft vel, output name=',trim(output_name)
1295  call add_field_to_bundle(trim(output_name),'Max hourly updraft velocity', 'm/s', "time: point", &
1296  axes(1:2), fcst_grid, kstt_wup,kend_wup, dyn_bundle, output_file, rcd=rc)
1297  if(rc==0) num_field_dyn=num_field_dyn+1
1298  endif
1299  if(id_wmaxdn > 0) then
1300  call find_outputname(trim(file_name),'wmaxdn',output_name)
1301  if(mpp_pe()==mpp_root_pe())print *,'max hourly downdraft vel, output name=',trim(output_name)
1302  call add_field_to_bundle(trim(output_name),'Max hourly downdraft velocity', 'm/s', "time: point", &
1303  axes(1:2), fcst_grid, kstt_wdn,kend_wdn, dyn_bundle, output_file, rcd=rc)
1304  if(rc==0) num_field_dyn=num_field_dyn+1
1305  endif
1306  if( .not.hydrostatico .and. id_uhmax03 > 0 ) then
1307  call find_outputname(trim(file_name),'uhmax03',output_name)
1308  if(mpp_pe()==mpp_root_pe())print *,'max hourly 0-3km updraft helicity, output name=',trim(output_name)
1309  call add_field_to_bundle(trim(output_name),'Max hourly 0-3km updraft helicity', 'm/s**2', "time: point", &
1310  axes(1:2), fcst_grid, kstt_uhmax03,kend_uhmax03, dyn_bundle, output_file, rcd=rc)
1311  if(rc==0) num_field_dyn=num_field_dyn+1
1312  endif
1313  if( .not.hydrostatico .and. id_uhmin03 > 0 ) then
1314  call find_outputname(trim(file_name),'uhmin03',output_name)
1315  if(mpp_pe()==mpp_root_pe())print *,'max hourly 0-3km updraft helicity, output name=',trim(output_name)
1316  call add_field_to_bundle(trim(output_name),'Max hourly 0-3km updraft helicity', 'm/s**2', "time: point", &
1317  axes(1:2), fcst_grid, kstt_uhmin03,kend_uhmin03, dyn_bundle, output_file, rcd=rc)
1318  if(rc==0) num_field_dyn=num_field_dyn+1
1319  endif
1320  if( .not.hydrostatico .and. id_uhmax25 > 0 ) then
1321  call find_outputname(trim(file_name),'uhmax25',output_name)
1322  if(mpp_pe()==mpp_root_pe())print *,'max hourly 2-5km updraft helicity, output name=',trim(output_name)
1323  call add_field_to_bundle(trim(output_name),'Max hourly 2-5km updraft helicity', 'm/s**2', "time: point", &
1324  axes(1:2), fcst_grid, kstt_uhmax25,kend_uhmax25, dyn_bundle, output_file, rcd=rc)
1325  if(rc==0) num_field_dyn=num_field_dyn+1
1326  endif
1327  if( .not.hydrostatico .and. id_uhmin25 > 0 ) then
1328  call find_outputname(trim(file_name),'uhmin25',output_name)
1329  if(mpp_pe()==mpp_root_pe())print *,'max hourly 2-5km updraft helicity, output name=',trim(output_name)
1330  call add_field_to_bundle(trim(output_name),'Max hourly 2-5km updraft helicity', 'm/s**2', "time: point", &
1331  axes(1:2), fcst_grid, kstt_uhmin25,kend_uhmin25, dyn_bundle, output_file, rcd=rc)
1332  if(rc==0) num_field_dyn=num_field_dyn+1
1333  endif
1334 
1335 !jwtest:
1336 ! call ESMF_FieldBundleGet(dyn_bundle, fieldCount=fieldCount, rc=rc)
1337 ! print *,'in dyn_bundle_setup, fieldCount=',fieldCount
1338 ! allocate(fieldnamelist(fieldCount),fieldlist(fieldCount))
1339 ! call ESMF_FieldBundleGet(dyn_bundle, fieldlist=fieldlist,fieldnamelist=fieldnamelist, rc=rc)
1340 ! do i=1,fieldCount
1341 ! call ESMF_AttributeGet(fieldlist(i), convention="NetCDF", purpose="FV3", &
1342 ! name="output_file", value=fld_outfilename, rc=rc)
1343 ! print *,'in dyn bundle setup, i=',i,' fieldname=',trim(fieldnamelist(i)),' out filename=',trim(fld_outfilename)
1344 ! enddo
1345 
1346  end subroutine fv_dyn_bundle_setup
1347 
1348  subroutine add_field_to_bundle(var_name,long_name,units,cell_methods, axes,dyn_grid, &
1349  kstt,kend,dyn_bundle,output_file, range, l3Dvector, rcd)
1350  use esmf
1351  implicit none
1352 
1353  character(*), intent(in) :: var_name, long_name, units, cell_methods
1354  integer, intent(in) :: axes(:)
1355  type(esmf_grid), intent(in) :: dyn_grid
1356  integer, intent(in) :: kstt,kend
1357  type(esmf_fieldbundle),intent(inout) :: dyn_bundle
1358  character(*), intent(in) :: output_file
1359  real, intent(in), optional :: range(2)
1360  logical, intent(in), optional :: l3Dvector
1361  integer, intent(out), optional :: rcd
1362 !
1363 !*** local variable
1364  type(esmf_field) :: field
1365  type(esmf_datacopy_flag) :: copyflag=esmf_datacopy_reference
1366  integer rc, i, j, idx
1367  real(4),dimension(:,:,:,:),pointer :: temp_r4d
1368  real(4),dimension(:,:,:), pointer :: temp_r3d
1369  real(4),dimension(:,:), pointer :: temp_r2d
1370  logical, save :: first=.true.
1371 !
1372 !*** create esmf field
1373  if( present(l3dvector) ) then
1374  temp_r4d => windvect(1:3,isco:ieco,jsco:jeco,kstt:kend)
1375  call esmf_logwrite('create winde vector esmf field', esmf_logmsg_info, rc=rc)
1376  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1377  line=__line__, &
1378  file=__file__)) &
1379  return ! bail out
1380 !jw field = ESMF_FieldCreate(dyn_grid, temp_r4d, datacopyflag=ESMF_DATACOPY_VALUE,
1381  field = esmf_fieldcreate(dyn_grid, temp_r4d, datacopyflag=esmf_datacopy_reference, &
1382  gridtofieldmap=(/2,3/), ungriddedlbound=(/1,kstt/), ungriddedubound=(/3,kend/), &
1383  name="windvector", indexflag=esmf_index_delocal, rc=rc)
1384  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1385  line=__line__, &
1386  file=__file__)) &
1387  return ! bail out
1388  call esmf_logwrite('create winde vector esmf field', esmf_logmsg_info, rc=rc)
1389  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1390  line=__line__, &
1391  file=__file__)) &
1392  return ! bail out
1393 
1394  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1395  attrlist=(/"output_file"/), rc=rc)
1396  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1397  name='output_file',value=trim(output_file),rc=rc)
1398 
1399  call esmf_fieldbundleadd(dyn_bundle,(/field/), rc=rc)
1400  if( present(rcd)) rcd=rc
1401  return
1402  else if( kend>kstt ) then
1403  temp_r3d => buffer_dyn(isco:ieco,jsco:jeco,kstt:kend)
1404  field = esmf_fieldcreate(dyn_grid, temp_r3d, datacopyflag=copyflag, &
1405  name=var_name, indexflag=esmf_index_delocal, rc=rc)
1406  else if(kend==kstt) then
1407  temp_r2d => buffer_dyn(isco:ieco,jsco:jeco,kstt)
1408  field = esmf_fieldcreate(dyn_grid, temp_r2d, datacopyflag=copyflag, &
1409  name=var_name, indexflag=esmf_index_delocal, rc=rc)
1410  endif
1411 !
1412 !*** add field attributes
1413  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1414  attrlist=(/"long_name"/), rc=rc)
1415  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1416  name='long_name',value=trim(long_name),rc=rc)
1417 
1418  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1419  attrlist=(/"units"/), rc=rc)
1420  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1421  name='units',value=trim(units),rc=rc)
1422 
1423  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1424  attrlist=(/"missing_value"/), rc=rc)
1425  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1426  name='missing_value',value=missing_value,rc=rc)
1427 
1428  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1429  attrlist=(/"_FillValue"/), rc=rc)
1430  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1431  name='_FillValue',value=missing_value,rc=rc)
1432 
1433  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1434  attrlist=(/"cell_methods"/), rc=rc)
1435  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1436  name='cell_methods',value=trim(cell_methods),rc=rc)
1437 
1438  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1439  attrlist=(/"output_file"/), rc=rc)
1440  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1441  name='output_file',value=trim(output_file),rc=rc)
1442 !
1443 !*** add vertical coord attribute:
1444  if( size(axes) > 2) then
1445  do i=3,size(axes)
1446  idx=0
1447  do j=1,size(all_axes)
1448  if (axes(i)==all_axes(j)) then
1449  idx=j
1450  exit
1451  endif
1452  enddo
1453  if (idx>0) then
1454  call esmf_attributeadd(field, convention="NetCDF", purpose="FV3", &
1455  attrlist=(/"ESMF:ungridded_dim_labels"/), rc=rc)
1456  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1457  line=__line__, &
1458  file=__file__)) &
1459  return ! bail out
1460  call esmf_attributeset(field, convention="NetCDF", purpose="FV3", &
1461  name="ESMF:ungridded_dim_labels", valuelist=(/trim(axis_name(idx))/), rc=rc)
1462 ! if( first ) then
1463 ! print *,'add axis_name to field,',trim(axis_name(idx))
1464 ! endif
1465  if (esmf_logfounderror(rctocheck=rc, msg=esmf_logerr_passthru, &
1466  line=__line__, &
1467  file=__file__)) &
1468  return ! bail out
1469  endif
1470  enddo
1471  first=.false.
1472  endif
1473 
1474 !*** add field into bundle
1475  call esmf_fieldbundleadd(dyn_bundle,(/field/), rc=rc)
1476  if( present(rcd)) rcd=rc
1477 !
1478  end subroutine add_field_to_bundle
1479 !-------------------------------------------------------------------------------------
1480  subroutine find_outputname(module_name, field_name, output_name)
1481  character(*), intent(in) :: module_name
1482  character(*), intent(in) :: field_name
1483  character(*), intent(out) :: output_name
1484 !
1485  integer i,j,in_num, out_num
1486  integer tile_count
1487 !
1488  tile_count=1
1489  in_num = find_input_field(module_name, field_name, tile_count)
1490 !
1491  output_name=''
1492  do i=1, max_output_fields
1493  if(output_fields(i)%input_field == in_num) then
1494  output_name=output_fields(i)%output_name
1495  exit
1496  endif
1497  enddo
1498  if(output_name=='') then
1499  print *,'Error, cant find out put name, field_name=',trim(field_name),'in_num=',in_num
1500  endif
1501 
1502  end subroutine find_outputname
1503 
1504 #endif
1505 
1506 end module fv_nggps_diags_mod
character(len=256) tlongname
real, dimension(:,:), allocatable maxvorthy1
subroutine store_data(id, work, Time, nstt, nend)
real, dimension(2) vrange
winds
subroutine, public max_uh(is, ie, js, je, ng, km, zvir, sphum, uphmax, uphmin, w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
real, dimension(:), allocatable bk
real, dimension(:,:), allocatable maxvort02
real, dimension(:,:), allocatable lat
subroutine, public fv_nggps_tavg(Atm, Time_step_atmos, avg_max_length, zvir)
The module &#39;multi_gases&#39; peforms multi constitutents computations.
Definition: multi_gases.F90:25
logical use_wrtgridcomp_output
real, dimension(:,:), allocatable uhmin25
real, parameter stndrd_atmos_lapse
real, dimension(:), allocatable ak
real, dimension(:,:), allocatable dn2
subroutine, public bunkers_vector(is, ie, js, je, ng, km, zvir, sphum, uc, vc, ua, va, delz, q, hydrostatic, pt, peln, phis, grav)
character(len=128) tname
integer ice_wat
GFDL physics.
integer, dimension(:), allocatable all_axes
subroutine, public max_vorticity_hy1(is, ie, js, je, km, vort, maxvorthy1)
integer, dimension(:), allocatable id_tracer
pure real function, public virq(q)
integer, dimension(:), allocatable kend_tracer
The module &#39;fv_nggps_diags&#39; computes output diagnostics entirely on 3D pressure levels.
real, dimension(:,:), allocatable lon
real(4), dimension(:,:,:,:), allocatable, target windvect
subroutine, public fv_nggps_diag(Atm, zvir, Time)
character(len=256) tunits
real, dimension(:,:), allocatable uhmax03
real, dimension(:,:), allocatable maxvort01
subroutine, public fv_nggps_diag_init(Atm, axes, Time)
The module &#39;fv_arrays&#39; contains the &#39;fv_atmos_type&#39; and associated datatypes.
Definition: fv_arrays.F90:24
real(4), dimension(:,:), allocatable, target psurf
real, dimension(:,:), allocatable uhmin03
character(20), dimension(:), allocatable axis_name
real, parameter stndrd_atmos_ps
real, dimension(:,:), allocatable up2
real, dimension(2) trange
temperature
integer, dimension(:), allocatable kstt_tracer
subroutine, public dbzcalc(q, pt, delp, peln, delz, dbz, maxdbz, allmax, bd, npz, ncnst, hydrostatic, zvir, in0r, in0s, in0g, iliqskin)
The subroutine &#39;dbzcalc&#39; computes equivalent reflectivity factor (in dBZ) at each model grid point...
character(len=64) file_name
@ The module &#39;fv_diagnostics&#39; contains routines to compute diagnosic fields.
real, dimension(2) skrange
dissipation estimate for SKEB
subroutine, public helicity_relative_caps(is, ie, js, je, ng, km, zvir, sphum, srh, uc, vc, ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
subroutine, public max_vorticity(is, ie, js, je, ng, km, zvir, sphum, delz, q, hydrostatic, pt, peln, phis, grav, vort, maxvort, z_bot, z_top)
real, dimension(2) wrange
vertical wind
real(4), dimension(:,:,:), allocatable, target buffer_dyn
real, dimension(:,:), allocatable uhmax25
character(20), dimension(:), allocatable axis_name_vert
subroutine, public max_vv(is, ie, js, je, npz, ng, up2, dn2, pe, w)
real, parameter missing_value
subroutine, public get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)