120 use constants_mod
, only: grav, rdgas, rvgas, pi=>pi_8, radius, kappa, wtmair, wtmco2, &
121 omega, hlv, cp_air, cp_vapor, tfreeze
122 use fms_mod
, only: write_version_number
123 use fms_io_mod
, only: set_domain, nullify_domain, write_version_number
124 use time_manager_mod
, only: time_type, get_date, get_time
125 use mpp_domains_mod
, only: domain2d, mpp_update_domains, dgrid_ne, east, north
126 use diag_manager_mod
, only: diag_axis_init, register_diag_field, &
127 register_static_field, send_data, diag_grid_init
131 use fv_mp_mod, only: mp_reduce_sum, mp_reduce_min, mp_reduce_max, is_master
138 use tracer_manager_mod
, only: get_tracer_names, get_number_tracers, get_tracer_index
139 use field_manager_mod
, only: model_atmos
140 use mpp_mod
, only: mpp_error, fatal, stdlog, mpp_pe, mpp_root_pe, mpp_sum, mpp_max, note, input_nml_file
141 use mpp_io_mod
, only: mpp_flush
142 use sat_vapor_pres_mod
, only: compute_qs, lookup_es
146 use gfdl_cloud_microphys_mod
, only: wqs1, qsmith_init
149 use column_diagnostics_mod
, only: column_diagnostics_init, &
150 initialize_diagnostic_columns, &
151 column_diagnostics_header, &
152 close_column_diagnostics_units
173 character(len=3) ::
gn =
'' 203 integer,
parameter ::
nplev = 10
242 #include<file_version.h> 246 subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
248 integer,
intent(out) :: axes(4)
249 type(time_type),
intent(in) :: Time
250 integer,
intent(in) :: npx, npy, npz
251 real,
intent(in):: p_ref
253 real,
allocatable :: grid_xt(:), grid_yt(:), grid_xe(:), grid_ye(:), grid_xn(:), grid_yn(:)
254 real,
allocatable :: grid_x(:), grid_y(:)
255 real,
allocatable :: a3(:,:,:)
257 real :: hyam(npz), hybm(npz)
260 integer :: id_bk, id_pk, id_area, id_lon, id_lat, id_lont, id_latt, id_phalf, id_pfull
261 integer :: id_hyam, id_hybm
263 integer :: i, j, k, m, n, ntileMe, id_xt, id_yt, id_x, id_y, id_xe, id_ye, id_xn, id_yn
264 integer :: isc, iec, jsc, jec
268 character(len=64) :: plev
269 character(len=64) :: field
276 character(len=64) :: errmsg
278 integer :: nlunit, ios
280 call write_version_number (
'FV_DIAGNOSTICS_MOD', version )
281 idiag => atm(1)%idiag
291 call set_domain(atm(1)%domain)
293 sphum = get_tracer_index(model_atmos,
'sphum')
294 liq_wat = get_tracer_index(model_atmos,
'liq_wat')
295 ice_wat = get_tracer_index(model_atmos,
'ice_wat')
297 rainwat = get_tracer_index(model_atmos,
'rainwat')
298 snowwat = get_tracer_index(model_atmos,
'snowwat')
299 graupel = get_tracer_index(model_atmos,
'graupel')
300 o3mr = get_tracer_index(model_atmos,
'o3mr')
301 cld_amt = get_tracer_index(model_atmos,
'cld_amt')
309 vrange = (/ -330., 330. /)
310 wrange = (/ -100., 100. /)
318 skrange = (/ -10000000.0, 10000000.0 /)
321 if (atm(1)%grid_number == 1)
fv_time = time
323 allocate (
idiag%phalf(npz+1) )
328 if ( pfull(k) > 30.e2 )
then 333 if ( is_master() )
write(*,*)
'mp_top=',
mp_top,
'pfull=', pfull(
mp_top)
336 allocate(grid_xt(npx-1), grid_yt(npy-1))
337 grid_xt = (/ (i, i=1,npx-1) /)
338 grid_yt = (/ (j, j=1,npy-1) /)
344 allocate(grid_x(npx), grid_y(npy))
345 grid_x = (/ (i, i=1,npx) /)
346 grid_y = (/ (j, j=1,npy) /)
349 isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
350 jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
353 call diag_grid_init(domain=atm(n)%domain, &
354 & glo_lon=
rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,1), &
355 & glo_lat=
rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,2), &
356 & aglo_lon=
rad2deg*atm(n)%gridstruct%agrid(isc-1:iec+1,jsc-1:jec+1,1), &
357 & aglo_lat=
rad2deg*atm(n)%gridstruct%agrid(isc-1:iec+1,jsc-1:jec+1,2))
359 ntileme =
size(atm(:))
360 if (ntileme > 1)
call mpp_error(fatal,
"fv_diag_init can only be called with one grid at a time.")
366 id_xt = diag_axis_init(
'grid_xt',grid_xt,
'degrees_E',
'x',
'T-cell longitude', &
367 set_name=trim(field),domain2=atm(n)%Domain, tile_count=n)
368 id_yt = diag_axis_init(
'grid_yt',grid_yt,
'degrees_N',
'y',
'T-cell latitude', &
369 set_name=trim(field), domain2=atm(n)%Domain, tile_count=n)
380 id_x = diag_axis_init(
'grid_x',grid_x,
'degrees_E',
'x',
'cell corner longitude', &
381 set_name=trim(field),domain2=atm(n)%Domain, tile_count=n, domain_position=east)
382 id_y = diag_axis_init(
'grid_y',grid_y,
'degrees_N',
'y',
'cell corner latitude', &
383 set_name=trim(field), domain2=atm(n)%Domain, tile_count=n, domain_position=north)
387 deallocate(grid_xt, grid_yt)
388 deallocate(grid_x, grid_y )
390 id_phalf = diag_axis_init(
'phalf',
idiag%phalf,
'mb',
'z', &
391 'ref half pressure level', direction=-1, set_name=
"dynamics")
392 id_pfull = diag_axis_init(
'pfull', pfull,
'mb',
'z', &
393 'ref full pressure level', direction=-1, set_name=
"dynamics", edges=id_phalf)
397 id_bk = register_static_field(
"dynamics",
'bk', (/id_phalf/), &
398 'vertical coordinate sigma value',
'none' )
400 id_pk = register_static_field(
"dynamics",
'pk', (/id_phalf/), &
401 'pressure part of the hybrid coordinate',
'pascal' )
403 id_hyam = register_static_field(
"dynamics",
'hyam', (/id_pfull/), &
404 'vertical coordinate A value',
'1E-5 Pa' )
406 id_hybm = register_static_field(
"dynamics",
'hybm', (/id_pfull/), &
407 'vertical coordinate B value',
'none' )
411 if ( id_bk > 0 ) used = send_data( id_bk,atm(1)%bk, time )
412 if ( id_pk > 0 ) used = send_data( id_pk,atm(1)%ak, time )
413 if ( id_hyam > 0 )
then 415 hyam(k) = 0.5 * ( atm(1)%ak(k) + atm(1)%ak(k+1) ) * 1.e-5
417 used = send_data( id_hyam, hyam, time )
419 if ( id_hybm > 0 )
then 421 hybm(k) = 0.5 * ( atm(1)%bk(k) + atm(1)%bk(k+1) )
423 used = send_data( id_hybm, hybm, time )
436 levs = (/50,100,200,250,300,500,750,850,925,1000/)
441 levs = (/1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,925,950,975,1000/)
448 id_plev = diag_axis_init(
'plev',
levs(:)*1.0,
'mb',
'z', &
449 'actual pressure level', direction=-1, set_name=
"dynamics")
460 id_lon = register_static_field( trim(field),
'grid_lon', (/id_x,id_y/), &
461 'longitude',
'degrees_E' )
462 id_lat = register_static_field( trim(field),
'grid_lat', (/id_x,id_y/), &
463 'latitude',
'degrees_N' )
464 id_lont = register_static_field( trim(field),
'grid_lont', (/id_xt,id_yt/), &
465 'longitude',
'degrees_E' )
466 id_latt = register_static_field( trim(field),
'grid_latt', (/id_xt,id_yt/), &
467 'latitude',
'degrees_N' )
468 id_area = register_static_field( trim(field),
'area', axes(1:2), &
469 'cell area',
'm**2' )
471 idiag%id_zsurf = register_static_field( trim(field),
'zsurf', axes(1:2), &
472 'surface height',
'm' )
474 idiag%id_zs = register_static_field( trim(field),
'zs', axes(1:2), &
475 'Original Mean Terrain',
'm' )
477 idiag%id_ze = register_static_field( trim(field),
'ze', axes(1:3), &
478 'Hybrid_Z_surface',
'm' )
479 idiag%id_oro = register_static_field( trim(field),
'oro', axes(1:2), &
480 'Land/Water Mask',
'none' )
481 idiag%id_sgh = register_static_field( trim(field),
'sgh', axes(1:2), &
482 'Terrain Standard deviation',
'm' )
489 idiag%ic_ps = register_static_field( trim(field),
'ps_ic', axes(1:2), &
490 'initial surface pressure',
'Pa' )
491 idiag%ic_ua = register_static_field( trim(field),
'ua_ic', axes(1:3), &
492 'zonal wind',
'm/sec' )
493 idiag%ic_va = register_static_field( trim(field),
'va_ic', axes(1:3), &
494 'meridional wind',
'm/sec' )
495 idiag%ic_ppt= register_static_field( trim(field),
'ppt_ic', axes(1:3), &
496 'potential temperature perturbation',
'K' )
497 idiag%ic_sphum = register_static_field( trim(field),
'sphum_ic', axes(1:2), &
498 'initial surface pressure',
'Pa' )
502 master = (mpp_pe()==mpp_root_pe())
505 isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
506 jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
508 allocate (
idiag%zsurf(isc:iec,jsc:jec) )
512 idiag%zsurf(i,j) =
ginv * atm(n)%phis(i,j)
520 isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
521 jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
522 if (id_lon > 0) used = send_data(id_lon,
rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,1), time)
523 if (id_lat > 0) used = send_data(id_lat,
rad2deg*atm(n)%gridstruct%grid(isc:iec+1,jsc:jec+1,2), time)
524 if (id_lont > 0) used = send_data(id_lont,
rad2deg*atm(n)%gridstruct%agrid(isc:iec,jsc:jec,1), time)
525 if (id_latt > 0) used = send_data(id_latt,
rad2deg*atm(n)%gridstruct%agrid(isc:iec,jsc:jec,2), time)
526 if (id_area > 0) used = send_data(id_area, atm(n)%gridstruct%area(isc:iec,jsc:jec), time)
528 if (
idiag%id_zsurf > 0) used = send_data(
idiag%id_zsurf,
idiag%zsurf, time)
530 if ( atm(n)%flagstruct%fv_land )
then 532 if (
idiag%id_oro > 0) used = send_data(
idiag%id_oro, atm(n)%oro(isc:iec,jsc:jec), time)
533 if (
idiag%id_sgh > 0) used = send_data(
idiag%id_sgh, atm(n)%sgh(isc:iec,jsc:jec), time)
536 if ( atm(n)%flagstruct%ncep_ic )
then 537 if (
idiag%id_ts > 0) used = send_data(
idiag%id_ts, atm(n)%ts(isc:iec,jsc:jec), time)
540 if ( atm(n)%flagstruct%hybrid_z .and.
idiag%id_ze > 0 ) &
541 used = send_data(
idiag%id_ze, atm(n)%ze0(isc:iec,jsc:jec,1:npz), time)
543 if (
idiag%ic_ps > 0) used = send_data(
idiag%ic_ps, atm(n)%ps(isc:iec,jsc:jec)*
ginv, time)
545 if(
idiag%ic_ua > 0) used=send_data(
idiag%ic_ua, atm(n)%ua(isc:iec,jsc:jec,:), time)
546 if(
idiag%ic_va > 0) used=send_data(
idiag%ic_va, atm(n)%va(isc:iec,jsc:jec,:), time)
548 pk0 = 1000.e2 ** kappa
549 if(
idiag%ic_ppt> 0)
then 551 allocate (
idiag%pt1(npz) )
552 allocate ( a3(isc:iec,jsc:jec,npz) )
554 call gw_1d(npz, 1000.e2, atm(n)%ak, atm(n)%ak, atm(n)%ak(1), 10.e3,
idiag%pt1)
561 a3(i,j,k) = (atm(n)%pt(i,j,k)/atm(n)%pkz(i,j,k) -
idiag%pt1(k)) *
pk0 565 used=send_data(
idiag%ic_ppt, a3, time)
567 deallocate (
idiag%pt1 )
575 allocate(
idiag%id_tracer(ncnst))
576 allocate(
idiag%id_tracer_dmmr(ncnst))
577 allocate(
idiag%id_tracer_dvmr(ncnst))
578 allocate(
idiag%w_mr(ncnst))
579 idiag%id_tracer(:) = 0
580 idiag%id_tracer_dmmr(:) = 0
581 idiag%id_tracer_dvmr(:) = 0
602 idiag%id_zsurf = register_diag_field( trim(field),
'zsurf', axes(1:2), time, &
603 'surface height',
'm')
608 idiag%id_ps = register_diag_field( trim(field),
'ps', axes(1:2), time, &
614 idiag%id_mq = register_diag_field( trim(field),
'mq', axes(1:2), time, &
619 idiag%id_aam = register_diag_field( trim(field),
'aam', axes(1:2), time, &
621 idiag%id_amdt = register_diag_field( trim(field),
'amdt', axes(1:2), time, &
627 if (atm(n)%flagstruct%write_3d_diags)
then 629 idiag%id_T_dt_phys = register_diag_field( trim(field),
'T_dt_phys', axes(1:3), time, &
631 if (
idiag%id_T_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_t_dt(isc:iec,jsc:jec,npz))
632 idiag%id_u_dt_phys = register_diag_field( trim(field),
'u_dt_phys', axes(1:3), time, &
634 if (
idiag%id_u_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_u_dt(isc:iec,jsc:jec,npz))
635 idiag%id_v_dt_phys = register_diag_field( trim(field),
'v_dt_phys', axes(1:3), time, &
637 if (
idiag%id_v_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_v_dt(isc:iec,jsc:jec,npz))
639 idiag%id_qv_dt_phys = register_diag_field( trim(field),
'qv_dt_phys', axes(1:3), time, &
641 if (
idiag%id_qv_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_qv_dt(isc:iec,jsc:jec,npz))
642 idiag%id_ql_dt_phys = register_diag_field( trim(field),
'ql_dt_phys', axes(1:3), time, &
644 if (
idiag%id_ql_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_ql_dt(isc:iec,jsc:jec,npz))
645 idiag%id_qi_dt_phys = register_diag_field( trim(field),
'qi_dt_phys', axes(1:3), time, &
647 if (
idiag%id_qi_dt_phys > 0)
allocate (atm(n)%phys_diag%phys_qi_dt(isc:iec,jsc:jec,npz))
652 write(plev,
'(I5)')
levs(i)
654 idiag%id_h(i) = register_diag_field(trim(field),
'z'//trim(adjustl(plev)), axes(1:2), time, &
657 idiag%id_u(i) = register_diag_field(trim(field),
'u'//trim(adjustl(plev)), axes(1:2), time, &
660 idiag%id_v(i) = register_diag_field(trim(field),
'v'//trim(adjustl(plev)), axes(1:2), time, &
663 idiag%id_t(i) = register_diag_field(trim(field),
't'//trim(adjustl(plev)), axes(1:2), time, &
666 idiag%id_q(i) = register_diag_field(trim(field),
'q'//trim(adjustl(plev)), axes(1:2), time, &
669 idiag%id_omg(i) = register_diag_field(trim(field),
'omg'//trim(adjustl(plev)), axes(1:2), time, &
673 if (atm(n)%flagstruct%write_3d_diags)
then 674 idiag%id_u_plev = register_diag_field( trim(field),
'u_plev', axe2(1:3), time, &
676 idiag%id_v_plev = register_diag_field( trim(field),
'v_plev', axe2(1:3), time, &
678 idiag%id_t_plev = register_diag_field( trim(field),
't_plev', axe2(1:3), time, &
680 idiag%id_h_plev = register_diag_field( trim(field),
'h_plev', axe2(1:3), time, &
682 idiag%id_q_plev = register_diag_field( trim(field),
'q_plev', axe2(1:3), time, &
684 idiag%id_omg_plev = register_diag_field( trim(field),
'omg_plev', axe2(1:3), time, &
690 if ( all(
idiag%id_h(minloc(abs(
levs-10)))>0) .or. all(
idiag%id_h(minloc(abs(
levs-50)))>0) .or. &
691 all(
idiag%id_h(minloc(abs(
levs-100)))>0) .or. all(
idiag%id_h(minloc(abs(
levs-200)))>0) .or. &
692 all(
idiag%id_h(minloc(abs(
levs-250)))>0) .or. all(
idiag%id_h(minloc(abs(
levs-300)))>0) .or. &
693 all(
idiag%id_h(minloc(abs(
levs-500)))>0) .or. all(
idiag%id_h(minloc(abs(
levs-700)))>0) .or. &
694 all(
idiag%id_h(minloc(abs(
levs-850)))>0) .or. all(
idiag%id_h(minloc(abs(
levs-1000)))>0) )
then 695 idiag%id_any_hght = 1
697 idiag%id_any_hght = 0
702 idiag%id_tm = register_diag_field(trim(field),
'tm', axes(1:2), time, &
708 idiag%id_slp = register_diag_field(trim(field),
'slp', axes(1:2), time, &
714 idiag%id_pmask = register_diag_field(trim(field),
'pmask', axes(1:2), time, &
715 'masking pressure at lowest level',
'mb', &
720 idiag%id_pmaskv2 = register_diag_field(trim(field),
'pmaskv2', axes(1:2), time,&
727 idiag%id_c15 = register_diag_field(trim(field),
'cat15', axes(1:2), time, &
729 idiag%id_c25 = register_diag_field(trim(field),
'cat25', axes(1:2), time, &
731 idiag%id_c35 = register_diag_field(trim(field),
'cat35', axes(1:2), time, &
733 idiag%id_c45 = register_diag_field(trim(field),
'cat45', axes(1:2), time, &
736 idiag%id_f15 = register_diag_field(trim(field),
'f15', axes(1:2), time, &
738 idiag%id_f25 = register_diag_field(trim(field),
'f25', axes(1:2), time, &
740 idiag%id_f35 = register_diag_field(trim(field),
'f35', axes(1:2), time, &
742 idiag%id_f45 = register_diag_field(trim(field),
'f45', axes(1:2), time, &
747 if (atm(n)%flagstruct%write_3d_diags)
then 748 idiag%id_ua = register_diag_field( trim(field),
'ucomp', axes(1:3), time, &
750 idiag%id_va = register_diag_field( trim(field),
'vcomp', axes(1:3), time, &
752 if ( .not. atm(n)%flagstruct%hydrostatic ) &
753 idiag%id_w = register_diag_field( trim(field),
'w', axes(1:3), time, &
756 idiag%id_pt = register_diag_field( trim(field),
'temp', axes(1:3), time, &
758 idiag%id_ppt = register_diag_field( trim(field),
'ppt', axes(1:3), time, &
760 idiag%id_theta_e = register_diag_field( trim(field),
'theta_e', axes(1:3), time, &
762 idiag%id_omga = register_diag_field( trim(field),
'omega', axes(1:3), time, &
764 idiag%id_divg = register_diag_field( trim(field),
'divg', axes(1:3), time, &
767 idiag%id_hght3d = register_diag_field( trim(field),
'hght', axes(1:3), time, &
771 idiag%id_diss = register_diag_field( trim(field),
'diss_est', axes(1:3), time, &
774 idiag%id_rh = register_diag_field( trim(field),
'rh', axes(1:3), time, &
777 idiag%id_delp = register_diag_field( trim(field),
'delp', axes(1:3), time, &
779 if ( .not. atm(n)%flagstruct%hydrostatic ) &
780 idiag%id_delz = register_diag_field( trim(field),
'delz', axes(1:3), time, &
782 if( atm(n)%flagstruct%hydrostatic )
then 783 idiag%id_pfhy = register_diag_field( trim(field),
'pfhy', axes(1:3), time, &
786 idiag%id_pfnh = register_diag_field( trim(field),
'pfnh', axes(1:3), time, &
789 idiag%id_zratio = register_diag_field( trim(field),
'zratio', axes(1:3), time, &
794 idiag%id_qn = register_diag_field( trim(field),
'qn', axes(1:3), time, &
796 idiag%id_qp = register_diag_field( trim(field),
'qp', axes(1:3), time, &
799 idiag%id_mdt = register_diag_field( trim(field),
'mdt', axes(1:3), time, &
801 idiag%id_qdt = register_diag_field( trim(field),
'qdt', axes(1:3), time, &
803 idiag%id_dbz = register_diag_field( trim(field),
'reflectivity', axes(1:3), time, &
809 idiag%id_vort = register_diag_field( trim(field),
'vort', axes(1:3), time, &
814 idiag%id_pv = register_diag_field( trim(field),
'pv', axes(1:3), time, &
820 idiag%id_uw = register_diag_field( trim(field),
'uw', axes(1:3), time, &
822 idiag%id_vw = register_diag_field( trim(field),
'vw', axes(1:3), time, &
824 idiag%id_hw = register_diag_field( trim(field),
'hw', axes(1:3), time, &
826 idiag%id_qvw = register_diag_field( trim(field),
'qvw', axes(1:3), time, &
828 idiag%id_qlw = register_diag_field( trim(field),
'qlw', axes(1:3), time, &
830 idiag%id_qiw = register_diag_field( trim(field),
'qiw', axes(1:3), time, &
832 idiag%id_o3w = register_diag_field( trim(field),
'o3w', axes(1:3), time, &
838 idiag%id_te = register_diag_field( trim(field),
'te', axes(1:2), time, &
841 idiag%id_ke = register_diag_field( trim(field),
'ke', axes(1:2), time, &
843 idiag%id_ws = register_diag_field( trim(field),
'ws', axes(1:2), time, &
845 idiag%id_maxdbz = register_diag_field( trim(field),
'max_reflectivity', axes(1:2), time, &
847 idiag%id_basedbz = register_diag_field( trim(field),
'base_reflectivity', axes(1:2), time, &
849 idiag%id_dbz4km = register_diag_field( trim(field),
'4km_reflectivity', axes(1:2), time, &
851 idiag%id_dbztop = register_diag_field( trim(field),
'echo_top', axes(1:2), time, &
853 idiag%id_dbz_m10C = register_diag_field( trim(field),
'm10C_reflectivity', axes(1:2), time, &
860 idiag%id_vorts = register_diag_field( trim(field),
'vorts', axes(1:2), time, &
862 idiag%id_us = register_diag_field( trim(field),
'us', axes(1:2), time, &
864 idiag%id_vs = register_diag_field( trim(field),
'vs', axes(1:2), time, &
866 idiag%id_tq = register_diag_field( trim(field),
'tq', axes(1:2), time, &
868 idiag%id_iw = register_diag_field( trim(field),
'iw', axes(1:2), time, &
870 idiag%id_lw = register_diag_field( trim(field),
'lw', axes(1:2), time, &
872 idiag%id_ts = register_diag_field( trim(field),
'ts', axes(1:2), time, &
873 'Skin temperature',
'K' )
874 idiag%id_tb = register_diag_field( trim(field),
'tb', axes(1:2), time, &
875 'lowest layer temperature',
'K' )
876 idiag%id_ctt = register_diag_field( trim(field),
'ctt', axes(1:2), time, &
878 idiag%id_ctp = register_diag_field( trim(field),
'ctp', axes(1:2), time, &
880 idiag%id_ctz = register_diag_field( trim(field),
'ctz', axes(1:2), time, &
882 idiag%id_cape = register_diag_field( trim(field),
'cape', axes(1:2), time, &
884 idiag%id_cin = register_diag_field( trim(field),
'cin', axes(1:2), time, &
889 idiag%id_intqv = register_diag_field( trim(field),
'intqv', axes(1:2), time, &
891 idiag%id_intql = register_diag_field( trim(field),
'intql', axes(1:2), time, &
893 idiag%id_intqi = register_diag_field( trim(field),
'intqi', axes(1:2), time, &
895 idiag%id_intqr = register_diag_field( trim(field),
'intqr', axes(1:2), time, &
897 idiag%id_intqs = register_diag_field( trim(field),
'intqs', axes(1:2), time, &
899 idiag%id_intqg = register_diag_field( trim(field),
'intqg', axes(1:2), time, &
903 idiag%id_acl = register_diag_field( trim(field),
'acl', axes(1:2), time, &
905 idiag%id_acl2 = register_diag_field( trim(field),
'acl2', axes(1:2), time, &
907 idiag%id_acly = register_diag_field( trim(field),
'acly', axes(1:2), time, &
914 idiag%id_vort850 = register_diag_field( trim(field),
'vort850', axes(1:2), time, &
917 idiag%id_vort200 = register_diag_field( trim(field),
'vort200', axes(1:2), time, &
922 idiag%id_s200 = register_diag_field( trim(field),
's200', axes(1:2), time, &
924 idiag%id_sl12 = register_diag_field( trim(field),
'sl12', axes(1:2), time, &
926 idiag%id_sl13 = register_diag_field( trim(field),
'sl13', axes(1:2), time, &
929 idiag%id_qn200 = register_diag_field( trim(field),
'qn200', axes(1:2), time, &
931 idiag%id_qn500 = register_diag_field( trim(field),
'qn500', axes(1:2), time, &
933 idiag%id_qn850 = register_diag_field( trim(field),
'qn850', axes(1:2), time, &
936 idiag%id_vort500 = register_diag_field( trim(field),
'vort500', axes(1:2), time, &
939 idiag%id_rain5km = register_diag_field( trim(field),
'rain5km', axes(1:2), time, &
944 if( .not. atm(n)%flagstruct%hydrostatic )
then 945 idiag%id_w200 = register_diag_field( trim(field),
'w200', axes(1:2), time, &
947 idiag%id_w500 = register_diag_field( trim(field),
'w500', axes(1:2), time, &
949 idiag%id_w700 = register_diag_field( trim(field),
'w700', axes(1:2), time, &
952 idiag%id_w850 = register_diag_field( trim(field),
'w850', axes(1:2), time, &
954 idiag%id_w5km = register_diag_field( trim(field),
'w5km', axes(1:2), time, &
956 idiag%id_w2500m = register_diag_field( trim(field),
'w2500m', axes(1:2), time, &
958 idiag%id_w1km = register_diag_field( trim(field),
'w1km', axes(1:2), time, &
960 idiag%id_wmaxup = register_diag_field( trim(field),
'wmaxup', axes(1:2), time, &
962 idiag%id_wmaxdn = register_diag_field( trim(field),
'wmaxdn', axes(1:2), time, &
968 idiag%id_x850 = register_diag_field( trim(field),
'x850', axes(1:2), time, &
976 idiag%id_srh1 = register_diag_field( trim(field),
'srh01', axes(1:2), time, &
978 idiag%id_srh3 = register_diag_field( trim(field),
'srh03', axes(1:2), time, &
980 idiag%id_ustm = register_diag_field( trim(field),
'ustm', axes(1:2), time, &
982 idiag%id_vstm = register_diag_field( trim(field),
'vstm', axes(1:2), time, &
985 idiag%id_srh25 = register_diag_field( trim(field),
'srh25', axes(1:2), time, &
988 if( .not. atm(n)%flagstruct%hydrostatic )
then 989 idiag%id_uh03 = register_diag_field( trim(field),
'uh03', axes(1:2), time, &
991 idiag%id_uh25 = register_diag_field( trim(field),
'uh25', axes(1:2), time, &
995 if( .not. atm(n)%flagstruct%hydrostatic ) &
996 idiag%id_w100m = register_diag_field( trim(field),
'w100m', axes(1:2), time, &
998 idiag%id_u100m = register_diag_field( trim(field),
'u100m', axes(1:2), time, &
1000 idiag%id_v100m = register_diag_field( trim(field),
'v100m', axes(1:2), time, &
1005 idiag%id_rh10 = register_diag_field( trim(field),
'rh10', axes(1:2), time, &
1007 idiag%id_rh50 = register_diag_field( trim(field),
'rh50', axes(1:2), time, &
1009 idiag%id_rh100 = register_diag_field( trim(field),
'rh100', axes(1:2), time, &
1011 idiag%id_rh200 = register_diag_field( trim(field),
'rh200', axes(1:2), time, &
1013 idiag%id_rh250 = register_diag_field( trim(field),
'rh250', axes(1:2), time, &
1015 idiag%id_rh300 = register_diag_field( trim(field),
'rh300', axes(1:2), time, &
1017 idiag%id_rh500 = register_diag_field( trim(field),
'rh500', axes(1:2), time, &
1019 idiag%id_rh700 = register_diag_field( trim(field),
'rh700', axes(1:2), time, &
1021 idiag%id_rh850 = register_diag_field( trim(field),
'rh850', axes(1:2), time, &
1023 idiag%id_rh925 = register_diag_field( trim(field),
'rh925', axes(1:2), time, &
1025 idiag%id_rh1000 = register_diag_field( trim(field),
'rh1000', axes(1:2), time, &
1030 idiag%id_dp10 = register_diag_field( trim(field),
'dp10', axes(1:2), time, &
1032 idiag%id_dp50 = register_diag_field( trim(field),
'dp50', axes(1:2), time, &
1034 idiag%id_dp100 = register_diag_field( trim(field),
'dp100', axes(1:2), time, &
1036 idiag%id_dp200 = register_diag_field( trim(field),
'dp200', axes(1:2), time, &
1038 idiag%id_dp250 = register_diag_field( trim(field),
'dp250', axes(1:2), time, &
1040 idiag%id_dp300 = register_diag_field( trim(field),
'dp300', axes(1:2), time, &
1042 idiag%id_dp500 = register_diag_field( trim(field),
'dp500', axes(1:2), time, &
1044 idiag%id_dp700 = register_diag_field( trim(field),
'dp700', axes(1:2), time, &
1046 idiag%id_dp850 = register_diag_field( trim(field),
'dp850', axes(1:2), time, &
1048 idiag%id_dp925 = register_diag_field( trim(field),
'dp925', axes(1:2), time, &
1050 idiag%id_dp1000 = register_diag_field( trim(field),
'dp1000', axes(1:2), time, &
1055 idiag%id_rh10_cmip = register_diag_field( trim(field),
'rh10_cmip', axes(1:2), time, &
1057 idiag%id_rh50_cmip = register_diag_field( trim(field),
'rh50_cmip', axes(1:2), time, &
1059 idiag%id_rh100_cmip = register_diag_field( trim(field),
'rh100_cmip', axes(1:2), time, &
1061 idiag%id_rh250_cmip = register_diag_field( trim(field),
'rh250_cmip', axes(1:2), time, &
1063 idiag%id_rh300_cmip = register_diag_field( trim(field),
'rh300_cmip', axes(1:2), time, &
1065 idiag%id_rh500_cmip = register_diag_field( trim(field),
'rh500_cmip', axes(1:2), time, &
1067 idiag%id_rh700_cmip = register_diag_field( trim(field),
'rh700_cmip', axes(1:2), time, &
1069 idiag%id_rh850_cmip = register_diag_field( trim(field),
'rh850_cmip', axes(1:2), time, &
1071 idiag%id_rh925_cmip = register_diag_field( trim(field),
'rh925_cmip', axes(1:2), time, &
1073 idiag%id_rh1000_cmip = register_diag_field( trim(field),
'rh1000_cmip', axes(1:2), time, &
1076 if (atm(n)%flagstruct%write_3d_diags)
then 1082 idiag%id_tracer(i) = register_diag_field( field, trim(
tname), &
1086 if (
idiag%id_tracer(i) > 0)
then 1088 write(unit,
'(a,a,a,a)') &
1089 &
'Diagnostics available for tracer ',trim(
tname), &
1090 ' in module ', trim(field)
1098 if (trim(
tname).eq.
'co2')
then 1099 idiag%w_mr(:) = wtmco2
1100 idiag%id_tracer_dmmr(i) = register_diag_field( field, trim(
tname)//
'_dmmr', &
1101 axes(1:3), time, trim(
tlongname)//
" (dry mmr)", &
1103 idiag%id_tracer_dvmr(i) = register_diag_field( field, trim(
tname)//
'_dvmr', &
1104 axes(1:3), time, trim(
tlongname)//
" (dry vmr)", &
1108 if (
idiag%id_tracer_dmmr(i) > 0)
then 1109 write(unit,
'(a,a,a,a)')
'Diagnostics available for '//trim(
tname)//
' dry mmr ', &
1110 trim(
tname)//
'_dmmr',
' in module ', trim(field)
1112 if (
idiag%id_tracer_dvmr(i) > 0)
then 1113 write(unit,
'(a,a,a,a)')
'Diagnostics available for '//trim(
tname)//
' dry vmr ', &
1114 trim(
tname)//
'_dvmr',
' in module ', trim(field)
1123 if ( atm(1)%flagstruct%consv_am .or.
idiag%id_mq > 0 .or.
idiag%id_amdt > 0 )
then 1124 allocate (
idiag%zxg(isc:iec,jsc:jec) )
1126 call init_mq(atm(n)%phis, atm(n)%gridstruct, &
1127 npx, npy, isc, iec, jsc, jec, atm(n)%ng)
1134 call prt_mass(npz, atm(n)%ncnst, isc, iec, jsc, jec, atm(n)%ng, max(1,atm(n)%flagstruct%nwat), &
1135 atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
1137 call prt_mass(npz, atm(n)%ncnst, isc, iec, jsc, jec, atm(n)%ng, atm(n)%flagstruct%nwat, &
1138 atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
1162 #ifdef INTERNAL_FILE_NML 1163 read(input_nml_file, nml=fv_diag_column_nml,iostat=ios)
1165 inquire (file=trim(atm(n)%nml_filename), exist=exists)
1166 if (.not. exists)
then 1167 write(errmsg,*)
'fv_diag_column_nml: namelist file ',trim(atm(n)%nml_filename),
' does not exist' 1168 call mpp_error(fatal, errmsg)
1170 open (unit=nlunit, file=atm(n)%nml_filename, readonly, status=
'OLD', iostat=ios)
1173 read (nlunit, nml=fv_diag_column_nml, iostat=ios)
1177 call column_diagnostics_init
1203 call initialize_diagnostic_columns(
"DEBUG", num_diag_pts_latlon=
num_diag_debug, num_diag_pts_ij=0, &
1204 global_i=(/1/), global_j=(/1/), &
1206 lonb_in=atm(n)%gridstruct%agrid(isc:iec,jsc:jec,1), latb_in=atm(n)%gridstruct%agrid(isc:iec,jsc:jec,2), &
1216 write(*,
'(A, 1x, I04, 1x, A, 4F7.2, 2I5)')
'DEBUG POINT: ', mpp_pe(),
diag_debug_names(m),
diag_debug_lon_in(m),
diag_debug_lat_in(m), &
1249 call initialize_diagnostic_columns(
"Sounding", num_diag_pts_latlon=
num_diag_sonde, num_diag_pts_ij=0, &
1250 global_i=(/1/), global_j=(/1/), &
1252 lonb_in=atm(n)%gridstruct%agrid(isc:iec,jsc:jec,1), latb_in=atm(n)%gridstruct%agrid(isc:iec,jsc:jec,2), &
1262 write(*,
'(A, 1x, I04, 1x, A, 4F7.2, 2I5)')
'SONDE POINT: ', mpp_pe(),
diag_sonde_names(m),
diag_sonde_lon_in(m),
diag_sonde_lat_in(m), &
1279 call nullify_domain()
1284 if(
idiag%id_theta_e >0 )
call qsmith_init
1289 subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
1290 integer,
intent(in):: npx, npy, is, ie, js, je, ng
1291 real,
intent(in):: phis(is-ng:ie+ng, js-ng:je+ng)
1295 real zs(is-ng:ie+ng, js-ng:je+ng)
1296 real zb(is-ng:ie+ng, js-ng:je+ng)
1297 real pdx(3,is:ie,js:je+1)
1298 real pdy(3,is:ie+1,js:je)
1301 real,
pointer :: rarea(:,:)
1302 real,
pointer,
dimension(:,:) :: dx, dy
1303 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: en1, en2, vlon, vlat
1304 real,
pointer,
dimension(:,:,:) :: agrid
1306 rarea => gridstruct%rarea
1309 en1 => gridstruct%en1
1310 en2 => gridstruct%en2
1311 agrid => gridstruct%agrid
1312 vlon => gridstruct%vlon
1313 vlat => gridstruct%vlat
1319 zs(i,j) = phis(i,j) / grav
1325 call a2b_ord4(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)
1330 pdx(n,i,j) = 0.5*(zb(i,j)+zb(i+1,j))*dx(i,j)*en1(n,i,j)
1337 pdy(n,i,j) = 0.5*(zb(i,j)+zb(i,j+1))*dy(i,j)*en2(n,i,j)
1345 idiag%zxg(i,j) = vlon(i,j,1)*(pdx(1,i,j+1)-pdx(1,i,j)-pdy(1,i,j)+pdy(1,i+1,j)) &
1346 + vlon(i,j,2)*(pdx(2,i,j+1)-pdx(2,i,j)-pdy(2,i,j)+pdy(2,i+1,j)) &
1347 + vlon(i,j,3)*(pdx(3,i,j+1)-pdx(3,i,j)-pdy(3,i,j)+pdy(3,i+1,j))
1350 idiag%zxg(i,j) =
idiag%zxg(i,j)*rarea(i,j) * radius*cos(agrid(i,j,2))
1357 subroutine fv_diag(Atm, zvir, Time, print_freq)
1360 type(time_type),
intent(in) :: Time
1361 real,
intent(in):: zvir
1362 integer,
intent(in):: print_freq
1364 integer :: isc, iec, jsc, jec, n, ntileMe
1365 integer :: isd, ied, jsd, jed, npz, itrac
1366 integer :: ngc, nwater
1368 real,
allocatable :: a2(:,:),a3(:,:,:), wk(:,:,:), wz(:,:,:), ucoor(:,:,:), vcoor(:,:,:)
1369 real,
allocatable :: ustm(:,:), vstm(:,:)
1370 real,
allocatable :: slp(:,:), depress(:,:), ws_max(:,:), tc_count(:,:)
1371 real,
allocatable :: u2(:,:), v2(:,:), x850(:,:), var1(:,:), var2(:,:), var3(:,:)
1372 real,
allocatable :: dmmr(:,:,:), dvmr(:,:,:)
1376 real :: tot_mq, tmp, sar, slon, slat
1377 real :: a1d(atm(1)%npz)
1379 logical :: do_cs_intp
1381 logical :: bad_range
1382 integer i,j,k, yr, mon, dd, hr, mn, days, seconds, nq, theta_d
1383 character(len=128) :: tname
1384 real,
parameter:: ws_0 = 16.
1385 real,
parameter:: ws_1 = 20.
1386 real,
parameter:: vort_c0= 2.2e-5
1387 logical,
allocatable :: storm(:,:), cat_crt(:,:)
1388 real :: tmp2, pvsum, e2, einf, qm, mm, maxdbz, allmax, rgrav, cv_vapor
1389 real,
allocatable :: cvm(:)
1390 integer :: Cl, Cl2, k1, k2
1403 pout(i) =
levs(i) * 1.e2
1404 plevs(i) = log( pout(i) )
1407 ntileme =
size(atm(:))
1409 isc = atm(n)%bd%isc; iec = atm(n)%bd%iec
1410 jsc = atm(n)%bd%jsc; jec = atm(n)%bd%jec
1414 nq =
size (atm(n)%q,4)
1416 isd = atm(n)%bd%isd; ied = atm(n)%bd%ied
1417 jsd = atm(n)%bd%jsd; jed = atm(n)%bd%jed
1419 if(
idiag%id_c15>0 )
then 1420 allocate ( storm(isc:iec,jsc:jec) )
1421 allocate ( depress(isc:iec,jsc:jec) )
1422 allocate ( ws_max(isc:iec,jsc:jec) )
1423 allocate ( cat_crt(isc:iec,jsc:jec) )
1424 allocate (tc_count(isc:iec,jsc:jec) )
1427 if(
idiag%id_x850>0 )
then 1428 allocate ( x850(isc:iec,jsc:jec) )
1432 call set_domain(atm(1)%domain)
1435 call get_date(
fv_time, yr, mon, dd, hr, mn, seconds)
1436 if( print_freq == 0 )
then 1438 elseif( print_freq < 0 )
then 1442 prt_minmax = mod(hr, print_freq) == 0 .and. mn==0 .and. seconds==0
1451 call get_time (
fv_time, seconds, days)
1452 if( print_freq == 0 )
then 1454 elseif( print_freq < 0 )
then 1458 prt_minmax = mod(seconds, 3600*print_freq) == 0
1471 if(
master)
write(*,*) yr, mon, dd, hr, mn, seconds
1473 if(
master)
write(*,*) days, seconds
1477 allocate ( a2(isc:iec,jsc:jec) )
1481 call prt_mxm(
'ZS',
idiag%zsurf, isc, iec, jsc, jec, 0, 1, 1.0, atm(n)%gridstruct%area_64, atm(n)%domain)
1482 call prt_maxmin(
'PS', atm(n)%ps, isc, iec, jsc, jec, ngc, 1, 0.01)
1485 allocate(var2(isc:iec,jsc:jec))
1489 slat =
rad2deg*atm(n)%gridstruct%agrid(i,j,2)
1490 if (slat >= 0.)
then 1491 a2(i,j) = atm(n)%ps(i,j)
1495 var2(i,j) = atm(n)%ps(i,j)
1499 call prt_maxmin(
'NH PS', a2, isc, iec, jsc, jec, 0, 1, 0.01)
1500 call prt_maxmin(
'SH PS', var2, isc, iec, jsc, jec, 0, 1, 0.01)
1506 call prt_mass(npz, nq, isc, iec, jsc, jec, ngc, max(1,atm(n)%flagstruct%nwat), &
1507 atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
1509 call prt_mass(npz, nq, isc, iec, jsc, jec, ngc, atm(n)%flagstruct%nwat, &
1510 atm(n)%ps, atm(n)%delp, atm(n)%q, atm(n)%gridstruct%area_64, atm(n)%domain)
1514 if (atm(n)%flagstruct%consv_te > 1.e-5)
then 1519 write(*,*)
'ENG Deficit (W/m**2)', trim(
gn),
'=',
e_flux 1524 if ( .not. atm(n)%flagstruct%hydrostatic ) &
1526 atm(n)%w, atm(n)%delz, atm(n)%pt, atm(n)%delp, &
1527 atm(n)%q, atm(n)%phis, atm(n)%gridstruct%area, atm(n)%domain, &
1529 atm(n)%ua, atm(n)%va, atm(n)%flagstruct%moist_phys, a2)
1531 call prt_maxmin(
'UA_top', atm(n)%ua(isc:iec,jsc:jec,1), &
1532 isc, iec, jsc, jec, 0, 1, 1.)
1533 call prt_maxmin(
'UA', atm(n)%ua, isc, iec, jsc, jec, ngc, npz, 1.)
1534 call prt_maxmin(
'VA', atm(n)%va, isc, iec, jsc, jec, ngc, npz, 1.)
1536 if ( .not. atm(n)%flagstruct%hydrostatic )
then 1537 call prt_maxmin(
'W ', atm(n)%w , isc, iec, jsc, jec, ngc, npz, 1.)
1538 call prt_maxmin(
'Bottom w', atm(n)%w(isc:iec,jsc:jec,npz), isc, iec, jsc, jec, 0, 1, 1.)
1541 a2(i,j) = -atm(n)%w(i,j,npz)/atm(n)%delz(i,j,npz)
1544 call prt_maxmin(
'Bottom: w/dz', a2, isc, iec, jsc, jec, 0, 1, 1.)
1546 if ( atm(n)%flagstruct%hybrid_z )
call prt_maxmin(
'Hybrid_ZTOP (km)', atm(n)%ze0(isc:iec,jsc:jec,1), &
1547 isc, iec, jsc, jec, 0, 1, 1.e-3)
1548 call prt_maxmin(
'DZ (m)', atm(n)%delz(isc:iec,jsc:jec,1:npz), &
1549 isc, iec, jsc, jec, 0, npz, 1.)
1550 call prt_maxmin(
'Bottom DZ (m)', atm(n)%delz(isc:iec,jsc:jec,npz), &
1551 isc, iec, jsc, jec, 0, 1, 1.)
1557 call prt_maxmin(
'TA', atm(n)%pt, isc, iec, jsc, jec, ngc, npz, 1.)
1560 call prt_maxmin(
'OM', atm(n)%omga, isc, iec, jsc, jec, ngc, npz, 1.)
1563 elseif ( atm(n)%flagstruct%range_warn )
then 1564 call range_check(
'DELP', atm(n)%delp, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid, &
1565 0.01*
ptop, 200.e2, bad_range, time)
1566 call range_check(
'UA', atm(n)%ua, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid, &
1567 -250., 250., bad_range, time)
1568 call range_check(
'VA', atm(n)%va, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid, &
1569 -250., 250., bad_range, time)
1571 call range_check(
'TA', atm(n)%pt, isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid, &
1573 130., 350., bad_range, time)
1575 150., 350., bad_range, time)
1578 call range_check(
'Qv', atm(n)%q(:,:,:,
sphum), isc, iec, jsc, jec, ngc, npz, atm(n)%gridstruct%agrid, &
1579 -1.e-8, 1.e20, bad_range, time)
1583 allocate ( u2(isc:iec,jsc:jec) )
1584 allocate ( v2(isc:iec,jsc:jec) )
1585 allocate ( wk(isc:iec,jsc:jec,npz) )
1586 if ( any(
idiag%id_tracer_dmmr > 0) .or. any(
idiag%id_tracer_dvmr > 0) )
then 1587 allocate ( dmmr(isc:iec,jsc:jec,1:npz) )
1588 allocate ( dvmr(isc:iec,jsc:jec,1:npz) )
1595 if(
idiag%id_zsurf > 0) used=send_data(
idiag%id_zsurf,
idiag%zsurf, time)
1597 if(
idiag%id_ps > 0) used=send_data(
idiag%id_ps, atm(n)%ps(isc:iec,jsc:jec), time)
1599 if (
idiag%id_qv_dt_phys > 0) used=send_data(
idiag%id_qv_dt_phys, atm(n)%phys_diag%phys_qv_dt(isc:iec,jsc:jec,1:npz), time)
1600 if (
idiag%id_ql_dt_phys > 0) used=send_data(
idiag%id_ql_dt_phys, atm(n)%phys_diag%phys_ql_dt(isc:iec,jsc:jec,1:npz), time)
1601 if (
idiag%id_qi_dt_phys > 0) used=send_data(
idiag%id_qi_dt_phys, atm(n)%phys_diag%phys_qi_dt(isc:iec,jsc:jec,1:npz), time)
1602 if (
idiag%id_t_dt_phys > 0) used=send_data(
idiag%id_t_dt_phys, atm(n)%phys_diag%phys_t_dt(isc:iec,jsc:jec,1:npz), time)
1603 if (
idiag%id_u_dt_phys > 0) used=send_data(
idiag%id_u_dt_phys, atm(n)%phys_diag%phys_u_dt(isc:iec,jsc:jec,1:npz), time)
1604 if (
idiag%id_v_dt_phys > 0) used=send_data(
idiag%id_v_dt_phys, atm(n)%phys_diag%phys_v_dt(isc:iec,jsc:jec,1:npz), time)
1607 call wind_max(isc, iec, jsc, jec ,isd, ied, jsd, jed, atm(n)%ua(isc:iec,jsc:jec,npz), &
1608 atm(n)%va(isc:iec,jsc:jec,npz), ws_max, atm(n)%domain)
1611 if( abs(atm(n)%gridstruct%agrid(i,j,2)*
rad2deg)<45.0 .and. &
1612 atm(n)%phis(i,j)*
ginv<500.0 .and. ws_max(i,j)>ws_0 )
then 1615 storm(i,j) = .false.
1621 if (
idiag%id_vort200>0 .or.
idiag%id_vort500>0 .or.
idiag%id_vort850>0 .or.
idiag%id_vorts>0 &
1624 call get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, atm(n)%u, atm(n)%v, wk, &
1625 atm(n)%gridstruct%dx, atm(n)%gridstruct%dy, atm(n)%gridstruct%rarea)
1627 if(
idiag%id_vort >0) used=send_data(
idiag%id_vort, wk, time)
1628 if(
idiag%id_vorts>0) used=send_data(
idiag%id_vorts, wk(isc:iec,jsc:jec,npz), time)
1630 if(
idiag%id_c15>0)
then 1634 storm(i,j) = (atm(n)%gridstruct%agrid(i,j,2)>0. .and. wk(i,j,npz)> vort_c0) .or. &
1635 (atm(n)%gridstruct%agrid(i,j,2)<0. .and. wk(i,j,npz)<-vort_c0)
1640 if(
idiag%id_vort200>0 )
then 1642 200.e2, atm(n)%peln, wk, a2)
1643 used=send_data(
idiag%id_vort200, a2, time)
1645 if(
idiag%id_vort500>0 )
then 1647 500.e2, atm(n)%peln, wk, a2)
1648 used=send_data(
idiag%id_vort500, a2, time)
1653 850.e2, atm(n)%peln, wk, a2)
1654 used=send_data(
idiag%id_vort850, a2, time)
1655 if (
idiag%id_x850>0 ) x850(:,:) = a2(:,:)
1657 if(
idiag%id_c15>0)
then 1661 storm(i,j) = (atm(n)%gridstruct%agrid(i,j,2)>0. .and. a2(i,j)> vort_c0) .or. &
1662 (atm(n)%gridstruct%agrid(i,j,2)<0. .and. a2(i,j)<-vort_c0)
1669 if( .not. atm(n)%flagstruct%hydrostatic )
then 1671 if (
idiag%id_uh03 > 0 )
then 1673 atm(n)%w, wk, atm(n)%delz, atm(n)%q, &
1674 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 3.e3)
1675 used = send_data(
idiag%id_uh03, a2, time )
1679 tmp =
rad2deg * atm(n)%gridstruct%agrid(i,j,1)
1680 tmp2 =
rad2deg * atm(n)%gridstruct%agrid(i,j,2)
1681 if ( tmp2<25. .or. tmp2>50. &
1682 .or. tmp<235. .or. tmp>300. )
then 1687 call prt_maxmin(
'UH over CONUS', a2, isc, iec, jsc, jec, 0, 1, 1.)
1690 if (
idiag%id_uh25 > 0 )
then 1692 atm(n)%w, wk, atm(n)%delz, atm(n)%q, &
1693 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 2.e3, 5.e3)
1694 used = send_data(
idiag%id_uh25, a2, time )
1699 if (
idiag%id_srh1 > 0 .or.
idiag%id_srh3 > 0 .or.
idiag%id_srh25 > 0 .or.
idiag%id_ustm > 0 .or.
idiag%id_vstm > 0)
then 1700 allocate(ustm(isc:iec,jsc:jec), vstm(isc:iec,jsc:jec))
1702 call bunkers_vector(isc, iec, jsc, jec, ngc, npz, zvir,
sphum, ustm, vstm, &
1703 atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q, &
1704 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav)
1706 if (
idiag%id_ustm > 0 )
then 1707 used = send_data(
idiag%id_ustm, ustm, time )
1709 if (
idiag%id_vstm > 0 )
then 1710 used = send_data(
idiag%id_vstm, vstm, time )
1713 if (
idiag%id_srh1 > 0 )
then 1714 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir,
sphum, a2, ustm, vstm, &
1715 atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q, &
1716 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 1.e3)
1717 used = send_data(
idiag%id_srh1, a2, time )
1721 tmp =
rad2deg * atm(n)%gridstruct%agrid(i,j,1)
1722 tmp2 =
rad2deg * atm(n)%gridstruct%agrid(i,j,2)
1723 if ( tmp2<25. .or. tmp2>50. &
1724 .or. tmp<235. .or. tmp>300. )
then 1729 call prt_maxmin(
'SRH (0-1 km) over CONUS', a2, isc, iec, jsc, jec, 0, 1, 1.)
1733 if (
idiag%id_srh3 > 0 )
then 1734 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir,
sphum, a2, ustm, vstm, &
1735 atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q, &
1736 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 0., 3e3)
1737 used = send_data(
idiag%id_srh3, a2, time )
1741 tmp =
rad2deg * atm(n)%gridstruct%agrid(i,j,1)
1742 tmp2 =
rad2deg * atm(n)%gridstruct%agrid(i,j,2)
1743 if ( tmp2<25. .or. tmp2>50. &
1744 .or. tmp<235. .or. tmp>300. )
then 1749 call prt_maxmin(
'SRH (0-3 km) over CONUS', a2, isc, iec, jsc, jec, 0, 1, 1.)
1753 if (
idiag%id_srh25 > 0 )
then 1754 call helicity_relative_caps(isc, iec, jsc, jec, ngc, npz, zvir,
sphum, a2, ustm, vstm, &
1755 atm(n)%ua, atm(n)%va, atm(n)%delz, atm(n)%q, &
1756 atm(n)%flagstruct%hydrostatic, atm(n)%pt, atm(n)%peln, atm(n)%phis, grav, 2.e3, 5e3)
1757 used = send_data(
idiag%id_srh25, a2, time )
1761 tmp =
rad2deg * atm(n)%gridstruct%agrid(i,j,1)
1762 tmp2 =
rad2deg * atm(n)%gridstruct%agrid(i,j,2)
1763 if ( tmp2<25. .or. tmp2>50. &
1764 .or. tmp<235. .or. tmp>300. )
then 1769 call prt_maxmin(
'SRH (2-5 km) over CONUS', a2, isc, iec, jsc, jec, 0, 1, 1.)
1773 deallocate(ustm, vstm)
1777 if (
idiag%id_pv > 0 )
then 1779 call pv_entropy(isc, iec, jsc, jec, ngc, npz, wk, &
1780 atm(n)%gridstruct%f0, atm(n)%pt, atm(n)%pkz, atm(n)%delp, grav)
1781 used = send_data(
idiag%id_pv, wk, time )
1818 if (
idiag%id_rh > 0 )
then 1823 a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
1827 call qsmith((iec-isc+1)*(jec-jsc+1), npz, &
1828 (iec-isc+1)*(jec-jsc+1), 1, atm(n)%pt(isc:iec,jsc:jec,k), &
1829 a2, atm(n)%q(isc:iec,jsc:jec,k,
sphum), wk(isc,jsc,k))
1831 call qsmith(iec-isc+1, jec-jsc+1, 1, atm(n)%pt(isc:iec,jsc:jec,k), &
1832 a2, atm(n)%q(isc:iec,jsc:jec,k,
sphum), wk(isc,jsc,k))
1836 wk(i,j,k) = 100.*atm(n)%q(i,j,k,
sphum)/wk(i,j,k)
1840 used = send_data(
idiag%id_rh, wk, time )
1842 call prt_maxmin(
'RH_sf (%)', wk(isc:iec,jsc:jec,npz), isc, iec, jsc, jec, 0, 1, 1.)
1843 call prt_maxmin(
'RH_3D (%)', wk, isc, iec, jsc, jec, 0, npz, 1.)
1850 idiag%id_rh925>0 .or.
idiag%id_rh1000>0 .or. &
1853 idiag%id_dp925>0 .or.
idiag%id_dp1000>0)
then 1858 a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
1861 call rh_calc (a2, atm(n)%pt(isc:iec,jsc:jec,k), &
1862 atm(n)%q(isc:iec,jsc:jec,k,
sphum), wk(isc:iec,jsc:jec,k))
1864 if (
idiag%id_rh50>0)
then 1865 call interpolate_vertical(isc, iec, jsc, jec, npz, 50.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1866 used=send_data(
idiag%id_rh50, a2, time)
1868 if (
idiag%id_rh100>0)
then 1869 call interpolate_vertical(isc, iec, jsc, jec, npz, 100.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1870 used=send_data(
idiag%id_rh100, a2, time)
1872 if (
idiag%id_rh200>0)
then 1873 call interpolate_vertical(isc, iec, jsc, jec, npz, 200.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1874 used=send_data(
idiag%id_rh200, a2, time)
1876 if (
idiag%id_rh250>0)
then 1877 call interpolate_vertical(isc, iec, jsc, jec, npz, 250.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1878 used=send_data(
idiag%id_rh250, a2, time)
1880 if (
idiag%id_rh300>0)
then 1881 call interpolate_vertical(isc, iec, jsc, jec, npz, 300.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1882 used=send_data(
idiag%id_rh300, a2, time)
1884 if (
idiag%id_rh500>0)
then 1885 call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1886 used=send_data(
idiag%id_rh500, a2, time)
1888 if (
idiag%id_rh700>0)
then 1889 call interpolate_vertical(isc, iec, jsc, jec, npz, 700.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1890 used=send_data(
idiag%id_rh700, a2, time)
1892 if (
idiag%id_rh850>0)
then 1893 call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1894 used=send_data(
idiag%id_rh850, a2, time)
1896 if (
idiag%id_rh925>0)
then 1897 call interpolate_vertical(isc, iec, jsc, jec, npz, 925.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1898 used=send_data(
idiag%id_rh925, a2, time)
1900 if (
idiag%id_rh1000>0)
then 1901 call interpolate_vertical(isc, iec, jsc, jec, npz, 1000.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1902 used=send_data(
idiag%id_rh1000, a2, time)
1907 idiag%id_dp925>0 .or.
idiag%id_dp1000>0 )
then 1909 if (
allocated(a3))
deallocate(a3)
1910 allocate(a3(isc:iec,jsc:jec,1:npz))
1916 tmp = ( log(max(wk(i,j,k)*1.e-2,1.e-2)) + 17.27 * ( atm(n)%pt(i,j,k) - 273.14 )/ ( -35.84 + atm(n)%pt(i,j,k)) ) / 17.27
1917 a3(i,j,k) = 273.14 + 237.3*tmp/ ( 1. - tmp )
1922 if (
idiag%id_dp50>0)
then 1924 used=send_data(
idiag%id_dp50, a2, time)
1926 if (
idiag%id_dp100>0)
then 1928 used=send_data(
idiag%id_dp100, a2, time)
1930 if (
idiag%id_dp200>0)
then 1932 used=send_data(
idiag%id_dp200, a2, time)
1934 if (
idiag%id_dp250>0)
then 1936 used=send_data(
idiag%id_dp250, a2, time)
1938 if (
idiag%id_dp300>0)
then 1940 used=send_data(
idiag%id_dp300, a2, time)
1942 if (
idiag%id_dp500>0)
then 1944 used=send_data(
idiag%id_dp500, a2, time)
1946 if (
idiag%id_dp700>0)
then 1948 used=send_data(
idiag%id_dp700, a2, time)
1950 if (
idiag%id_dp850>0)
then 1952 used=send_data(
idiag%id_dp850, a2, time)
1954 if (
idiag%id_dp925>0)
then 1956 used=send_data(
idiag%id_dp925, a2, time)
1958 if (
idiag%id_dp1000>0)
then 1960 used=send_data(
idiag%id_dp1000, a2, time)
1969 if (
idiag%id_rh10_cmip>0 .or.
idiag%id_rh50_cmip>0 .or.
idiag%id_rh100_cmip>0 .or. &
1970 idiag%id_rh250_cmip>0 .or.
idiag%id_rh300_cmip>0 .or.
idiag%id_rh500_cmip>0 .or. &
1971 idiag%id_rh700_cmip>0 .or.
idiag%id_rh850_cmip>0 .or.
idiag%id_rh925_cmip>0 .or. &
1972 idiag%id_rh1000_cmip>0)
then 1977 a2(i,j) = atm(n)%delp(i,j,k)/(atm(n)%peln(i,k+1,j)-atm(n)%peln(i,k,j))
1980 call rh_calc (a2, atm(n)%pt(isc:iec,jsc:jec,k), &
1981 atm(n)%q(isc:iec,jsc:jec,k,
sphum), wk(isc:iec,jsc:jec,k), do_cmip=.true.)
1983 if (
idiag%id_rh10_cmip>0)
then 1984 call interpolate_vertical(isc, iec, jsc, jec, npz, 10.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1985 used=send_data(
idiag%id_rh10_cmip, a2, time)
1987 if (
idiag%id_rh50_cmip>0)
then 1988 call interpolate_vertical(isc, iec, jsc, jec, npz, 50.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1989 used=send_data(
idiag%id_rh50_cmip, a2, time)
1991 if (
idiag%id_rh100_cmip>0)
then 1992 call interpolate_vertical(isc, iec, jsc, jec, npz, 100.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1993 used=send_data(
idiag%id_rh100_cmip, a2, time)
1995 if (
idiag%id_rh250_cmip>0)
then 1996 call interpolate_vertical(isc, iec, jsc, jec, npz, 250.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
1997 used=send_data(
idiag%id_rh250_cmip, a2, time)
1999 if (
idiag%id_rh300_cmip>0)
then 2000 call interpolate_vertical(isc, iec, jsc, jec, npz, 300.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2001 used=send_data(
idiag%id_rh300_cmip, a2, time)
2003 if (
idiag%id_rh500_cmip>0)
then 2004 call interpolate_vertical(isc, iec, jsc, jec, npz, 500.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2005 used=send_data(
idiag%id_rh500_cmip, a2, time)
2007 if (
idiag%id_rh700_cmip>0)
then 2008 call interpolate_vertical(isc, iec, jsc, jec, npz, 700.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2009 used=send_data(
idiag%id_rh700_cmip, a2, time)
2011 if (
idiag%id_rh850_cmip>0)
then 2012 call interpolate_vertical(isc, iec, jsc, jec, npz, 850.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2013 used=send_data(
idiag%id_rh850_cmip, a2, time)
2015 if (
idiag%id_rh925_cmip>0)
then 2016 call interpolate_vertical(isc, iec, jsc, jec, npz, 925.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2017 used=send_data(
idiag%id_rh925_cmip, a2, time)
2019 if (
idiag%id_rh1000_cmip>0)
then 2020 call interpolate_vertical(isc, iec, jsc, jec, npz, 1000.e2, atm(n)%peln, wk(isc:iec,jsc:jec,:), a2)
2021 used=send_data(
idiag%id_rh1000_cmip, a2, time)
2028 if ( storm(i,j) .and. ws_max(i,j)>ws_1 )
then 2029 cat_crt(i,j) = .true.
2031 cat_crt(i,j) = .false.
2041 allocate ( wz(isc:iec,jsc:jec,npz+1) )
2042 call get_height_field(isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%delz, &
2043 wz, atm(n)%pt, atm(n)%q, atm(n)%peln, zvir)
2045 call prt_mxm(
'ZTOP',wz(isc:iec,jsc:jec,1), isc, iec, jsc, jec, 0, 1, 1.e-3, atm(n)%gridstruct%area_64, atm(n)%domain)
2048 if (
idiag%id_hght3d > 0)
then 2049 used = send_data(
idiag%id_hght3d, 0.5*(wz(isc:iec,jsc:jec,1:npz)+wz(isc:iec,jsc:jec,2:npz+1)), time)
2052 if(
idiag%id_slp > 0)
then 2054 allocate ( slp(isc:iec,jsc:jec) )
2056 atm(n)%pt(:,:,npz), atm(n)%peln, slp, 0.01)
2058 if ( atm(n)%flagstruct%range_warn )
then 2059 call range_check(
'SLP', slp, isc, iec, jsc, jec, 0, atm(n)%gridstruct%agrid, &
2062 used = send_data(
idiag%id_slp, slp, time)
2064 call prt_maxmin(
'SLP', slp, isc, iec, jsc, jec, 0, 1, 1.)
2069 slon =
rad2deg*atm(n)%gridstruct%agrid(i,j,1)
2070 slat =
rad2deg*atm(n)%gridstruct%agrid(i,j,2)
2071 if ( slat>15. .and. slat<40. .and. slon>270. .and. slon<290. )
then 2076 call prt_maxmin(
'ATL SLP', a2, isc, iec, jsc, jec, 0, 1, 1.)
2083 allocate( a3(isc:iec,jsc:jec,
nplev) )
2085 idg(:) =
idiag%id_h(:)
2087 if (
idiag%id_tm>0 )
then 2088 idg(minloc(abs(
levs-300))) = 1
2089 idg(minloc(abs(
levs-500))) = 1
2091 idg(minloc(abs(
levs-300))) =
idiag%id_h(minloc(abs(
levs-300)))
2092 idg(minloc(abs(
levs-500))) =
idiag%id_h(minloc(abs(
levs-500)))
2095 call get_height_given_pressure(isc, iec, jsc, jec, npz, wz,
nplev, idg, plevs, atm(n)%peln, a3)
2097 idg(minloc(abs(
levs-300))) =
idiag%id_h(minloc(abs(
levs-300)))
2098 idg(minloc(abs(
levs-500))) =
idiag%id_h(minloc(abs(
levs-500)))
2101 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
2104 if (
idiag%id_h_plev>0)
then 2106 call get_height_given_pressure(isc, iec, jsc, jec, npz, wz,
nplev, id1, plevs, atm(n)%peln, a3)
2107 used=send_data(
idiag%id_h_plev, a3(isc:iec,jsc:jec,:), time)
2112 if(all(
idiag%id_h(minloc(abs(
levs-100)))>0)) &
2113 call prt_mxm(
'Z100',a3(isc:iec,jsc:jec,
k100),isc,iec,jsc,jec,0,1,1.e-3,atm(n)%gridstruct%area_64,atm(n)%domain)
2115 if(all(
idiag%id_h(minloc(abs(
levs-500)))>0))
then 2116 if (atm(n)%gridstruct%bounded_domain)
then 2117 call prt_mxm(
'Z500',a3(isc:iec,jsc:jec,
k500),isc,iec,jsc,jec,0,1,1.,atm(n)%gridstruct%area_64,atm(n)%domain)
2119 call prt_gb_nh_sh(
'fv_GFS Z500', isc,iec, jsc,jec, a3(isc,jsc,
k500), atm(n)%gridstruct%area_64(isc:iec,jsc:jec), &
2120 atm(n)%gridstruct%agrid_64(isc:iec,jsc:jec,2))
2127 if(
idiag%id_tm>0 )
then 2131 if (abs(
levs(k)-500.) < 1.)
then 2137 if (abs(
levs(k)-300.) < 1.)
then 2142 if (k1 <= 0 .or. k2 <= 0)
then 2143 call mpp_error(note,
"Could not find levs for 300--500 mb mean temperature, setting to -1")
2148 a2(i,j) = grav*(a3(i,j,k2)-a3(i,j,k1))/(rdgas*(plevs(k1)-plevs(k2)))
2152 used = send_data(
idiag%id_tm, a2, time )
2159 if ( storm(i,j) )
then 2160 if( a2(i,j)<254.0 .or. atm(n)%pt(i,j,npz)<281.0 )
Then 2161 storm(i,j) = .false.
2162 cat_crt(i,j) = .false.
2170 if ( storm(i,j) .and. slp(i,j)<1000.0 )
then 2171 depress(i,j) = 1000. - slp(i,j)
2179 used = send_data(
idiag%id_c15, depress, time)
2180 if(
idiag%id_f15>0) used = send_data(
idiag%id_f15, tc_count, time)
2184 slon =
rad2deg*atm(n)%gridstruct%agrid(i,j,1)
2185 slat =
rad2deg*atm(n)%gridstruct%agrid(i,j,2)
2187 if ( slat>0. .and. slat<40. .and. slon>110. .and. slon<180. )
then 2188 depress(i,j) = -depress(i,j)
2192 call prt_maxmin(
'Depress', depress, isc, iec, jsc, jec, 0, 1, 1.)
2195 if ( atm(n)%gridstruct%agrid(i,j,2)<0.)
then 2201 call prt_maxmin(
'NH Deps', depress, isc, iec, jsc, jec, 0, 1, 1.)
2206 tmp =
rad2deg * atm(n)%gridstruct%agrid(i,j,1)
2207 if ( tmp<280. )
then 2212 call prt_maxmin(
'ATL Deps', depress, isc, iec, jsc, jec, 0, 1, 1.)
2217 if(
idiag%id_c25>0)
then 2220 if ( cat_crt(i,j) .and. slp(i,j)<980.0 )
then 2221 depress(i,j) = 980. - slp(i,j)
2229 used = send_data(
idiag%id_c25, depress, time)
2230 if(
idiag%id_f25>0) used = send_data(
idiag%id_f25, tc_count, time)
2234 if(
idiag%id_c35>0)
then 2237 if ( cat_crt(i,j) .and. slp(i,j)<964.0 )
then 2238 depress(i,j) = 964. - slp(i,j)
2246 used = send_data(
idiag%id_c35, depress, time)
2247 if(
idiag%id_f35>0) used = send_data(
idiag%id_f35, tc_count, time)
2251 if(
idiag%id_c45>0)
then 2254 if ( cat_crt(i,j) .and. slp(i,j)<944.0 )
then 2255 depress(i,j) = 944. - slp(i,j)
2263 used = send_data(
idiag%id_c45, depress, time)
2264 if(
idiag%id_f45>0) used = send_data(
idiag%id_f45, tc_count, time)
2267 if (
idiag%id_c15>0)
then 2272 deallocate(tc_count)
2275 if(
idiag%id_slp>0 )
deallocate( slp )
2284 idg(:) =
idiag%id_t(:)
2286 do_cs_intp = .false.
2288 if ( idg(i)>0 )
then 2294 if ( do_cs_intp )
then 2295 if(.not.
allocated (a3) )
allocate( a3(isc:iec,jsc:jec,
nplev) )
2297 plevs, wz, atm(n)%peln, idg, a3, 1)
2299 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
2302 call prt_mxm(
'T100:', a3(isc:iec,jsc:jec,11), isc, iec, jsc, jec, 0, 1, 1., &
2303 atm(n)%gridstruct%area_64, atm(n)%domain)
2304 if (.not. atm(n)%gridstruct%bounded_domain)
then 2310 slat = atm(n)%gridstruct%agrid(i,j,2)*
rad2deg 2311 if( (slat>-10.0 .and. slat<10.) )
then 2312 sar = sar + atm(n)%gridstruct%area(i,j)
2313 tmp = tmp + a3(i,j,11)*atm(n)%gridstruct%area(i,j)
2317 call mp_reduce_sum(sar)
2318 call mp_reduce_sum(tmp)
2319 if ( sar > 0. )
then 2320 if (
master)
write(*,*)
'Tropical [10s,10n] mean T100 =', tmp/sar
2322 if (
master)
write(*,*)
'Warning: problem computing tropical mean T100' 2327 call prt_mxm(
'T200:', a3(isc:iec,jsc:jec,
k200), isc, iec, jsc, jec, 0, 1, 1., &
2328 atm(n)%gridstruct%area_64, atm(n)%domain)
2329 if (.not. atm(n)%gridstruct%bounded_domain)
then 2334 slat = atm(n)%gridstruct%agrid(i,j,2)*
rad2deg 2335 if( (slat>-20 .and. slat<20) )
then 2336 sar = sar + atm(n)%gridstruct%area(i,j)
2337 tmp = tmp + a3(i,j,
k200)*atm(n)%gridstruct%area(i,j)
2341 call mp_reduce_sum(sar)
2342 call mp_reduce_sum(tmp)
2343 if ( sar > 0. )
then 2344 if (
master)
write(*,*)
'Tropical [-20.,20.] mean T200 =', tmp/sar
2351 if (
idiag%id_t_plev>0)
then 2352 if(.not.
allocated (a3) )
allocate( a3(isc:iec,jsc:jec,
nplev) )
2355 plevs, wz, atm(n)%peln, id1, a3, 1)
2356 used=send_data(
idiag%id_t_plev, a3(isc:iec,jsc:jec,:), time)
2360 if(
idiag%id_mq > 0)
then 2365 a2(i,j) = -1.e-18 * atm(n)%ps(i,j)*
idiag%zxg(i,j)
2368 used = send_data(
idiag%id_mq, a2, time)
2370 tot_mq =
g_sum( atm(n)%domain, a2, isc, iec, jsc, jec, ngc, atm(n)%gridstruct%area_64, 0)
2373 if(
master)
write(*,*)
'Total (global) mountain torque (Hadleys)=', tot_mq
2377 if (
idiag%id_ts > 0) used = send_data(
idiag%id_ts, atm(n)%ts(isc:iec,jsc:jec), time)
2379 if (
idiag%id_tq>0 )
then 2380 nwater = atm(1)%flagstruct%nwat
2386 a2(i,j) = a2(i,j) + sum(atm(n)%q(i,j,k,1:nwater))*atm(n)%delp(i,j,k)
2390 used = send_data(
idiag%id_tq, a2*
ginv, time)
2393 cl = get_tracer_index(model_atmos,
'Cl')
2394 cl2 = get_tracer_index(model_atmos,
'Cl2')
2395 if (cl > 0 .and. cl2 > 0)
then 2396 allocate(var2(isc:iec,jsc:jec))
2401 var2(i,j) = var2(i,j) + atm(n)%delp(i,j,k)
2406 if (
idiag%id_acl > 0 )
then 2413 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,cl)*atm(n)%delp(i,j,k)
2420 a2(i,j) = a2(i,j) / var2(i,j)
2423 used = send_data(
idiag%id_acl, a2, time)
2425 if (
idiag%id_acl2 > 0 )
then 2432 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,cl2)*atm(n)%delp(i,j,k)
2439 a2(i,j) = a2(i,j) / var2(i,j)
2442 used = send_data(
idiag%id_acl2, a2, time)
2444 if (
idiag%id_acly > 0 )
then 2452 mm = (atm(n)%q(i,j,k,cl)+2.*atm(n)%q(i,j,k,cl2))*atm(n)%delp(i,j,k)
2453 a2(i,j) = a2(i,j) + mm
2454 qm = qm + mm*atm(n)%gridstruct%area_64(i,j)
2461 a2(i,j) = a2(i,j) / var2(i,j)
2464 used = send_data(
idiag%id_acly, a2, time)
2467 e2 = e2 + ((a2(i,j) -
qcly0)**2)*atm(n)%gridstruct%area_64(i,j)
2468 einf = max(einf, abs(a2(i,j) -
qcly0))
2471 if (
prt_minmax .and. .not. atm(n)%gridstruct%bounded_domain)
then 2472 call mp_reduce_sum(qm)
2473 call mp_reduce_max(einf)
2474 call mp_reduce_sum(e2)
2476 write(*,*)
' TERMINATOR TEST: ' 2477 write(*,*)
' chlorine mass: ', qm/(4.*pi*radius*radius)
2478 write(*,*)
' L2 err: ', sqrt(e2)/sqrt(4.*pi*radius*radius)/
qcly0 2479 write(*,*)
' max err: ', einf/
qcly0 2488 if (
idiag%id_iw>0 )
then 2494 a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) * &
2504 a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) * &
2514 a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k) * &
2520 used = send_data(
idiag%id_iw, a2*
ginv, time)
2522 if (
idiag%id_lw>0 )
then 2528 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
liq_wat)*atm(n)%delp(i,j,k)
2537 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
rainwat)*atm(n)%delp(i,j,k)
2542 used = send_data(
idiag%id_lw, a2*
ginv, time)
2548 if (
idiag%id_intqv>0 )
then 2554 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
sphum)*atm(n)%delp(i,j,k)
2559 used = send_data(
idiag%id_intqv, a2*
ginv, time)
2561 if (
idiag%id_intql>0 )
then 2567 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
liq_wat)*atm(n)%delp(i,j,k)
2572 used = send_data(
idiag%id_intql, a2*
ginv, time)
2574 if (
idiag%id_intqi>0 )
then 2580 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
ice_wat)*atm(n)%delp(i,j,k)
2585 used = send_data(
idiag%id_intqi, a2*
ginv, time)
2587 if (
idiag%id_intqr>0 )
then 2593 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
rainwat)*atm(n)%delp(i,j,k)
2598 used = send_data(
idiag%id_intqr, a2*
ginv, time)
2600 if (
idiag%id_intqs>0 )
then 2606 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
snowwat)*atm(n)%delp(i,j,k)
2611 used = send_data(
idiag%id_intqs, a2*
ginv, time)
2613 if (
idiag%id_intqg>0 )
then 2619 a2(i,j) = a2(i,j) + atm(n)%q(i,j,k,
graupel)*atm(n)%delp(i,j,k)
2624 used = send_data(
idiag%id_intqg, a2*
ginv, time)
2628 if ( (
idiag%id_ctt>0 .or.
idiag%id_ctp>0 .or.
idiag%id_ctz>0).and. atm(n)%flagstruct%nwat==6)
then 2629 allocate ( var1(isc:iec,jsc:jec) )
2630 allocate ( var2(isc:iec,jsc:jec) )
2637 if( tmp>5.e-6 )
then 2638 a2(i,j) = atm(n)%pt(i,j,k)
2639 var1(i,j) = 0.01*atm(n)%pe(i,k,j)
2640 var2(i,j) = wz(i,j,k) - wz(i,j,npz+1)
2642 elseif( k==npz )
then 2652 if (
idiag%id_ctt>0 )
then 2653 used = send_data(
idiag%id_ctt, a2, time)
2656 if (
idiag%id_ctp>0 )
then 2657 used = send_data(
idiag%id_ctp, var1, time)
2661 if (
idiag%id_ctz>0 )
then 2662 used = send_data(
idiag%id_ctz, var2, time)
2683 wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,
liq_wat)*atm(n)%delp(i,j,k)
2693 wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,
ice_wat)*atm(n)%delp(i,j,k)
2698 if (
idiag%id_qn>0 ) used = send_data(
idiag%id_qn, wk, time)
2699 if (
idiag%id_qn200>0 )
then 2701 used=send_data(
idiag%id_qn200, a2, time)
2703 if (
idiag%id_qn500>0 )
then 2705 used=send_data(
idiag%id_qn500, a2, time)
2707 if (
idiag%id_qn850>0 )
then 2709 used=send_data(
idiag%id_qn850, a2, time)
2713 if (
idiag%id_qp>0 )
then 2727 wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,
rainwat)*atm(n)%delp(i,j,k)
2737 wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,
snowwat)*atm(n)%delp(i,j,k)
2747 wk(i,j,k) = wk(i,j,k) + atm(n)%q(i,j,k,
graupel)*atm(n)%delp(i,j,k)
2752 used = send_data(
idiag%id_qp, wk, time)
2755 if(
idiag%id_us > 0 .and.
idiag%id_vs > 0)
then 2756 u2(:,:) = atm(n)%ua(isc:iec,jsc:jec,npz)
2757 v2(:,:) = atm(n)%va(isc:iec,jsc:jec,npz)
2760 a2(i,j) = sqrt(u2(i,j)**2 + v2(i,j)**2)
2763 used=send_data(
idiag%id_us, u2, time)
2764 used=send_data(
idiag%id_vs, v2, time)
2768 if(
idiag%id_tb > 0)
then 2769 a2(:,:) = atm(n)%pt(isc:iec,jsc:jec,npz)
2770 used=send_data(
idiag%id_tb, a2, time)
2772 call prt_mxm(
'T_bot:', a2, isc, iec, jsc, jec, 0, 1, 1., atm(n)%gridstruct%area_64, atm(n)%domain)
2775 if(
idiag%id_ua > 0) used=send_data(
idiag%id_ua, atm(n)%ua(isc:iec,jsc:jec,:), time)
2776 if(
idiag%id_va > 0) used=send_data(
idiag%id_va, atm(n)%va(isc:iec,jsc:jec,:), time)
2780 allocate( a3(isc:iec,jsc:jec,npz) )
2785 wk(i,j,k) = atm(n)%w(i,j,k)*atm(n)%delp(i,j,k)*
ginv 2790 if (
idiag%id_uw > 0)
then 2794 a3(i,j,k) = atm(n)%ua(i,j,k)*wk(i,j,k)
2798 used = send_data(
idiag%id_uw, a3, time)
2800 if (
idiag%id_vw > 0)
then 2804 a3(i,j,k) = atm(n)%va(i,j,k)*wk(i,j,k)
2808 used = send_data(
idiag%id_vw, a3, time)
2811 if (
idiag%id_hw > 0)
then 2812 allocate(cvm(isc:iec))
2816 call moist_cv(isc,iec,isd,ied,jsd,jed,npz,j,k,atm(n)%flagstruct%nwat,
sphum,
liq_wat,
rainwat, &
2819 a3(i,j,k) = atm(n)%pt(i,j,k)*cvm(i)*wk(i,j,k)
2822 cv_vapor = cp_vapor - rvgas
2824 a3(i,j,k) = atm(n)%pt(i,j,k)*cv_vapor*wk(i,j,k)
2829 used = send_data(
idiag%id_hw, a3, time)
2833 if (
idiag%id_qvw > 0)
then 2837 a3(i,j,k) = atm(n)%q(i,j,k,
sphum)*wk(i,j,k)
2841 used = send_data(
idiag%id_qvw, a3, time)
2843 if (
idiag%id_qlw > 0)
then 2844 if (
liq_wat < 0 .or.
rainwat < 0)
call mpp_error(fatal,
'qlw does not work without liq_wat and rainwat defined')
2848 a3(i,j,k) = (atm(n)%q(i,j,k,
liq_wat)+atm(n)%q(i,j,k,
rainwat))*wk(i,j,k)
2852 used = send_data(
idiag%id_qlw, a3, time)
2854 if (
idiag%id_qiw > 0)
then 2856 call mpp_error(fatal,
'qiw does not work without ice_wat, snowwat, and graupel defined')
2861 a3(i,j,k) = (atm(n)%q(i,j,k,
ice_wat)+atm(n)%q(i,j,k,
snowwat)+atm(n)%q(i,j,k,
graupel))*wk(i,j,k)
2865 used = send_data(
idiag%id_qiw, a3, time)
2867 if (
idiag%id_o3w > 0)
then 2869 call mpp_error(fatal,
'o3w does not work without o3mr defined')
2874 a3(i,j,k) = atm(n)%q(i,j,k,
o3mr)*wk(i,j,k)
2878 used = send_data(
idiag%id_o3w, a3, time)
2884 if(
idiag%id_ke > 0)
then 2889 a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k)*(atm(n)%ua(i,j,k)**2+atm(n)%va(i,j,k)**2)
2896 a2(i,j) = 0.5*a2(i,j)/(atm(n)%ps(i,j)-
ptop)
2899 used=send_data(
idiag%id_ke, a2, time)
2901 tot_mq =
g_sum( atm(n)%domain, a2, isc, iec, jsc, jec, ngc, atm(n)%gridstruct%area_64, 1)
2902 if (
master)
write(*,*)
'SQRT(2.*KE; m/s)=', sqrt(2.*tot_mq)
2908 if(
idiag%id_delp > 0 .or.
idiag%id_cape > 0 .or.
idiag%id_cin > 0 .or. ((.not. atm(n)%flagstruct%hydrostatic) .and.
idiag%id_pfnh > 0))
then 2912 wk(i,j,k) = atm(n)%delp(i,j,k)*(1.-sum(atm(n)%q(i,j,k,2:atm(n)%flagstruct%nwat)))
2916 if (
idiag%id_delp > 0) used=send_data(
idiag%id_delp, wk, time)
2919 if( ( (.not. atm(n)%flagstruct%hydrostatic) .and.
idiag%id_pfnh > 0) .or.
idiag%id_cape > 0 .or.
idiag%id_cin > 0)
then 2924 wk(i,j,k) = -wk(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas* &
2927 wk(i,j,k) = -wk(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas* &
2928 atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,
sphum))
2936 used=send_data(
idiag%id_pfnh, wk, time)
2939 if(
idiag%id_delp > 0) used=send_data(
idiag%id_delp, atm(n)%delp(isc:iec,jsc:jec,:), time)
2941 if( (.not. atm(n)%flagstruct%hydrostatic) .and. (
idiag%id_pfnh > 0 .or.
idiag%id_cape > 0 .or.
idiag%id_cin > 0))
then 2946 wk(i,j,k) = -atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas* &
2949 wk(i,j,k) = -atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas* &
2950 atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,
sphum))
2955 used=send_data(
idiag%id_pfnh, wk, time)
2959 if( atm(n)%flagstruct%hydrostatic .and. (
idiag%id_pfhy > 0 .or.
idiag%id_cape > 0 .or.
idiag%id_cin > 0) )
then 2963 wk(i,j,k) = 0.5 *(atm(n)%pe(i,k,j)+atm(n)%pe(i,k+1,j))
2967 used=send_data(
idiag%id_pfhy, wk, time)
2970 if (
idiag%id_cape > 0 .or.
idiag%id_cin > 0)
then 2973 allocate(var2(isc:iec,jsc:jec))
2974 allocate(a3(isc:iec,jsc:jec,npz))
2976 call eqv_pot(a3, atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%peln, atm(n)%pkz, atm(n)%q(isd,jsd,1,
sphum), &
2977 isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%flagstruct%moist_phys)
2985 call getcape(npz, wk(i,j,:), atm(n)%pt(i,j,:), -atm(n)%delz(i,j,:), atm(n)%q(i,j,:,
sphum), a3(i,j,:), a2(i,j), var2(i,j), source_in=1)
2989 if (
idiag%id_cape > 0)
then 2991 call prt_maxmin(
' CAPE (J/kg)', a2, isc,iec,jsc,jec, 0, 1, 1.)
2993 used=send_data(
idiag%id_cape, a2, time)
2995 if (
idiag%id_cin > 0)
then 2997 call prt_maxmin(
' CIN (J/kg)', var2, isc,iec,jsc,jec, 0, 1, 1.)
2999 used=send_data(
idiag%id_cin, var2, time)
3008 if((.not. atm(n)%flagstruct%hydrostatic) .and.
idiag%id_delz > 0)
then 3012 wk(i,j,k) = -atm(n)%delz(i,j,k)
3016 used=send_data(
idiag%id_delz, wk, time)
3022 if (
idiag%id_pmask>0)
then 3025 a2(i,j) = exp((atm(n)%peln(i,npz+1,j)+atm(n)%peln(i,npz+1,j))*0.5)*0.01
3029 used=send_data(
idiag%id_pmask, a2, time)
3034 if (
idiag%id_pmaskv2>0)
then 3037 a2(i,j) = exp((atm(n)%peln(i,npz,j)+atm(n)%peln(i,npz+1,j))*0.5)*0.01
3040 used=send_data(
idiag%id_pmaskv2, a2, time)
3044 & .or.
idiag%id_w1km>0 .or.
idiag%id_basedbz>0 .or.
idiag%id_dbz4km>0)
then 3045 if (.not.
allocated(wz))
allocate ( wz(isc:iec,jsc:jec,npz+1) )
3046 if ( atm(n)%flagstruct%hydrostatic)
then 3056 wz(i,j,k) = wz(i,j,k+1) - (rdgas*rgrav)*atm(n)%pt(i,j,k)*(atm(n)%peln(i,k,j) - atm(n)%peln(i,k+1,j))
3069 wz(i,j,k) = wz(i,j,k+1) - atm(n)%delz(i,j,k)
3075 call prt_maxmin(
'ZTOP', wz(isc:iec,jsc:jec,1)+atm(n)%phis(isc:iec,jsc:jec)/grav, isc, iec, jsc, jec, 0, 1, 1.e-3)
3078 if (
idiag%id_rain5km>0 )
then 3079 rainwat = get_tracer_index(model_atmos,
'rainwat')
3080 call interpolate_z(isc, iec, jsc, jec, npz, 5.e3, wz, atm(n)%q(isc:iec,jsc:jec,:,
rainwat), a2)
3081 used=send_data(
idiag%id_rain5km, a2, time)
3084 if (
idiag%id_w5km>0 )
then 3085 call interpolate_z(isc, iec, jsc, jec, npz, 5.e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
3086 used=send_data(
idiag%id_w5km, a2, time)
3089 if (
idiag%id_w2500m>0 )
then 3090 call interpolate_z(isc, iec, jsc, jec, npz, 2.5e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
3091 used=send_data(
idiag%id_w2500m, a2, time)
3094 if (
idiag%id_w1km>0 )
then 3095 call interpolate_z(isc, iec, jsc, jec, npz, 1.e3, wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
3096 used=send_data(
idiag%id_w1km, a2, time)
3099 if (
idiag%id_w100m>0 )
then 3100 call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%w(isc:iec,jsc:jec,:), a2)
3101 used=send_data(
idiag%id_w100m, a2, time)
3104 if (
idiag%id_u100m>0 )
then 3105 call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%ua(isc:iec,jsc:jec,:), a2)
3106 used=send_data(
idiag%id_u100m, a2, time)
3109 if (
idiag%id_v100m>0 )
then 3110 call interpolate_z(isc, iec, jsc, jec, npz, 100., wz, atm(n)%va(isc:iec,jsc:jec,:), a2)
3111 used=send_data(
idiag%id_v100m, a2, time)
3116 & .or.
idiag%id_dbztop>0 .or.
idiag%id_dbz_m10C>0))
then 3118 if (.not.
allocated(a3))
allocate(a3(isc:iec,jsc:jec,npz))
3121 call dbzcalc(atm(n)%q, atm(n)%pt, atm(n)%delp, atm(n)%peln, atm(n)%delz, &
3122 a3, a2, allmax, atm(n)%bd, npz, atm(n)%ncnst, atm(n)%flagstruct%hydrostatic, &
3123 zvir, .false., .false., .false., .true. )
3125 if (
idiag%id_dbz > 0) used=send_data(
idiag%id_dbz, a3, time)
3126 if (
idiag%id_maxdbz > 0) used=send_data(
idiag%id_maxdbz, a2, time)
3128 if (
idiag%id_basedbz > 0)
then 3130 call cs_interpolator(isc, iec, jsc, jec, npz, a3, 1000., wz, a2, -20.)
3131 used=send_data(
idiag%id_basedbz, a2, time)
3135 if (
idiag%id_dbz4km > 0)
then 3137 call cs_interpolator(isc, iec, jsc, jec, npz, a3, 4000., wz, a2, -20.)
3138 used=send_data(
idiag%id_dbz4km, a2, time)
3140 if (
idiag%id_dbztop > 0)
then 3145 if (wz(i,j,k) >= 25000. )
continue 3146 if (a3(i,j,k) >= 18.5 )
then 3153 used=send_data(
idiag%id_dbztop, a2, time)
3155 if (
idiag%id_dbz_m10C > 0)
then 3160 if (wz(i,j,k) >= 25000. )
exit 3161 if (atm(n)%pt(i,j,k) <= 263.14 )
then 3168 used=send_data(
idiag%id_dbz_m10C, a2, time)
3172 call mpp_max(allmax)
3173 if (
master)
write(*,*)
'max reflectivity = ', allmax,
' dBZ' 3182 if(.not.
allocated(a3))
allocate( a3(isc:iec,jsc:jec,
nplev) )
3184 idg(:) =
idiag%id_u(:)
3186 do_cs_intp = .false.
3188 if ( idg(i)>0 )
then 3194 if ( do_cs_intp )
then 3196 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
3199 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
3203 if (
idiag%id_u_plev>0)
then 3206 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
3207 used=send_data(
idiag%id_u_plev, a3(isc:iec,jsc:jec,:), time)
3211 idg(:) =
idiag%id_v(:)
3213 do_cs_intp = .false.
3215 if ( idg(i)>0 )
then 3221 if ( do_cs_intp )
then 3223 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
3226 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
3230 if (
idiag%id_v_plev>0)
then 3233 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
3234 used=send_data(
idiag%id_v_plev, a3(isc:iec,jsc:jec,:), time)
3238 idg(:) =
idiag%id_q(:)
3240 do_cs_intp = .false.
3242 if ( idg(i)>0 )
then 3248 if ( do_cs_intp )
then 3249 call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%q(isc:iec,jsc:jec,:,
sphum),
nplev, &
3250 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, 0)
3253 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
3257 if (
idiag%id_q_plev>0)
then 3259 call cs3_interpolator(isc,iec,jsc,jec,npz, atm(n)%q(isc:iec,jsc:jec,:,
sphum),
nplev, &
3260 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, 0)
3261 used=send_data(
idiag%id_q_plev, a3(isc:iec,jsc:jec,:), time)
3265 idg(:) =
idiag%id_omg(:)
3267 do_cs_intp = .false.
3269 if ( idg(i)>0 )
then 3274 if ( do_cs_intp )
then 3276 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), idg, a3, -1)
3279 if (idg(i)>0) used=send_data(idg(i), a3(isc:iec,jsc:jec,i), time)
3283 if (
idiag%id_omg_plev>0)
then 3286 pout, wz, atm(n)%pe(isc:iec,1:npz+1,jsc:jec), id1, a3, -1)
3287 used=send_data(
idiag%id_omg_plev, a3(isc:iec,jsc:jec,:), time)
3290 if(
allocated(a3) )
deallocate (a3)
3293 if (
idiag%id_sl12>0 )
then 3296 a2(i,j) = sqrt(atm(n)%ua(i,j,12)**2 + atm(n)%va(i,j,12)**2)
3299 used=send_data(
idiag%id_sl12, a2, time)
3301 if (
idiag%id_sl13>0 )
then 3304 a2(i,j) = sqrt(atm(n)%ua(i,j,13)**2 + atm(n)%va(i,j,13)**2)
3307 used=send_data(
idiag%id_sl13, a2, time)
3310 if ( (.not.atm(n)%flagstruct%hydrostatic) .and.
idiag%id_w200>0 )
then 3312 200.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
3313 used=send_data(
idiag%id_w200, a2, time)
3316 if ( (.not.atm(n)%flagstruct%hydrostatic) .and.
idiag%id_w500>0 )
then 3318 500.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
3319 used=send_data(
idiag%id_w500, a2, time)
3321 if ( (.not.atm(n)%flagstruct%hydrostatic) .and.
idiag%id_w700>0 )
then 3323 700.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
3324 used=send_data(
idiag%id_w700, a2, time)
3326 if ( (.not.atm(n)%flagstruct%hydrostatic) .and.
idiag%id_w850>0 .or.
idiag%id_x850>0)
then 3328 850.e2, atm(n)%peln, atm(n)%w(isc:iec,jsc:jec,:), a2)
3329 used=send_data(
idiag%id_w850, a2, time)
3331 if (
idiag%id_x850>0 .and.
idiag%id_vort850>0 )
then 3332 x850(:,:) = x850(:,:)*a2(:,:)
3333 used=send_data(
idiag%id_x850, x850, time)
3339 if ( .not.atm(n)%flagstruct%hydrostatic .and.
idiag%id_w>0 )
then 3340 used=send_data(
idiag%id_w, atm(n)%w(isc:iec,jsc:jec,:), time)
3342 if ( .not. atm(n)%flagstruct%hydrostatic .and. (
idiag%id_wmaxup>0 .or.
idiag%id_wmaxdn>0) )
then 3343 allocate(var2(isc:iec,jsc:jec))
3349 if (atm(n)%pe(i,k,j) <= 400.e2)
continue 3350 a2(i,j) = max(a2(i,j),atm(n)%w(i,j,k))
3351 var2(i,j) = min(var2(i,j),atm(n)%w(i,j,k))
3355 if (
idiag%id_wmaxup > 0)
then 3356 used=send_data(
idiag%id_wmaxup, a2, time)
3358 if (
idiag%id_wmaxdn > 0)
then 3359 used=send_data(
idiag%id_wmaxdn, var2, time)
3364 if(
idiag%id_pt > 0) used=send_data(
idiag%id_pt , atm(n)%pt (isc:iec,jsc:jec,:), time)
3365 if(
idiag%id_omga > 0) used=send_data(
idiag%id_omga, atm(n)%omga(isc:iec,jsc:jec,:), time)
3366 if(
idiag%id_diss > 0) used=send_data(
idiag%id_diss, atm(n)%diss_est(isc:iec,jsc:jec,:), time)
3368 allocate( a3(isc:iec,jsc:jec,npz) )
3369 if(
idiag%id_theta_e > 0 )
then 3371 if ( atm(n)%flagstruct%adiabatic .and. atm(n)%flagstruct%kord_tm>0 )
then 3375 a3(i,j,k) = atm(n)%pt(i,j,k)
3380 call eqv_pot(a3, atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%peln, atm(n)%pkz, atm(n)%q(isd,jsd,1,
sphum), &
3381 isc, iec, jsc, jec, ngc, npz, atm(n)%flagstruct%hydrostatic, atm(n)%flagstruct%moist_phys)
3384 if (
idiag%id_theta_e > 0)
then 3386 used=send_data(
idiag%id_theta_e, a3, time)
3388 theta_d = get_tracer_index(model_atmos,
'theta_d')
3389 if ( theta_d>0 )
then 3396 a2(i,j) = a2(i,j) + atm(n)%delp(i,j,k)*(atm(n)%q(i,j,k,theta_d)-a3(i,j,k))**2
3400 call prt_mxm(
'PT_SUM', a2, isc, iec, jsc, jec, 0, 1, 1.e-5, atm(n)%gridstruct%area_64, atm(n)%domain)
3405 a3(i,j,k) = atm(n)%q(i,j,k,theta_d)/a3(i,j,k) - 1.
3409 call prt_maxmin(
'Theta_Err (%)', a3, isc, iec, jsc, jec, 0, npz, 100.)
3416 if(
idiag%id_ppt> 0)
then 3418 allocate (
idiag%pt1(npz) )
3419 if( .not.
allocated(a3) )
allocate ( a3(isc:iec,jsc:jec,npz) )
3421 call gw_1d(npz, 1000.e2, atm(n)%ak, atm(n)%ak, atm(n)%ak(1), 10.e3,
idiag%pt1)
3425 if (.not. atm(n)%flagstruct%hydrostatic)
then 3429 wk(i,j,k) = (atm(n)%pt(i,j,k)*exp(-kappa*log(-atm(n)%delp(i,j,k)/(atm(n)%delz(i,j,k)*grav)*rdgas* &
3433 atm(n)%pt(i,j,k)*(1.+zvir*atm(n)%q(i,j,k,
sphum)))) -
idiag%pt1(k)) *
pk0 3445 wk(i,j,k) = (atm(n)%pt(i,j,k)/atm(n)%pkz(i,j,k) -
idiag%pt1(k)) *
pk0 3450 used=send_data(
idiag%id_ppt, wk, time)
3453 call prt_maxmin(
'PoTemp', wk, isc, iec, jsc, jec, 0, npz, 1.)
3456 if(
allocated(a3) )
deallocate ( a3 )
3457 deallocate (
idiag%pt1 )
3462 do itrac=1, atm(n)%ncnst
3463 call get_tracer_names (model_atmos, itrac, tname)
3464 if (
idiag%id_tracer(itrac) > 0 .and. itrac.gt.nq)
then 3465 used = send_data(
idiag%id_tracer(itrac), atm(n)%qdiag(isc:iec,jsc:jec,:,itrac), time )
3467 used = send_data(
idiag%id_tracer(itrac), atm(n)%q(isc:iec,jsc:jec,:,itrac), time )
3469 if (itrac .le. nq)
then 3471 isc, iec, jsc, jec, ngc, npz, 1.)
3474 isc, iec, jsc, jec, ngc, npz, 1.)
3483 if (
idiag%id_tracer_dmmr(itrac) > 0 .or.
idiag%id_tracer_dvmr(itrac) > 0)
then 3484 if (itrac .gt. nq)
then 3485 dmmr(:,:,:) = atm(n)%qdiag(isc:iec,jsc:jec,1:npz,itrac) &
3486 /(1.0-atm(n)%q(isc:iec,jsc:jec,1:npz,1))
3488 dmmr(:,:,:) = atm(n)%q(isc:iec,jsc:jec,1:npz,itrac) &
3489 /(1.0-atm(n)%q(isc:iec,jsc:jec,1:npz,1))
3491 dvmr(:,:,:) = dmmr(isc:iec,jsc:jec,1:npz) * wtmair/
idiag%w_mr(itrac)
3492 used = send_data(
idiag%id_tracer_dmmr(itrac), dmmr, time )
3493 used = send_data(
idiag%id_tracer_dvmr(itrac), dvmr, time )
3495 call prt_maxmin(trim(tname)//
'_dmmr', dmmr, &
3496 isc, iec, jsc, jec, 0, npz, 1.)
3497 call prt_maxmin(trim(tname)//
'_dvmr', dvmr, &
3498 isc, iec, jsc, jec, 0, npz, 1.)
3504 if ( .not. atm(n)%gridstruct%bounded_domain )
then 3510 a2(i,j) = max(a2(i,j), atm(n)%q(i,j,k,
cld_amt) )
3514 call prt_gb_nh_sh(
'Max_cld GB_NH_SH_EQ',isc,iec, jsc,jec, a2, atm(n)%gridstruct%area_64(isc:iec,jsc:jec), &
3515 atm(n)%gridstruct%agrid_64(isc:iec,jsc:jec,2))
3522 call debug_column(atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%u, atm(n)%v, atm(n)%w, atm(n)%q, &
3523 atm(n)%npz, atm(n)%ncnst,
sphum, atm(n)%flagstruct%nwat, atm(n)%flagstruct%hydrostatic, atm(n)%bd, time)
3527 call sounding_column(atm(n)%pt, atm(n)%delp, atm(n)%delz, atm(n)%u, atm(n)%v, atm(n)%q, atm(n)%peln, atm(n)%pkz, atm(n)%phis, &
3528 atm(n)%npz, atm(n)%ncnst,
sphum, atm(n)%flagstruct%nwat, atm(n)%flagstruct%hydrostatic, atm(n)%flagstruct%moist_phys, &
3529 zvir, atm(n)%ng, atm(n)%bd, time)
3540 if (
allocated(a3))
deallocate(a3)
3541 if (
allocated(wz))
deallocate(wz)
3542 if (
allocated(dmmr))
deallocate(dmmr)
3543 if (
allocated(dvmr))
deallocate(dvmr)
3545 call nullify_domain()
3550 subroutine wind_max(isc, iec, jsc, jec ,isd, ied, jsd, jed, us, vs, ws_max, domain)
3551 integer isc, iec, jsc, jec
3552 integer isd, ied, jsd, jed
3553 real,
intent(in),
dimension(isc:iec,jsc:jec):: us, vs
3554 real,
intent(out) :: ws_max(isc:iec,jsc:jec)
3555 type(domain2d),
intent(INOUT) :: domain
3557 real :: wx(isc:iec,jsd:jed), ws(isd:ied,jsd:jed)
3563 ws(i,j) = sqrt(us(i,j)**2 + vs(i,j)**2)
3567 call mpp_update_domains( ws, domain )
3571 wx(i,j) = max(ws(i-3,j), ws(i-2,j), ws(i-1,j), ws(i,j), ws(i+1,j), ws(i+2,j), ws(i+3,j))
3577 ws_max(i,j) = max(wx(i,j-3), wx(i,j-2), wx(i,j-1), wx(i,j), wx(i,j+1), wx(i,j+2), wx(i,j+3))
3584 subroutine get_vorticity(isc, iec, jsc, jec ,isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
3585 integer isd, ied, jsd, jed, npz
3586 integer isc, iec, jsc, jec
3587 real,
intent(in) :: u(isd:ied, jsd:jed+1, npz), v(isd:ied+1, jsd:jed, npz)
3588 real,
intent(out) :: vort(isc:iec, jsc:jec, npz)
3589 real,
intent(IN) :: dx(isd:ied,jsd:jed+1)
3590 real,
intent(IN) :: dy(isd:ied+1,jsd:jed)
3591 real,
intent(IN) :: rarea(isd:ied,jsd:jed)
3593 real :: utmp(isc:iec, jsc:jec+1), vtmp(isc:iec+1, jsc:jec)
3599 utmp(i,j) = u(i,j,k)*dx(i,j)
3604 vtmp(i,j) = v(i,j,k)*dy(i,j)
3610 vort(i,j,k) = rarea(i,j)*(utmp(i,j)-utmp(i,j+1)-vtmp(i,j)+vtmp(i+1,j))
3618 subroutine get_height_field(is, ie, js, je, ng, km, hydrostatic, delz, wz, pt, q, peln, zvir)
3619 integer,
intent(in):: is, ie, js, je, km, ng
3620 real,
intent(in):: peln(is:ie,km+1,js:je)
3621 real,
intent(in):: pt(is-ng:ie+ng,js-ng:je+ng,km)
3622 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
3623 real,
intent(in):: delz(is:,js:,1:)
3624 real,
intent(in):: zvir
3625 logical,
intent(in):: hydrostatic
3626 real,
intent(out):: wz(is:ie,js:je,km+1)
3635 wz(i,j,km+1) =
idiag%zsurf(i,j)
3637 if (hydrostatic )
then 3641 wz(i,j,k) = wz(i,j,k+1) + gg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas)) &
3642 *(peln(i,k+1,j)-peln(i,k,j))
3644 wz(i,j,k) = wz(i,j,k+1) + gg*pt(i,j,k)*(1.+zvir*q(i,j,k,
sphum)) &
3645 *(peln(i,k+1,j)-peln(i,k,j))
3652 wz(i,j,k) = wz(i,j,k+1) - delz(i,j,k)
3660 subroutine range_check_3d(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_range, Time)
3661 character(len=*),
intent(in):: qname
3662 integer,
intent(in):: is, ie, js, je
3663 integer,
intent(in):: n_g, km
3664 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3665 real,
intent(in):: pos(is-n_g:ie+n_g, js-n_g:je+n_g,2)
3666 real,
intent(in):: q_low, q_hi
3667 logical,
optional,
intent(out):: bad_range
3668 type(time_type),
optional,
intent(IN) :: Time
3672 integer year, month, day, hour, minute, second
3674 if (
present(bad_range) ) bad_range = .false.
3681 if( q(i,j,k) < qmin )
then 3683 elseif( q(i,j,k) > qmax )
then 3690 call mp_reduce_min(qmin)
3691 call mp_reduce_max(qmax)
3693 if( qmin<q_low .or. qmax>q_hi )
then 3694 if(
master)
write(*,*)
'Range_check Warning:', qname,
' max = ', qmax,
' min = ', qmin
3695 if (
present(time))
then 3696 call get_date(time, year, month, day, hour, minute, second)
3697 if (
master)
write(*,999) year, month, day, hour, minute, second
3698 999
format(
' Range violation on: ', i4,
'/', i02,
'/', i02,
' ', i02,
':', i02,
':', i02)
3700 if (
present(bad_range) )
then 3705 if (
present(bad_range) )
then 3707 if ( bad_range .EQV. .false. )
return 3711 if( q(i,j,k)<q_low .or. q(i,j,k)>q_hi )
then 3712 write(*,998) k,i,j, pos(i,j,1)*
rad2deg, pos(i,j,2)*
rad2deg, qname, q(i,j,k)
3714 998
format(
'Warn_K=',i4,
' (i,j)=',2i5,
' (lon,lat)=',f7.3,1x,f7.3,1x, a,
' =',f10.5)
3715 997
format(
' K=',i4,3x,f10.5)
3716 if ( k/= 1 )
write(*,997) k-1, q(i,j,k-1)
3717 if ( k/=km )
write(*,997) k+1, q(i,j,k+1)
3722 call mpp_error(note,
'==> Error from range_check: data out of bound')
3727 subroutine range_check_2d(qname, q, is, ie, js, je, n_g, pos, q_low, q_hi, bad_range, Time)
3728 character(len=*),
intent(in):: qname
3729 integer,
intent(in):: is, ie, js, je
3730 integer,
intent(in):: n_g
3731 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g)
3732 real,
intent(in):: pos(is-n_g:ie+n_g, js-n_g:je+n_g,2)
3733 real,
intent(in):: q_low, q_hi
3734 logical,
optional,
intent(out):: bad_range
3735 type(time_type),
optional,
intent(IN) :: Time
3739 integer year, month, day, hour, minute, second
3741 if (
present(bad_range) ) bad_range = .false.
3747 if( q(i,j) < qmin )
then 3749 elseif( q(i,j) > qmax )
then 3755 call mp_reduce_min(qmin)
3756 call mp_reduce_max(qmax)
3758 if( qmin<q_low .or. qmax>q_hi )
then 3759 if(
master)
write(*,*)
'Range_check Warning:', qname,
' max = ', qmax,
' min = ', qmin
3760 if (
present(time))
then 3761 call get_date(time, year, month, day, hour, minute, second)
3762 if (
master)
write(*,999) year, month, day, hour, minute, second
3763 999
format(
' Range violation on: ', i4,
'/', i02,
'/', i02,
' ', i02,
':', i02,
':', i02)
3765 if (
present(bad_range) )
then 3770 if (
present(bad_range) )
then 3772 if ( bad_range .EQV. .false. )
return 3775 if( q(i,j)<q_low .or. q(i,j)>q_hi )
then 3776 write(*,*)
'Warn_(i,j)=',i,j, pos(i,j,1)*
rad2deg, pos(i,j,2)*
rad2deg, q(i,j)
3780 call mpp_error(note,
'==> Error from range_check: data out of bound')
3785 subroutine prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
3786 character(len=*),
intent(in):: qname
3787 integer,
intent(in):: is, ie, js, je
3788 integer,
intent(in):: n_g, km
3789 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3790 real,
intent(in):: fac
3795 master = (mpp_pe()==mpp_root_pe()) .or. is_master()
3805 if( q(i,j,k) < qmin )
then 3807 elseif( q(i,j,k) > qmax )
then 3814 call mp_reduce_min(qmin)
3815 call mp_reduce_max(qmax)
3818 write(*,*) qname//trim(
gn),
' max = ', qmax*fac,
' min = ', qmin*fac
3823 subroutine prt_mxm(qname, q, is, ie, js, je, n_g, km, fac, area, domain)
3824 character(len=*),
intent(in):: qname
3825 integer,
intent(in):: is, ie, js, je
3826 integer,
intent(in):: n_g, km
3827 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3828 real,
intent(in):: fac
3831 real(kind=R_GRID),
intent(IN):: area(is-3:ie+3, js-3:je+3)
3832 type(domain2d),
intent(INOUT) :: domain
3834 real qmin, qmax, gmean
3838 master = (mpp_pe()==mpp_root_pe()) .or. is_master()
3848 if( q(i,j,k) < qmin )
then 3850 elseif( q(i,j,k) > qmax )
then 3857 call mp_reduce_min(qmin)
3858 call mp_reduce_max(qmax)
3862 gmean =
g_sum(domain, q(is:ie,js:je,km), is, ie, js, je, 3, area, 1)
3864 if(
master)
write(6,*) qname,
gn, qmax*fac, qmin*fac, gmean*fac
3869 subroutine prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain)
3871 integer,
intent(in):: is, ie, js, je
3872 integer,
intent(in):: nq, n_g, km, nwat
3873 real,
intent(in):: ps(is-n_g:ie+n_g, js-n_g:je+n_g)
3874 real,
intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3875 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km, nq)
3876 real(kind=R_GRID),
intent(IN):: area(is-n_g:ie+n_g,js-n_g:je+n_g)
3877 type(domain2d),
intent(INOUT) :: domain
3879 real psq(is:ie,js:je,nwat), psqv(is:ie,js:je)
3880 real q_strat(is:ie,js:je)
3881 real qtot(nwat), qwat
3882 real psmo, totw, psdry
3883 integer k, n, kstrat
3886 sphum = get_tracer_index(model_atmos,
'sphum')
3887 liq_wat = get_tracer_index(model_atmos,
'liq_wat')
3888 ice_wat = get_tracer_index(model_atmos,
'ice_wat')
3890 rainwat = get_tracer_index(model_atmos,
'rainwat')
3891 snowwat = get_tracer_index(model_atmos,
'snowwat')
3892 graupel = get_tracer_index(model_atmos,
'graupel')
3895 psmo =
g_sum(domain, ps(is:ie,js:je), is, ie, js, je, n_g, area, 1)
3896 if(
master )
write(*,*)
'Total surface pressure (mb)', trim(
gn),
' = ', 0.01*psmo
3897 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,1 ), psqv(is,js))
3902 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
sphum ), psq(is,js,
sphum ))
3905 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
liq_wat), psq(is,js,
liq_wat))
3908 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
rainwat), psq(is,js,
rainwat))
3912 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
ice_wat), psq(is,js,
ice_wat))
3915 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
snowwat), psq(is,js,
snowwat))
3917 call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,
graupel), psq(is,js,
graupel))
3921 if (
idiag%phalf(2)< 75. )
then 3924 if (
idiag%phalf(k+1) > 75. )
exit 3927 call z_sum(is, ie, js, je, kstrat, n_g, delp, q(is-n_g,js-n_g,1,
sphum), q_strat(is,js))
3928 psmo =
g_sum(domain, q_strat(is,js), is, ie, js, je, n_g, area, 1) * 1.e6 &
3929 /
p_sum(is, ie, js, je, kstrat, n_g, delp, area, domain)
3930 if(
master)
write(*,*)
'Mean specific humidity (mg/kg) above 75 mb', trim(
gn),
'=', psmo
3937 psmo =
g_sum(domain, ps(is:ie,js:je), is, ie, js, je, n_g, area, 1)
3940 qtot(n) =
g_sum(domain, psq(is,js,n), is, ie, js, je, n_g, area, 1)
3943 totw = sum(qtot(1:nwat))
3947 write(*,*)
'Total surface pressure (mb)', trim(
gn),
' = ', 0.01*psmo
3948 write(*,*)
'mean dry surface pressure', trim(
gn),
' = ', 0.01*psdry
3949 write(*,*)
'Total Water Vapor (kg/m**2)', trim(
gn),
' =', qtot(
sphum)*
ginv 3951 write(*,*)
'--- Micro Phys water substances (kg/m**2) ---' 3952 write(*,*)
'Total cloud water', trim(
gn),
'=', qtot(
liq_wat)*
ginv 3954 write(*,*)
'Total rain water', trim(
gn),
'=', qtot(
rainwat)*
ginv 3956 write(*,*)
'Total cloud ice ', trim(
gn),
'=', qtot(
ice_wat)*
ginv 3960 write(*,*)
'Total graupel ', trim(
gn),
'=', qtot(
graupel)*
ginv 3961 write(*,*)
'---------------------------------------------' 3962 elseif ( nwat==2 )
then 3963 write(*,*)
'GFS condensate (kg/m^2)', trim(
gn),
'=', qtot(
liq_wat)*
ginv 3970 subroutine z_sum(is, ie, js, je, km, n_g, delp, q, sum2)
3971 integer,
intent(in):: is, ie, js, je, n_g, km
3972 real,
intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3973 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3974 real,
intent(out):: sum2(is:ie,js:je)
3980 sum2(i,j) = delp(i,j,1)*q(i,j,1)
3984 sum2(i,j) = sum2(i,j) + delp(i,j,k)*q(i,j,k)
3989 end subroutine z_sum 3992 real function p_sum(is, ie, js, je, km, n_g, delp, area, domain)
3993 integer,
intent(in):: is, ie, js, je, n_g, km
3994 real,
intent(in):: delp(is-n_g:ie+n_g, js-n_g:je+n_g, km)
3995 real(kind=R_GRID),
intent(IN) :: area(is-n_g:ie+n_g, js-n_g:je+n_g)
3996 real :: sum2(is:ie,js:je)
3998 type(domain2d),
intent(INOUT) :: domain
4003 sum2(i,j) = delp(i,j,1)
4007 sum2(i,j) = sum2(i,j) + delp(i,j,k)
4011 p_sum =
g_sum(domain, sum2, is, ie, js, je, n_g, area, 1)
4020 integer,
intent(in):: is, ie, js, je, km, ng
4021 integer,
intent(in):: kd
4022 real,
intent(in):: wz(is:ie,js:je,km+1)
4023 real,
intent(in):: ts(is-ng:ie+ng,js-ng:je+ng)
4024 real,
intent(in):: peln(is:ie,km+1,js:je)
4025 real,
intent(in):: height(kd)
4026 real,
intent(out):: a2(is:ie,js:je,kd)
4027 real,
optional,
intent(in):: fac
4042 if ( height(n) >= wz(i,j,km+1) )
then 4047 if( height(n) < wz(i,j,k) .and. height(n) >= wz(i,j,k+1) )
then 4049 ptmp = peln(i,k,j) + (peln(i,k+1,j)-peln(i,k,j)) * &
4050 (wz(i,j,k)-height(n)) / (wz(i,j,k)-wz(i,j,k+1))
4051 a2(i,j,n) = exp(ptmp)
4060 tm = rdgas*
ginv*(ts(i,j) + 3.25e-3*(wz(i,j,km)-height(n)))
4061 a2(i,j,n) = exp( peln(i,km+1,j) + (wz(i,j,km+1) - height(n))/tm )
4063 500
if (
present(fac) ) a2(i,j,n) = fac * a2(i,j,n)
4071 subroutine get_height_given_pressure(is, ie, js, je, km, wz, kd, id, log_p, peln, a2)
4072 integer,
intent(in):: is, ie, js, je, km
4073 integer,
intent(in):: kd
4074 integer,
intent(in):: id(kd)
4075 real,
intent(in):: log_p(kd)
4077 real,
intent(in):: wz(is:ie,js:je,km+1)
4078 real,
intent(in):: peln(is:ie,km+1,js:je)
4079 real,
intent(out):: a2(is:ie,js:je,kd)
4081 real,
dimension(2*km+1):: pn, gz
4082 integer n,i,j,k, k1, k2, l
4084 k2 = max(12, km/2+1)
4099 gz(k) = 2.*gz(km+1) - gz(l)
4100 pn(k) = 2.*pn(km+1) - pn(l)
4104 if( id(n)<0 )
goto 1000
4106 if( log_p(n) <= pn(k+1) .and. log_p(n) >= pn(k) )
then 4107 a2(i,j,n) = gz(k) + (gz(k+1)-gz(k))*(log_p(n)-pn(k))/(pn(k+1)-pn(k))
4118 subroutine prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
4119 character(len=*),
intent(in):: qname
4120 integer,
intent(in):: is, ie, js, je, ng, km
4121 real,
intent(in):: press
4122 real,
intent(in):: peln(is:ie,km+1,js:je)
4123 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
4124 real,
intent(in):: delz(is:,js:,1:)
4125 real(kind=R_GRID),
intent(in),
dimension(is:ie, js:je):: area, lat
4127 real:: a2(is:ie,js:je)
4128 real(kind=R_GRID),
dimension(2*km+1):: pn, gz
4129 real(kind=R_GRID):: log_p
4130 integer i,j,k, k2, l
4133 k2 = max(12, km/2+1)
4145 gz(km+1) = phis(i,j)/grav
4147 gz(k) = gz(k+1) - delz(i,j,k)
4151 gz(k) = 2.*gz(km+1) - gz(l)
4152 pn(k) = 2.*pn(km+1) - pn(l)
4156 if( log_p <= pn(k+1) .and. log_p >= pn(k) )
then 4157 a2(i,j) = gz(k) + (gz(k+1)-gz(k))*(log_p-pn(k))/(pn(k+1)-pn(k))
4167 subroutine prt_gb_nh_sh(qname, is,ie, js,je, a2, area, lat)
4168 character(len=*),
intent(in):: qname
4169 integer,
intent(in):: is, ie, js, je
4170 real,
intent(in),
dimension(is:ie, js:je):: a2
4171 real(kind=R_GRID),
intent(in),
dimension(is:ie, js:je):: area, lat
4173 real(R_GRID),
parameter:: rad2deg = 180./pi
4175 real:: t_eq, t_nh, t_sh, t_gb
4176 real:: area_eq, area_nh, area_sh, area_gb
4179 t_eq = 0. ; t_nh = 0.; t_sh = 0.; t_gb = 0.
4180 area_eq = 0.; area_nh = 0.; area_sh = 0.; area_gb = 0.
4183 slat = lat(i,j)*rad2deg
4184 area_gb = area_gb + area(i,j)
4185 t_gb = t_gb + a2(i,j)*area(i,j)
4186 if( (slat>-20. .and. slat<20.) )
then 4187 area_eq = area_eq + area(i,j)
4188 t_eq = t_eq + a2(i,j)*area(i,j)
4189 elseif( slat>=20. .and. slat<80. )
then 4190 area_nh = area_nh + area(i,j)
4191 t_nh = t_nh + a2(i,j)*area(i,j)
4192 elseif( slat<=-20. .and. slat>-80. )
then 4193 area_sh = area_sh + area(i,j)
4194 t_sh = t_sh + a2(i,j)*area(i,j)
4198 call mp_reduce_sum(area_gb)
4199 call mp_reduce_sum( t_gb)
4200 call mp_reduce_sum(area_nh)
4201 call mp_reduce_sum( t_nh)
4202 call mp_reduce_sum(area_sh)
4203 call mp_reduce_sum( t_sh)
4204 call mp_reduce_sum(area_eq)
4205 call mp_reduce_sum( t_eq)
4207 if (area_gb <= 1.) area_gb = -1.0
4208 if (area_nh <= 1.) area_nh = -1.0
4209 if (area_sh <= 1.) area_sh = -1.0
4210 if (area_eq <= 1.) area_eq = -1.0
4211 if (is_master())
write(*,*) qname, t_gb/area_gb, t_nh/area_nh, t_sh/area_sh, t_eq/area_eq
4215 subroutine cs3_interpolator(is, ie, js, je, km, qin, kd, pout, wz, pe, id, qout, iv)
4219 integer,
intent(in):: is, ie, js, je, km, iv
4220 integer,
intent(in):: kd
4221 integer,
intent(in):: id(kd)
4222 real,
intent(in):: pout(kd)
4223 real,
intent(in):: pe(is:ie,km+1,js:je)
4224 real,
intent(in):: qin(is:ie,js:je,km)
4225 real,
intent(in):: wz(is:ie,js:je,km+1)
4226 real,
intent(out):: qout(is:ie,js:je,kd)
4228 real,
parameter:: gcp = grav / cp_air
4229 real:: qe(is:ie,km+1)
4230 real,
dimension(is:ie,km):: q2, dp
4232 integer:: i,j,k, n, k1
4240 dp(i,k) = pe(i,k+1,j) - pe(i,k,j)
4241 q2(i,k) = qin(i,j,k)
4245 call cs_prof(q2, dp, qe, km, is, ie, iv)
4250 if ( id(n) > 0 )
then 4251 if( pout(n) <= pe(i,1,j) )
then 4253 qout(i,j,n) = qe(i,1)
4254 elseif ( pout(n) >= pe(i,km+1,j) )
then 4260 qout(i,j,n) = gcp*exp(kappa*pout(n)) * (wz(i,j,km-2) - wz(i,j,km)) / &
4261 ( exp(kappa*pe(i,km,j)) - exp(kappa*pe(i,km-2,j)) )
4263 qout(i,j,n) = qe(i,km+1)
4267 if ( pout(n)>=pe(i,k,j) .and. pout(n) <= pe(i,k+1,j) )
then 4269 a6 = 3.*(2.*q2(i,k) - (qe(i,k)+qe(i,k+1)))
4270 s0 = (pout(n)-pe(i,k,j)) / dp(i,k)
4271 qout(i,j,n) = qe(i,k) + s0*(qe(i,k+1)-qe(i,k)+a6*(1.-s0))
4277 500
if ( iv==0 ) qout(i,j,n) = max(0., qout(i,j,n))
4287 subroutine cs_interpolator(is, ie, js, je, km, qin, zout, wz, qout, qmin)
4288 integer,
intent(in):: is, ie, js, je, km
4289 real,
intent(in):: zout, qmin
4290 real,
intent(in):: qin(is:ie,js:je,km)
4291 real,
intent(in):: wz(is:ie,js:je,km+1)
4292 real,
intent(out):: qout(is:ie,js:je)
4294 real:: qe(is:ie,km+1)
4295 real,
dimension(is:ie,km):: q2, dz
4305 dz(i,k) = wz(i,j,k) - wz(i,j,k+1)
4306 q2(i,k) = qin(i,j,k)
4310 call cs_prof(q2, dz, qe, km, is, ie, 1)
4313 if( zout >= wz(i,j,1) )
then 4316 elseif ( zout <= wz(i,j,km+1) )
then 4317 qout(i,j) = qe(i,km+1)
4320 if ( zout<=wz(i,j,k) .and. zout >= wz(i,j,k+1) )
then 4322 a6 = 3.*(2.*q2(i,k) - (qe(i,k)+qe(i,k+1)))
4323 s0 = (wz(i,j,k)-zout) / dz(i,k)
4324 qout(i,j) = qe(i,k) + s0*(qe(i,k+1)-qe(i,k)+a6*(1.-s0))
4329 500 qout(i,j) = max(qmin, qout(i,j))
4337 subroutine cs_prof(q2, delp, q, km, i1, i2, iv)
4339 integer,
intent(in):: i1, i2, km
4340 integer,
intent(in):: iv
4341 real,
intent(in) :: q2(i1:i2,km)
4342 real,
intent(in) :: delp(i1:i2,km)
4343 real,
intent(out):: q(i1:i2,km+1)
4347 real bet, a_bot, grat
4351 grat = delp(i,2) / delp(i,1)
4352 bet = grat*(grat+0.5)
4353 q(i,1) = ( (grat+grat)*(grat+1.)*q2(i,1) + q2(i,2) ) / bet
4354 gam(i,1) = ( 1. + grat*(grat+1.5) ) / bet
4359 d4(i) = delp(i,k-1) / delp(i,k)
4360 bet = 2. + d4(i) + d4(i) - gam(i,k-1)
4361 q(i,k) = ( 3.*(q2(i,k-1)+d4(i)*q2(i,k)) - q(i,k-1) )/bet
4362 gam(i,k) = d4(i) / bet
4367 a_bot = 1. + d4(i)*(d4(i)+1.5)
4368 q(i,km+1) = (2.*d4(i)*(d4(i)+1.)*q2(i,km)+q2(i,km-1)-a_bot*q(i,km)) &
4369 / ( d4(i)*(d4(i)+0.5) - a_bot*gam(i,km) )
4374 q(i,k) = q(i,k) - gam(i,k)*q(i,k+1)
4380 q(i,2) = min( q(i,2), max(q2(i,1), q2(i,2)) )
4381 q(i,2) = max( q(i,2), min(q2(i,1), q2(i,2)) )
4386 gam(i,k) = q2(i,k) - q2(i,k-1)
4393 if ( gam(i,k-1)*gam(i,k+1)>0. )
then 4395 q(i,k) = min( q(i,k), max(q2(i,k-1),q2(i,k)) )
4396 q(i,k) = max( q(i,k), min(q2(i,k-1),q2(i,k)) )
4398 if ( gam(i,k-1) > 0. )
then 4400 q(i,k) = max(q(i,k), min(q2(i,k-1),q2(i,k)))
4403 q(i,k) = min(q(i,k), max(q2(i,k-1),q2(i,k)))
4404 if ( iv==0 ) q(i,k) = max(0., q(i,k))
4412 q(i,km) = min( q(i,km), max(q2(i,km-1), q2(i,km)) )
4413 q(i,km) = max( q(i,km), min(q2(i,km-1), q2(i,km)) )
4421 integer,
intent(in):: is, ie, js, je, km
4422 real,
intent(in):: peln(is:ie,km+1,js:je)
4423 real,
intent(in):: a3(is:ie,js:je,km)
4424 real,
intent(in):: plev
4425 real,
intent(out):: a2(is:ie,js:je)
4439 pm(k) = 0.5*(peln(i,k,j)+peln(i,k+1,j))
4442 if( logp <= pm(1) )
then 4444 elseif ( logp >= pm(km) )
then 4445 a2(i,j) = a3(i,j,km)
4448 if( logp <= pm(k+1) .and. logp >= pm(k) )
then 4449 a2(i,j) = a3(i,j,k) + (a3(i,j,k+1)-a3(i,j,k))*(logp-pm(k))/(pm(k+1)-pm(k))
4460 subroutine interpolate_z(is, ie, js, je, km, zl, hght, a3, a2)
4462 integer,
intent(in):: is, ie, js, je, km
4463 real,
intent(in):: hght(is:ie,js:je,km+1)
4464 real,
intent(in):: a3(is:ie,js:je,km)
4465 real,
intent(in):: zl
4466 real,
intent(out):: a2(is:ie,js:je)
4476 zm(k) = 0.5*(hght(i,j,k)+hght(i,j,k+1))
4478 if( zl >= zm(1) )
then 4480 elseif ( zl <= zm(km) )
then 4481 a2(i,j) = a3(i,j,km)
4484 if( zl <= zm(k) .and. zl >= zm(k+1) )
then 4485 a2(i,j) = a3(i,j,k) + (a3(i,j,k+1)-a3(i,j,k))*(zm(k)-zl)/(zm(k)-zm(k+1))
4496 ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
4498 integer,
intent(in):: is, ie, js, je, ng, km, sphum
4499 real,
intent(in):: grav, zvir, z_bot, z_top
4500 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
4501 real,
intent(in):: delz(is:ie,js:je,km)
4502 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
4503 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
4504 real,
intent(in):: peln(is:ie,km+1,js:je)
4505 logical,
intent(in):: hydrostatic
4506 real,
intent(out):: srh(is:ie,js:je)
4517 real,
dimension(is:ie):: zh, uc, vc, dz, zh0
4518 integer i, j, k, k0, k1
4519 logical below(is:ie)
4541 if ( hydrostatic )
then 4543 dz(i) = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
4545 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
4548 dz(i) = - delz(i,j,k)
4550 zh(i) = zh(i) + dz(i)
4551 if (zh(i) <= z_bot )
continue 4552 if (zh(i) > z_bot .and. below(i))
then 4553 uc(i) = ua(i,j,k)*dz(i)
4554 vc(i) = va(i,j,k)*dz(i)
4555 zh0(i) = zh(i) - dz(i)
4559 elseif ( zh(i) < z_top )
then 4560 uc(i) = uc(i) + ua(i,j,k)*dz(i)
4561 vc(i) = vc(i) + va(i,j,k)*dz(i)
4564 uc(i) = uc(i) / (zh(i)-dz(i) - zh0(i))
4565 vc(i) = vc(i) / (zh(i)-dz(i) - zh0(i))
4573 srh(i,j) = 0.5*(va(i,j,k1)-vc(i))*(ua(i,j,k1-1)-ua(i,j,k1)) - &
4574 0.5*(ua(i,j,k1)-uc(i))*(va(i,j,k1-1)-va(i,j,k1))
4576 srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i))*(ua(i,j,k-1)-ua(i,j,k+1)) - &
4577 0.5*(ua(i,j,k)-uc(i))*(va(i,j,k-1)-va(i,j,k+1))
4585 subroutine helicity_relative_caps(is, ie, js, je, ng, km, zvir, sphum, srh, uc, vc, &
4586 ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
4588 integer,
intent(in):: is, ie, js, je, ng, km, sphum
4589 real,
intent(in):: grav, zvir, z_bot, z_top
4590 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
4591 real,
intent(in):: delz(is:ie,js:je,km)
4592 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
4593 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
4594 real,
intent(in):: peln(is:ie,km+1,js:je)
4595 real,
intent(in):: uc(is:ie,js:je), vc(is:ie,js:je)
4596 logical,
intent(in):: hydrostatic
4597 real,
intent(out):: srh(is:ie,js:je)
4607 real,
dimension(is:ie):: zh, dz, zh0
4608 integer i, j, k, k0, k1, n
4628 if ( hydrostatic )
then 4630 dz(i) = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
4632 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
4635 dz(i) = -delz(i,j,k)
4638 zh(i) = zh(i) + dz(i)
4639 if (zh(i) <= z_bot )
continue 4640 if (zh(i) > z_bot .and. below)
then 4641 zh0(i) = zh(i) - dz(i)
4645 elseif ( zh(i) < z_top )
then 4656 srh(i,j) = 0.5*(va(i,j,k1)-vc(i,j))*(ua(i,j,k1-1)-ua(i,j,k1)) - &
4657 0.5*(ua(i,j,k1)-uc(i,j))*(va(i,j,k1-1)-va(i,j,k1))
4659 srh(i,j) = srh(i,j) + 0.5*(va(i,j,k)-vc(i,j))*(ua(i,j,k-1)-ua(i,j,k+1)) - &
4660 0.5*(ua(i,j,k)-uc(i,j))*(va(i,j,k-1)-va(i,j,k+1))
4668 subroutine bunkers_vector(is, ie, js, je, ng, km, zvir, sphum, uc, vc, &
4669 ua, va, delz, q, hydrostatic, pt, peln, phis, grav)
4671 integer,
intent(in):: is, ie, js, je, ng, km, sphum
4672 real,
intent(in):: grav, zvir
4673 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, ua, va
4674 real,
intent(in):: delz(is:ie,js:je,km)
4675 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
4676 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
4677 real,
intent(in):: peln(is:ie,km+1,js:je)
4678 logical,
intent(in):: hydrostatic
4679 real,
intent(out):: uc(is:ie,js:je), vc(is:ie,js:je)
4682 real :: zh, dz, usfc, vsfc, u6km, v6km, umn, vmn
4683 real :: ushr, vshr, shrmag
4685 real,
parameter :: bunkers_d = 7.5
4686 logical :: has_sfc, has_6km
4711 if ( hydrostatic )
then 4713 dz = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
4715 dz = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
4726 umn = umn + ua(i,j,k)*dz
4727 vmn = vmn + va(i,j,k)*dz
4735 u6km = u6km + (ua(i,j,k) - u6km) / dz * (6000. - (zh - dz))
4736 v6km = v6km + (va(i,j,k) - v6km) / dz * (6000. - (zh - dz))
4738 umn = umn / (zh - dz)
4739 vmn = vmn / (zh - dz)
4743 shrmag = sqrt(ushr * ushr + vshr * vshr)
4744 uc(i,j) = umn + bunkers_d * vshr / shrmag
4745 vc(i,j) = vmn - bunkers_d * ushr / shrmag
4754 w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
4756 integer,
intent(in):: is, ie, js, je, ng, km, sphum
4757 real,
intent(in):: grav, zvir, z_bot, z_top
4758 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, w
4759 real,
intent(in),
dimension(is:ie,js:je,km):: vort
4760 real,
intent(in):: delz(is:ie,js:je,km)
4761 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
4762 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
4763 real,
intent(in):: peln(is:ie,km+1,js:je)
4764 logical,
intent(in):: hydrostatic
4765 real,
intent(out):: uh(is:ie,js:je)
4770 real,
dimension(is:ie):: zh, dz, zh0
4772 logical below(is:ie)
4792 if ( hydrostatic )
then 4794 dz(i) = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
4796 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
4799 dz(i) = - delz(i,j,k)
4801 zh(i) = zh(i) + dz(i)
4802 if (zh(i) <= z_bot )
continue 4803 if (zh(i) > z_bot .and. below(i))
then 4804 uh(i,j) = vort(i,j,k)*w(i,j,k)*(zh(i) - z_bot)
4807 elseif ( zh(i) < z_top )
then 4808 uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*dz(i)
4810 uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*(z_top - (zh(i)-dz(i)) )
4823 subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
4826 integer,
intent(in):: is, ie, js, je, ng, km
4827 real,
intent(in):: grav
4828 real,
intent(in):: pt(is-ng:ie+ng,js-ng:je+ng,km)
4829 real,
intent(in):: pkz(is:ie,js:je,km)
4830 real,
intent(in):: delp(is-ng:ie+ng,js-ng:je+ng,km)
4831 real,
intent(in):: f_d(is-ng:ie+ng,js-ng:je+ng)
4834 real,
intent(inout):: vort(is:ie,js:je,km)
4855 real w3d(is:ie,js:je,km)
4856 real te(is:ie,js:je,km+1), t2(is:ie,km), delp2(is:ie,km)
4857 real te2(is:ie,km+1)
4864 vort(i,j,1) = grav * (vort(i,j,1)+f_d(i,j)) / delp(i,j,1)
4874 t2(i,k) = pt(i,j,k) / pkz(i,j,k)
4875 w3d(i,j,k) = t2(i,k)
4876 delp2(i,k) = delp(i,j,k)
4880 call ppme(t2, te2, delp2, ie-is+1, km)
4884 te(i,j,k) = te2(i,k)
4894 vort(i,j,k) = (vort(i,j,k)+f_d(i,j)) * ( te(i,j,k)-te(i,j,k+1) ) &
4895 / ( w3d(i,j,k)*delp(i,j,k) ) * grav
4904 subroutine ppme(p,qe,delp,im,km)
4906 integer,
intent(in):: im, km
4907 real,
intent(in):: p(im,km)
4908 real,
intent(in):: delp(im,km)
4909 real,
intent(out)::qe(im,km+1)
4912 real dc(im,km),delq(im,km), a6(im,km)
4913 real c1, c2, c3, tmp, qmax, qmin
4914 real a1, a2, s1, s2, s3, s4, ss3, s32, s34, s42
4915 real a3, b2, sc, dm, d1, d2, f1, f2, f3, f4
4923 500 a6(i,k) = delp(i,k-1) + delp(i,k)
4927 delq(i,k) = p(i,k+1) - p(i,k)
4932 c1 = (delp(i,k-1)+0.5*delp(i,k))/a6(i,k+1)
4933 c2 = (delp(i,k+1)+0.5*delp(i,k))/a6(i,k)
4934 tmp = delp(i,k)*(c1*delq(i,k) + c2*delq(i,k-1)) / &
4935 (a6(i,k)+delp(i,k+1))
4936 qmax = max(p(i,k-1),p(i,k),p(i,k+1)) - p(i,k)
4937 qmin = p(i,k) - min(p(i,k-1),p(i,k),p(i,k+1))
4938 dc(i,k) = sign(min(abs(tmp),qmax,qmin), tmp)
4947 c1 = delq(i,k-1)*delp(i,k-1) / a6(i,k)
4948 a1 = a6(i,k-1) / (a6(i,k) + delp(i,k-1))
4949 a2 = a6(i,k+1) / (a6(i,k) + delp(i,k))
4950 qe(i,k) = p(i,k-1) + c1 + 2./(a6(i,k-1)+a6(i,k+1)) * &
4951 ( delp(i,k)*(c1*(a1 - a2)+a2*dc(i,k-1)) - &
4952 delp(i,k-1)*a1*dc(i,k ) )
4966 s3 = delp(i,2) + delp(i,3)
4973 a3 = (delq(i,2) - delq(i,1)*s3/s2) / (s3*ss3)
4975 if(abs(a3) .gt. 1.e-14)
then 4976 b2 = delq(i,1)/s2 - a3*(s1+s2)
4978 if(sc .lt. 0. .or. sc .gt. s1)
then 4979 qe(i,1) = p(i,1) - s1*(a3*s1 + b2)
4981 qe(i,1) = p(i,1) - delq(i,1)*s1/s2
4985 qe(i,1) = p(i,1) - delq(i,1)*s1/s2
4987 dc(i,1) = p(i,1) - qe(i,1)
4989 dm = delp(i,1) / (s34*ss3*(delp(i,2)+s3)*(s4+delp(i,1)))
4990 f1 = delp(i,2)*s34 / ( s2*ss3*(s4+delp(i,1)) )
4991 f2 = (delp(i,2)+s3) * (ss3*(delp(i,2)*s3+s34+delp(i,2)*s4) &
4992 + s42*(delp(i,2)+s3+s32/s2))
4993 f3 = -delp(i,2)*( ss3*(s32*(s3+s4)/(s4-delp(i,2)) &
4994 + (delp(i,2)*s3+s34+delp(i,2)*s4)) &
4995 + s42*(delp(i,2)+s3) )
4996 f4 = ss3*delp(i,2)*s32*(delp(i,2)+s3) / (s4-delp(i,2))
4997 qe(i,2) = f1*p(i,1)+(f2*p(i,2)+f3*p(i,3)+f4*p(i,4))*dm
5005 qm = (d2*p(i,km)+d1*p(i,km1)) / (d1+d2)
5006 dq = 2.*(p(i,km1)-p(i,km)) / (d1+d2)
5007 c1 = (qe(i,km1)-qm-d2*dq) / (d2*(2.*d2*d2+d1*(d2+3.*d1)))
5008 c3 = dq - 2.0*c1*(d2*(5.*d1+d2)-3.*d1**2)
5009 qe(i,km ) = qm - c1*d1*d2*(d2+3.*d1)
5010 qe(i,km+1) = d1*(8.*c1*d1**2-c3) + qe(i,km)
5015 subroutine rh_calc (pfull, t, qv, rh, do_cmip)
5017 real,
intent (in),
dimension(:,:) :: pfull, t, qv
5018 real,
intent (out),
dimension(:,:) :: rh
5019 real,
dimension(size(t,1),size(t,2)) :: esat
5020 logical,
intent(in),
optional :: do_cmip
5022 real,
parameter :: d622 = rdgas/rvgas
5023 real,
parameter :: d378 = 1.-d622
5025 logical :: do_simple = .false.
5030 call lookup_es (t, esat)
5031 rh(:,:) = pfull(:,:)
5032 rh(:,:) = max(rh(:,:),esat(:,:))
5033 rh(:,:)=100.*qv(:,:)/(d622*esat(:,:)/rh(:,:))
5035 if (
present(do_cmip))
then 5036 call compute_qs (t, pfull, rh, q=qv, es_over_liq_and_ice = .true.)
5037 rh(:,:)=100.*qv(:,:)/rh(:,:)
5039 call compute_qs (t, pfull, rh, q=qv)
5040 rh(:,:)=100.*qv(:,:)/rh(:,:)
5046 #ifdef SIMPLIFIED_THETA_E 5052 subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, qi, is, ie, js, je, ng, npz, &
5055 subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, &
5058 integer,
intent(in):: is,ie,js,je,ng,npz
5060 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz,*):: qi
5062 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz):: pt, delp, q
5063 real,
intent(in),
dimension(is: ,js: ,1: ):: delz
5064 real,
intent(in),
dimension(is:ie,npz+1,js:je):: peln
5065 real,
intent(in):: pkz(is:ie,js:je,npz)
5066 logical,
intent(in):: hydrostatic, moist
5068 real,
dimension(is:ie,js:je,npz),
intent(out) :: theta_e
5070 real,
parameter:: tice = 273.16
5071 real,
parameter:: c_liq = 4190.
5073 real,
parameter:: dc_vap = 0.
5075 real,
parameter:: dc_vap = cp_vapor - c_liq
5077 real(kind=R_GRID),
dimension(is:ie):: pd, rq
5078 real(kind=R_GRID) :: wfac
5090 q(i,j,k) = qi(i,j,k,1)
5101 if ( hydrostatic )
then 5103 rq(i) = max(0., wfac*q(i,j,k))
5104 pd(i) = (1.-rq(i))*delp(i,j,k) / (peln(i,k+1,j) - peln(i,k,j))
5109 rq(i) = max(0., wfac*q(i,j,k))
5110 pd(i) = -rdgas*pt(i,j,k)*(1.-rq(i))*delp(i,j,k)/(grav*delz(i,j,k))
5116 rq(i) = max(0., q(i,j,k))
5123 theta_e(i,j,k) = pt(i,j,k)*exp(kappa * (
virqd(qi(i,j,k,:))/
vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i))) * &
5124 exp(rq(i)*hlv/(cp_air*
vicpqd(qi(i,j,k,:))*pt(i,j,k)))
5126 theta_e(i,j,k) = pt(i,j,k)*exp(kappa*log(1.e5/pd(i))) * &
5127 exp(rq(i)*hlv/(cp_air*pt(i,j,k)))
5131 theta_e(i,j,k) = pt(i,j,k)*exp( rq(i)/(cp_air*
vicpqd(q(i,j,k,:))*pt(i,j,k))*(hlv+dc_vap*(pt(i,j,k)-tice)) &
5132 + rdgas*
virqd(qi(i,j,k,:)) / (cp_air*
vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i)) )
5134 theta_e(i,j,k) = pt(i,j,k)*exp( rq(i)/(cp_air*pt(i,j,k))*(hlv+dc_vap*(pt(i,j,k)-tice)) &
5135 + kappa*log(1.e5/pd(i)) )
5140 if ( hydrostatic )
then 5142 theta_e(i,j,k) = pt(i,j,k)*
pk0/pkz(i,j,k)
5148 theta_e(i,j,k) = pt(i,j,k)*exp( kappa * (
virqd(qi(i,j,k,:))/
vicpqd(qi(i,j,k,:)))*log(1.e5/pd(i)) )
5150 theta_e(i,j,k) = pt(i,j,k)*exp( kappa*log(1.e5/pd(i)) )
5167 subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, qi, is, ie, js, je, ng, npz, &
5169 subroutine eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, &
5176 integer,
intent(in):: is,ie,js,je,ng,npz
5178 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz,*):: qi
5180 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz):: q
5182 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz):: pt, delp
5183 real,
intent(in),
dimension(is: ,js: ,1: ):: delz
5184 real,
intent(in),
dimension(is:ie,npz+1,js:je):: peln
5185 real,
intent(in):: pkz(is:ie,js:je,npz)
5186 logical,
intent(in):: hydrostatic, moist
5188 real,
dimension(is:ie,js:je,npz),
intent(out) :: theta_e
5191 real,
dimension(is-ng:ie+ng,js-ng:je+ng,npz):: q
5193 real,
parameter:: cv_vap = cp_vapor - rvgas
5194 real,
parameter:: cappa_b = 0.2854
5195 real(kind=R_GRID):: cv_air, cappa, zvir
5196 real(kind=R_GRID):: p_mb(is:ie)
5197 real(kind=R_GRID) :: r, e, t_l, rdg, capa
5201 q(:,:,:) = qi(:,:,:,1)
5203 cv_air = cp_air - rdgas
5206 zvir = rvgas/rdgas - 1.
5221 if ( hydrostatic )
then 5223 p_mb(i) = 0.01*delp(i,j,k) / (peln(i,k+1,j) - peln(i,k,j))
5228 p_mb(i) = 0.01*rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)*
virq(qi(i,j,k,1:
num_gas))
5230 p_mb(i) = 0.01*rdg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)*(1.+zvir*q(i,j,k))
5237 cappa = rdgas/(rdgas+((1.-q(i,j,k))*cv_air+q(i,j,k)*cv_vap)/
virq(qi(i,j,k,1:
num_gas)))
5239 cappa = rdgas/(rdgas+((1.-q(i,j,k))*cv_air+q(i,j,k)*cv_vap)/(1.+zvir*q(i,j,k)))
5242 r = q(i,j,k)/(1.-q(i,j,k))*1000.
5245 e = p_mb(i)*r/(622.+r)
5248 t_l = 2840./(3.5*log(pt(i,j,k))-log(e)-4.805)+55.
5252 capa = cappa*(1. - r*0.28e-3)
5253 theta_e(i,j,k) = exp( (3.376/t_l-0.00254)*r*(1.+r*0.81e-3) )*pt(i,j,k)*(1000./p_mb(i))**capa
5256 if ( hydrostatic )
then 5258 theta_e(i,j,k) = pt(i,j,k)*
pk0/pkz(i,j,k)
5263 theta_e(i,j,k) = pt(i,j,k)*exp( kappa *
virqd(qi(i,j,k,1:
num_gas))/
vicpqd(qi(i,j,k,1:
num_gas)) *log(1000./p_mb(i)) )
5265 theta_e(i,j,k) = pt(i,j,k)*exp( kappa*log(1000./p_mb(i)) )
5280 w, delz, pt, delp, q, hs, area, domain, &
5281 sphum, liq_wat, rainwat, ice_wat, &
5282 snowwat, graupel, nwat, ua, va, moist_phys, te)
5287 integer,
intent(in):: km, is, ie, js, je, isd, ied, jsd, jed
5288 integer,
intent(in):: nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
5289 real,
intent(in),
dimension(isd:ied,jsd:jed,km):: ua, va, pt, delp, w
5290 real,
intent(in),
dimension(is:ie,js:je,km) :: delz
5291 real,
intent(in),
dimension(isd:ied,jsd:jed,km,nwat):: q
5292 real,
intent(in):: hs(isd:ied,jsd:jed)
5293 real,
intent(in):: area(isd:ied, jsd:jed)
5294 logical,
intent(in):: moist_phys
5295 type(domain2d),
intent(INOUT) :: domain
5296 real,
intent(out):: te(is:ie,js:je)
5298 real,
parameter:: c_liq = 4190.
5299 real(kind=R_Grid) :: area_l(isd:ied, jsd:jed)
5300 real,
parameter:: cv_vap = cp_vapor - rvgas
5301 real phiz(is:ie,km+1)
5302 real,
dimension(is:ie):: cvm, qc
5307 cv_air = cp_air - rdgas
5319 phiz(i,km+1) = hs(i,j)
5324 phiz(i,k) = phiz(i,k+1) - grav*delz(i,j,k)
5328 if ( moist_phys )
then 5330 call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
5331 ice_wat, snowwat, graupel, q, qc, cvm)
5333 te(i,j) = te(i,j) + delp(i,j,k)*( cvm(i)*pt(i,j,k) + hlv*q(i,j,k,sphum) + &
5334 0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
5341 te(i,j) = te(i,j) + delp(i,j,k)*( cv_air*
vicvqd(q(i,j,k,1:
num_gas))*pt(i,j,k) + &
5342 0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
5344 te(i,j) = te(i,j) + delp(i,j,k)*( cv_air*pt(i,j,k) + &
5345 0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )
5352 te(i,j) = te(i,j)/grav
5356 psm =
g_sum(domain, te, is, ie, js, je, 3, area_l, 1)
5357 if(
master )
write(*,*)
'TE ( Joule/m^2 * E9) =', psm * 1.e-9
5371 subroutine dbzcalc(q, pt, delp, peln, delz, &
5372 dbz, maxdbz, allmax, bd, npz, ncnst, &
5373 hydrostatic, zvir, in0r, in0s, in0g, iliqskin)
5417 integer,
intent(IN) :: npz, ncnst
5418 real,
intent(IN),
dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz) :: pt, delp
5419 real,
intent(IN),
dimension(bd%is:, bd%js:, 1:) :: delz
5420 real,
intent(IN),
dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst) :: q
5421 real,
intent(IN),
dimension(bd%is :bd%ie, npz+1, bd%js:bd%je) :: peln
5422 real,
intent(OUT),
dimension(bd%is :bd%ie, bd%js :bd%je , npz) :: dbz
5423 real,
intent(OUT),
dimension(bd%is :bd%ie, bd%js :bd%je) :: maxdbz
5424 logical,
intent(IN) :: hydrostatic, in0r, in0s, in0g, iliqskin
5425 real,
intent(IN) :: zvir
5426 real,
intent(OUT) :: allmax
5430 real(kind=R_GRID),
parameter:: vconr = 2503.23638966667
5431 real(kind=R_GRID),
parameter:: vcong = 87.2382675
5432 real(kind=R_GRID),
parameter:: vcons = 6.6280504
5433 real(kind=R_GRID),
parameter:: normr = 25132741228.7183
5434 real(kind=R_GRID),
parameter:: normg = 5026548245.74367
5435 real(kind=R_GRID),
parameter:: norms = 942477796.076938
5439 real,
parameter :: r1=1.e-15
5440 real,
parameter :: ron=8.e6
5441 real,
parameter :: ron2=1.e10
5442 real,
parameter :: son=2.e7
5443 real,
parameter :: gon=5.e7
5444 real,
parameter :: ron_min = 8.e6
5445 real,
parameter :: ron_qr0 = 0.00010
5446 real,
parameter :: ron_delqr0 = 0.25*ron_qr0
5447 real,
parameter :: ron_const1r = (ron2-ron_min)*0.5
5448 real,
parameter :: ron_const2r = (ron2+ron_min)*0.5
5449 real,
parameter :: rnzs = 3.0e6
5452 real,
parameter :: gamma_seven = 720.
5454 real,
parameter :: rhor = 1.0e3
5455 real,
parameter :: rhos = 100.
5456 real,
parameter :: rhog0 = 400.
5457 real,
parameter :: rhog = 500.
5459 real,
parameter :: alpha = 0.224
5460 real(kind=R_GRID),
parameter :: factor_s = gamma_seven * 1.e18 * (1./(pi*rhos))**1.75 &
5461 * (rhos/rhor)**2 * alpha
5462 real,
parameter :: qmin = 1.e-12
5463 real,
parameter :: tice = 273.16
5466 real(kind=R_GRID),
dimension(bd%is:bd%ie) :: rhoair, denfac, z_e
5467 real(kind=R_GRID):: qr1, qs1, qg1, t1, t2, t3, rwat, vtr, vtg, vts
5468 real(kind=R_GRID):: factorb_s, factorb_g
5469 real(kind=R_GRID):: temp_c, pres, sonv, gonv, ronv
5472 integer :: is, ie, js, je
5488 if (hydrostatic)
then 5491 rhoair(i) = delp(i,j,k)/( (peln(i,k+1,j)-peln(i,k,j)) * rdgas * pt(i,j,k) *
virq(q(i,j,k,1:
num_gas)) )
5493 rhoair(i) = delp(i,j,k)/( (peln(i,k+1,j)-peln(i,k,j)) * rdgas * pt(i,j,k) * ( 1. + zvir*q(i,j,k,
sphum) ) )
5495 denfac(i) = sqrt(min(10., 1.2/rhoair(i)))
5500 rhoair(i) = -delp(i,j,k)/(grav*delz(i,j,k))
5501 denfac(i) = sqrt(min(10., 1.2/rhoair(i)))
5512 t1 = rhoair(i)*max(qmin, q(i,j,k,
rainwat)+dim(q(i,j,k,
liq_wat), 1.0e-3))
5513 vtr = max(1.e-3, vconr*denfac(i)*exp(0.2 *log(t1/normr)))
5514 z_e(i) = 200.*exp(1.6*log(3.6e6*t1/rhor*vtr))
5519 t3 = rhoair(i)*max(qmin, q(i,j,k,
graupel))
5520 vtg = max(1.e-3, vcong*denfac(i)*exp(0.125 *log(t3/normg)))
5521 z_e(i) = z_e(i) + 200.*exp(1.6*log(3.6e6*t3/rhog*vtg))
5526 t2 = rhoair(i)*max(qmin, q(i,j,k,
snowwat))
5528 z_e(i) = z_e(i) + (factor_s/alpha)*t2*exp(0.75*log(t2/rnzs))
5532 dbz(i,j,k) = 10.*log10( max(0.01, z_e(i)) )
5541 maxdbz(i,j) = max(dbz(i,j,k), maxdbz(i,j))
5548 allmax = max(maxdbz(i,j), allmax)
5555 integer,
intent(in):: is, ie, js, je, km
5556 real,
intent(in),
dimension(is:ie,js:je,km):: vort
5557 real,
intent(inout),
dimension(is:ie,js:je):: maxvorthy1
5562 maxvorthy1(i,j)=max(maxvorthy1(i,j),vort(i,j,km))
5569 subroutine max_vorticity(is, ie, js, je, ng, km, zvir, sphum, delz, q, hydrostatic, &
5570 pt, peln, phis, grav, vort, maxvort, z_bot, z_top)
5571 integer,
intent(in):: is, ie, js, je, ng, km, sphum
5572 real,
intent(in):: grav, zvir, z_bot, z_top
5573 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt
5574 real,
intent(in),
dimension(is:ie,js:je,km):: vort
5575 real,
intent(in):: delz(is:ie,js:je,km)
5576 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
5577 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
5578 real,
intent(in):: peln(is:ie,km+1,js:je)
5579 logical,
intent(in):: hydrostatic
5581 real,
intent(inout),
dimension(is:ie,js:je):: maxvort
5584 real,
dimension(is:ie):: zh, dz, zh0
5585 integer i, j, k,klevel
5586 logical below(is:ie)
5598 if ( hydrostatic )
then 5600 dz(i) = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
5602 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
5605 dz(i) = - delz(i,j,k)
5607 zh(i) = zh(i) + dz(i)
5608 if (zh(i) <= z_bot )
continue 5609 if (zh(i) > z_bot .and. below(i))
then 5610 maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
5612 elseif ( zh(i) < z_top )
then 5613 maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
5615 maxvort(i,j) = max(maxvort(i,j),vort(i,j,k))
5626 subroutine max_uh(is, ie, js, je, ng, km, zvir, sphum, uphmax,uphmin, &
5627 w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
5629 integer,
intent(in):: is, ie, js, je, ng, km, sphum
5630 real,
intent(in):: grav, zvir, z_bot, z_top
5631 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, w
5632 real,
intent(in),
dimension(is:ie,js:je,km):: vort
5633 real,
intent(in):: delz(is:ie,js:je,km)
5634 real,
intent(in):: q(is-ng:ie+ng,js-ng:je+ng,km,*)
5635 real,
intent(in):: phis(is-ng:ie+ng,js-ng:je+ng)
5636 real,
intent(in):: peln(is:ie,km+1,js:je)
5637 logical,
intent(in):: hydrostatic
5638 real :: uh(is:ie,js:je)
5639 real,
intent(inout),
dimension(is:ie,js:je):: uphmax,uphmin
5644 real,
dimension(is:ie):: zh, dz, zh0
5646 logical below(is:ie)
5659 if ( hydrostatic )
then 5661 dz(i) = rdg*pt(i,j,k)*
virq(q(i,j,k,1:
num_gas))*(peln(i,k+1,j)-peln(i,k,j))
5663 dz(i) = rdg*pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))*(peln(i,k+1,j)-peln(i,k,j))
5666 dz(i) = - delz(i,j,k)
5668 zh(i) = zh(i) + dz(i)
5669 if (zh(i) <= z_bot )
continue 5670 if (zh(i) > z_bot .and. below(i))
then 5671 if(w(i,j,k).lt.0)
then 5675 uh(i,j) = vort(i,j,k)*w(i,j,k)*(zh(i) - z_bot)
5678 elseif ( zh(i) < z_top )
then 5679 if(w(i,j,k).lt.0)
then 5683 uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*dz(i)
5685 if(w(i,j,k).lt.0)
then 5689 uh(i,j) = uh(i,j) + vort(i,j,k)*w(i,j,k)*(z_top - (zh(i)-dz(i)) )
5693 if (uh(i,j) > uphmax(i,j))
then 5694 uphmax(i,j) = uh(i,j)
5695 elseif (uh(i,j) < uphmin(i,j))
then 5696 uphmin(i,j) = uh(i,j)
5702 subroutine max_vv(is,ie,js,je,npz,ng,up2,dn2,pe,w)
5704 integer,
intent(in):: is, ie, js, je, ng, npz
5706 real,
intent(in),
dimension(is-ng:ie+ng,js-ng:je+ng,npz):: w
5707 real,
intent(in):: pe(is-1:ie+1,npz+1,js-1:je+1)
5708 real,
intent(inout),
dimension(is:ie,js:je):: up2,dn2
5712 if (pe(i,k,j) >= 100.e2)
then 5713 up2(i,j) = max(up2(i,j),w(i,j,k))
5714 dn2(i,j) = min(dn2(i,j),w(i,j,k))
5724 if (atm%grid_Number > 1)
then 5725 write(
gn,
"(A2,I1)")
" g", atm%grid_number
5740 subroutine getcape( nk , p , t , dz, q, the, cape , cin, source_in )
5743 integer,
intent(in) :: nk
5744 real,
dimension(nk),
intent(in) :: p,t,dz,q,the
5745 real,
intent(out) :: cape,cin
5746 integer,
intent(IN),
OPTIONAL :: source_in
5790 real,
parameter :: pinc = 10000.0
5796 real,
parameter :: ml_depth = 200.0
5799 integer,
parameter :: adiabat = 1
5811 integer :: source = 1
5812 logical :: doit,ice,cloud,not_converged
5813 integer :: k,kmin,n,nloop,i,orec
5814 real,
dimension(nk) :: pi,th,thv,z,pt,pb,pc,pn,ptv
5816 real :: maxthe,parea,narea,lfc
5817 real :: th1,p1,t1,qv1,ql1,qi1,b1,pi1,thv1,qt,dp,frac
5818 real :: th2,p2,t2,qv2,ql2,qi2,b2,pi2,thv2
5819 real :: thlast,fliq,fice,tbar,qvbar,qlbar,qibar,lhv,lhs,lhf,rm,cpm
5820 real*8 :: avgth,avgqv
5824 real,
parameter :: g = 9.81
5825 real,
parameter :: p00 = 100000.0
5826 real,
parameter :: cp = 1005.7
5827 real,
parameter :: rd = 287.04
5828 real,
parameter :: rv = 461.5
5829 real,
parameter :: xlv = 2501000.0
5830 real,
parameter :: xls = 2836017.0
5831 real,
parameter :: t0 = 273.15
5832 real,
parameter :: cpv = 1875.0
5833 real,
parameter :: cpl = 4190.0
5834 real,
parameter :: cpi = 2118.636
5835 real,
parameter :: lv1 = xlv+(cpl-cpv)*t0
5836 real,
parameter :: lv2 = cpl-cpv
5837 real,
parameter :: ls1 = xls+(cpi-cpv)*t0
5838 real,
parameter :: ls2 = cpi-cpv
5840 real,
parameter :: rp00 = 1.0/p00
5841 real,
parameter :: eps = rd/rv
5842 real,
parameter :: reps = rv/rd
5843 real,
parameter :: rddcp = rd/cp
5844 real,
parameter :: cpdrd = cp/rd
5845 real,
parameter :: cpdg = cp/g
5847 real,
parameter :: converge = 0.1
5849 integer,
parameter :: debug_level = 0
5851 if (
present(source_in)) source = source_in
5858 pi(k) = (p(k)*rp00)**rddcp
5860 thv(k) = th(k)*(1.0+reps*q(k))/(1.0+q(k))
5867 z(k) = z(k+1) + 0.5*(dz(k+1)+dz(k))
5876 ELSEIF(source.eq.2)
THEN 5879 IF(p(1).lt.50000.0)
THEN 5887 if(p(k).ge.50000.0)
then 5888 if( the(nk).gt.maxthe )
then 5895 if(debug_level.ge.100) print *,
' kmin,maxthe = ',kmin,maxthe
5959 print *,
' Unknown value for source' 5961 print *,
' source = ',source
5963 call mpp_error(fatal,
" Unknown CAPE source")
5970 if( (source.eq.1).or.(source.eq.2) )
then 5979 elseif( source.eq.3 )
then 5983 thv2 = th2*(1.0+reps*qv2)/(1.0+qv2)
5987 b2 = g*( thv2-thv(kmin) )/thv(kmin)
6000 if(adiabat.eq.1.or.adiabat.eq.2)
then 6011 if(debug_level.ge.100)
then 6012 print *,
' Start loop:' 6013 print *,
' p2,th2,qv2 = ',p2,th2,qv2
6016 do while( doit .and. (k.gt.1) )
6023 if( dp.lt.pinc )
then 6026 nloop = 1 + int( dp/pinc )
6027 dp = dp/float(nloop)
6042 pi2 = (p2*rp00)**rddcp
6046 not_converged = .true.
6048 do while( not_converged )
6052 fliq = max(min((t2-233.15)/(273.15-233.15),1.0),0.0)
6058 qv2 = min( qt , fliq*
getqvs(p2,t2) + fice*
getqvi(p2,t2) )
6059 qi2 = max( fice*(qt-qv2) , 0.0 )
6060 ql2 = max( qt-qv2-qi2 , 0.0 )
6063 qvbar = 0.5*(qv1+qv2)
6064 qlbar = 0.5*(ql1+ql2)
6065 qibar = 0.5*(qi1+qi2)
6072 cpm=cp+cpv*qvbar+cpl*qlbar+cpi*qibar
6073 th2=th1*exp( lhv*(ql2-ql1)/(cpm*tbar) &
6074 +lhs*(qi2-qi1)/(cpm*tbar) &
6075 +(rm/cpm-rd/cp)*alog(p2/p1) )
6077 if(i.gt.90) print *,i,th2,thlast,th2-thlast
6079 print *,
' getcape() error: lack of convergence, stopping iteration' 6080 not_converged = .false.
6082 if( abs(th2-thlast).gt.converge )
then 6083 thlast=thlast+0.3*(th2-thlast)
6085 not_converged = .false.
6092 if( ql2.ge.1.0e-10 ) cloud = .true.
6094 IF(adiabat.eq.1.or.adiabat.eq.3)
THEN 6099 ELSEIF(adiabat.le.0.or.adiabat.ge.5)
THEN 6101 print *,
' Undefined adiabat' 6108 thv2 = th2*(1.0+reps*qv2)/(1.0+qv2+ql2+qi2)
6109 b2 = g*( thv2-thv(k) )/thv(k)
6115 if( (b2.ge.0.0) .and. (b1.lt.0.0) )
then 6119 parea = 0.5*b2*dz(k)*frac
6120 narea = narea-0.5*b1*dz(k)*(1.0-frac)
6121 if(debug_level.ge.200)
then 6122 print *,
' b1,b2 = ',b1,b2
6124 print *,
' frac = ',frac
6125 print *,
' parea = ',parea
6126 print *,
' narea = ',narea
6130 elseif( (b2.lt.0.0) .and. (b1.gt.0.0) )
then 6134 parea = 0.5*b1*dz(k)*frac
6135 narea = -0.5*b2*dz(k)*(1.0-frac)
6136 if(debug_level.ge.200)
then 6137 print *,
' b1,b2 = ',b1,b2
6139 print *,
' frac = ',frac
6140 print *,
' parea = ',parea
6141 print *,
' narea = ',narea
6143 elseif( b2.lt.0.0 )
then 6146 narea = narea-0.5*dz(k)*(b1+b2)
6149 parea = 0.5*dz(k)*(b1+b2)
6153 cape = cape + max(0.0,parea)
6155 if(debug_level.ge.200)
then 6156 write(6,102) p2,b1,b2,cape,cin,cloud
6157 102
format(5(f13.4),2x,l1)
6160 if( (p(k).le.10000.0).and.(b2.lt.0.0) )
then 6189 real function getqvs(p,t)
6194 real,
parameter :: eps = 287.04/461.5
6196 es = 611.2*exp(17.67*(t-273.15)/(t-29.65))
6204 real function getqvi(p,t)
6209 real,
parameter :: eps = 287.04/461.5
6211 es = 611.2*exp(21.8745584*(t-273.15)/(t-7.66))
6218 subroutine debug_column(pt, delp, delz, u, v, w, q, npz, ncnst, sphum, nwat, hydrostatic, bd, Time)
6221 integer,
intent(IN) :: npz, ncnst, sphum, nwat
6222 logical,
intent(IN) :: hydrostatic
6223 real,
dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),
intent(IN) :: pt, delp, w
6224 real,
dimension(bd%is:, bd%js:,1:),
intent(IN) :: delz
6225 real,
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz),
intent(IN) :: u
6226 real,
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz),
intent(IN) :: v
6227 real,
dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst),
intent(IN) :: q
6230 type(time_type),
intent(IN) :: Time
6231 integer :: i,j,k,n,l
6239 if (i < bd%is .or. i > bd%ie) cycle
6240 if (j < bd%js .or. j > bd%je) cycle
6246 write(
diag_debug_units(n),
'(A4, A7, A8, A6, A8, A8, A8, A8, A9)')
'k',
'T',
'delp',
'delz',
'u',
'v',
'w',
'sphum',
'cond' 6247 write(
diag_debug_units(n),
'(A4, A7, A8, A6, A8, A8, A8, A8, A9)')
'',
'K',
'mb',
'm',
'm/s',
'm/s',
'm/s',
'g/kg',
'g/kg' 6248 if (hydrostatic)
then 6249 call mpp_error(note,
'Hydrostatic debug sounding not yet supported')
6254 cond = cond + q(i,j,k,l)
6256 write(
diag_debug_units(n),
'(I4, F7.2, F8.3, I6, F8.3, F8.3, F8.3, F8.3, F9.5 )') &
6257 k, pt(i,j,k), delp(i,j,k)*0.01, -int(delz(i,j,k)), u(i,j,k), v(i,j,k), w(i,j,k), &
6258 q(i,j,k,sphum)*1000., cond*1000.
6270 subroutine sounding_column( pt, delp, delz, u, v, q, peln, pkz, phis, &
6271 npz, ncnst, sphum, nwat, hydrostatic, moist_phys, zvir, ng, bd, Time )
6274 integer,
intent(IN) :: npz, ncnst, sphum, nwat, ng
6275 real,
intent(IN) :: zvir
6276 logical,
intent(IN) :: hydrostatic, moist_phys
6277 real,
dimension(bd%isd:bd%ied,bd%jsd:bd%jed,npz),
intent(IN) :: pt, delp
6278 real,
dimension(bd%is:, bd%js:, 1:),
intent(IN) :: delz
6279 real,
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz),
intent(IN) :: u
6280 real,
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz),
intent(IN) :: v
6281 real,
dimension(bd%isd:bd%ied, bd%jsd:bd%jed, npz, ncnst),
intent(IN) :: q
6282 real,
dimension(bd%is:bd%ie,npz+1,bd%js:bd%je),
intent(in):: peln
6283 real,
dimension(bd%is:bd%ie,bd%js:bd%je,npz),
intent(in):: pkz
6284 real,
dimension(bd%isd:bd%ied,bd%jsd:bd%jed),
intent(IN) :: phis
6285 type(time_type),
intent(IN) :: Time
6287 real :: Tv, pres, hght(npz), dewpt, rh, mixr, tmp, qs(1), wspd, wdir, rpk, theta, thetav
6288 real :: thetae(bd%is:bd%ie,bd%js:bd%je,npz)
6290 real,
PARAMETER :: rgrav = 1./grav
6291 real,
PARAMETER :: rdg = -rdgas*rgrav
6292 real,
PARAMETER :: sounding_top = 10.e2
6293 real,
PARAMETER :: ms_to_knot = 1.9438445
6294 real,
PARAMETER :: p0 = 1000.e2
6296 integer :: i, j, k, n
6297 integer :: yr_v, mo_v, dy_v, hr_v, mn_v, sec_v
6300 call get_date(time, yr_v, mo_v, dy_v, hr_v, mn_v, sec_v)
6301 call eqv_pot(thetae, pt, delp, delz, peln, pkz, q(bd%isd,bd%jsd,1,sphum), &
6302 bd%is, bd%ie, bd%js, bd%je, ng, npz, hydrostatic, moist_phys)
6309 if (i < bd%is .or. i > bd%ie) cycle
6310 if (j < bd%js .or. j > bd%je) cycle
6318 600
format(a,
'.v', i4, i2.2, i2.2, i2.2,
'.i', i4, i2.2, i2.2, i2.2,
'.', a,
'.dat########################################################')
6319 write(
diag_sonde_units(n),601) trim(
diag_sonde_names(n)), yr_v, mo_v, dy_v, hr_v,
yr_init,
mo_init,
dy_init,
hr_init, &
6321 601
format(3x, a16,
' Valid ', i4, i2.2, i2.2,
'.', i2.2,
'Z Init ', i4, i2.2, i2.2,
'.', i2.2,
'Z \n', a, 2f8.3)
6323 write(
diag_sonde_units(n),*)
'-------------------------------------------------------------------------------' 6324 write(
diag_sonde_units(n),
'(11A7)')
'PRES',
'HGHT',
"TEMP",
"DWPT",
"RELH",
"MIXR",
"DRCT",
"SKNT",
"THTA",
"THTE",
"THTV" 6325 write(
diag_sonde_units(n),
'(11A7)')
'hPa',
'm',
'C',
'C',
'%',
'g/kg',
'deg',
'knot',
'K',
'K',
'K' 6326 write(
diag_sonde_units(n),*)
'-------------------------------------------------------------------------------' 6328 if (hydrostatic)
then 6329 call mpp_error(note,
'Hydrostatic diagnostic sounding not yet supported')
6331 hght(npz) = phis(i,j)*rgrav - 0.5*delz(i,j,npz)
6333 hght(k) = hght(k+1) - 0.5*(delz(i,j,k)+delz(i,j,k+1))
6341 tv = pt(i,j,k)*(1.+zvir*q(i,j,k,sphum))
6343 pres = delp(i,j,k)/delz(i,j,k)*rdg*tv
6347 call qsmith((bd%ie-bd%is+1)*(bd%je-bd%js+1), npz, &
6348 1, 1, pt(i,j,k:k), &
6349 (/pres/), q(i,j,k:k,sphum), qs)
6351 call qsmith(1, 1, 1, pt(i,j,k:k), &
6352 (/pres/), q(i,j,k:k,sphum), qs)
6355 mixr = q(i,j,k,sphum)/(1.-sum(q(i,j,k,1:nwat)))
6356 rh = q(i,j,k,sphum)/qs(1)
6357 tmp = ( log(max(rh,1.e-2))/ 17.27 + ( pt(i,j,k) - 273.14 )/ ( -35.84 + pt(i,j,k)) )
6358 dewpt = 237.3* tmp/ ( 1. - tmp )
6359 wspd = 0.5*sqrt((u(i,j,k)+u(i,j+1,k))*(u(i,j,k)+u(i,j+1,k)) + (v(i,j,k)+v(i+1,j,k))*(v(i,j,k)+v(i+1,j,k)))*ms_to_knot
6360 if (wspd > 0.01)
then 6362 wdir = atan2(u(i,j,k)+u(i,j+1,k),v(i,j,k)+v(i+1,j,k)) *
rad2deg 6366 rpk = exp(-kappa*log(pres/p0))
6367 theta = pt(i,j,k)*rpk
6370 write(
diag_sonde_units(n),
'(F7.1, I7, F7.1, F7.1, I7, F7.2, I7, F7.2, F7.1, F7.1, F7.1)') &
6371 pres*1.e-2, int(hght(k)), pt(i,j,k)-tfreeze, dewpt, int(rh*100.), mixr*1.e3, int(wdir), wspd, theta, thetae(i,j,k), thetav
real, dimension(max_diag_column) diag_debug_lat_in
character(len=256) tlongname
subroutine, public prt_height(qname, is, ie, js, je, ng, km, press, phis, delz, peln, area, lat)
pure real function, public vicpqd(q)
real, parameter missing_value3
for variables where we look for smallest values
real, dimension(:), allocatable diag_debug_lat
subroutine cs_prof(q2, delp, q, km, i1, i2, iv)
The module 'fv_mp_mod' is a single program multiple data (SPMD) parallel decompostion/communication m...
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)
The type 'fv_grid_type' is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
real, dimension(2) rhrange
real, parameter missing_value2
for variables with many missing values
real, dimension(max_diag_column) diag_debug_lon_in
real, dimension(:,:), allocatable, public zs_g
subroutine, public moist_cv(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cvm, t1)
The subroutine 'moist_cv' computes the FV3-consistent moist heat capacity under constant volume...
real, public sphum_ll_fix
real, dimension(2) vrange
subroutine getcape(nk, p, t, dz, q, the, cape, cin, source_in)
The subroutine 'getcape' calculates the Convective Available Potential Energy (CAPE) from a Sounding...
The module 'multi_gases' peforms multi constitutents computations.
subroutine, public a2b_ord2(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
subroutine, public eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, hydrostatic, moist)
The subroutine 'eqv_pot' calculates the equivalent potential temperature.
subroutine, public get_eta_level(npz, p_s, pf, ph, ak, bk, pscale)
The subroutine 'get_eta_level' returns the interface and layer-mean pressures for reference...
real, dimension(2) skrange
logical, dimension(:,:), allocatable do_sonde_diag_column
real, dimension(2) trange
real, parameter missing_value
subroutine get_height_given_pressure(is, ie, js, je, km, wz, kd, id, log_p, peln, a2)
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function 'g_sum' is the fast version of 'globalsum'.
integer, parameter max_step
subroutine wind_max(isc, iec, jsc, jec, isd, ied, jsd, jed, us, vs, ws_max, domain)
integer, dimension(:), allocatable diag_debug_i
character(16), dimension(max_diag_column) diag_sonde_names
real, dimension(2) wrange
subroutine get_pressure_given_height(is, ie, js, je, ng, km, wz, kd, height, ts, peln, a2, fac)
integer, parameter, public r_grid
character(len=256) tunits
subroutine, public gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
logical, dimension(:,:), allocatable do_debug_diag_column
The module 'fv_sg' performs FV sub-grid mixing.
subroutine, public bunkers_vector(is, ie, js, je, ng, km, zvir, sphum, uc, vc, ua, va, delz, q, hydrostatic, pt, peln, phis, grav)
subroutine, public fv_diag_init_gn(Atm)
pure real function, public virqd(q)
real, dimension(max_diag_column) diag_sonde_lon_in
subroutine, public prt_mxm(qname, q, is, ie, js, je, n_g, km, fac, area, domain)
subroutine, public max_vorticity_hy1(is, ie, js, je, km, vort, maxvorthy1)
real, dimension(:), allocatable diag_sonde_lon
pure real function, public virq(q)
subroutine range_check_2d(qname, q, is, ie, js, je, n_g, pos, q_low, q_hi, bad_range, Time)
subroutine interpolate_z(is, ie, js, je, km, zl, hght, a3, a2)
The module 'a2b_edge' performs FV-consistent interpolation of pressure to corners.
subroutine, public z_sum(is, ie, js, je, km, n_g, delp, q, sum2)
subroutine, public fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
subroutine cs_interpolator(is, ie, js, je, km, qin, zout, wz, qout, qmin)
subroutine cs3_interpolator(is, ie, js, je, km, qin, kd, pout, wz, pe, id, qout, iv)
subroutine, public ppme(p, qe, delp, im, km)
subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
The module 'fv_mapz' contains the vertical mapping routines .
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
logical, public prt_minmax
The module 'fv_arrays' contains the 'fv_atmos_type' and associated datatypes.
integer, dimension(:), allocatable diag_sonde_i
subroutine, public get_height_field(is, ie, js, je, ng, km, hydrostatic, delz, wz, pt, q, peln, zvir)
The module 'fv_eta' contains routine to set up the reference (Eulerian) pressure coordinate.
subroutine debug_column(pt, delp, delz, u, v, w, q, npz, ncnst, sphum, nwat, hydrostatic, bd, Time)
subroutine, public fv_diag(Atm, zvir, Time, print_freq)
real, dimension(:), allocatable diag_sonde_lat
subroutine updraft_helicity(is, ie, js, je, ng, km, zvir, sphum, uh, w, vort, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
subroutine pv_entropy(is, ie, js, je, ng, km, vort, f_d, pt, pkz, delp, grav)
The subroutine 'pv_entropy' computes potential vorticity.
subroutine, public moist_cp(is, ie, isd, ied, jsd, jed, km, j, k, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, q, qd, cpm, t1)
The subroutine 'moist_cp' computes the FV3-consistent moist heat capacity under constant pressure...
type(time_type), public fv_time
logical module_is_initialized
real function getqvi(p, t)
real function getqvs(p, t)
integer, dimension(:), allocatable diag_sonde_units
integer, dimension(nplev) levs
subroutine, public dbzcalc(q, pt, delp, peln, delz, dbz, maxdbz, allmax, bd, npz, ncnst, hydrostatic, zvir, in0r, in0s, in0g, iliqskin)
The subroutine 'dbzcalc' computes equivalent reflectivity factor (in dBZ) at each model grid point...
@ The module 'fv_diagnostics' contains routines to compute diagnosic fields.
The module 'fv_grid_utils' contains routines for setting up and computing grid-related quantities...
character(len=3), public gn
subroutine, public prt_gb_nh_sh(qname, is, ie, js, je, a2, area, lat)
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
real, dimension(2) vsrange
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 rh_calc(pfull, t, qv, rh, do_cmip)
subroutine sounding_column(pt, delp, delz, u, v, q, peln, pkz, phis, npz, ncnst, sphum, nwat, hydrostatic, moist_phys, zvir, ng, bd, Time)
character(16), dimension(max_diag_column) diag_debug_names
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)
subroutine, public a2b_ord4(qin, qout, gridstruct, npx, npy, is, ie, js, je, ng, replace)
subroutine range_check_3d(qname, q, is, ie, js, je, n_g, km, pos, q_low, q_hi, bad_range, Time)
real, dimension(max_diag_column) diag_sonde_lat_in
real function p_sum(is, ie, js, je, km, n_g, delp, area, domain)
integer, parameter max_diag_column
type(fv_diag_type), pointer idiag
subroutine, public interpolate_vertical(is, ie, js, je, km, plev, peln, a3, a2)
real(kind=4), public e_flux
integer, dimension(:), allocatable diag_debug_units
real, dimension(:), allocatable diag_debug_lon
subroutine, public prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain)
subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, w, delz, pt, delp, q, hs, area, domain, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, nwat, ua, va, moist_phys, te)
The subroutine 'nh_total_energy computes vertically-integrated total energy per column.
subroutine, public max_vv(is, ie, js, je, npz, ng, up2, dn2, pe, w)
subroutine helicity_relative(is, ie, js, je, ng, km, zvir, sphum, srh, ua, va, delz, q, hydrostatic, pt, peln, phis, grav, z_bot, z_top)
real, dimension(2) slprange
integer, dimension(:), allocatable diag_debug_j
pure real function, public vicvqd(q)
integer, dimension(:), allocatable diag_sonde_j
subroutine, public get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)