100 use constants_mod,
only: cnst_radius=>
radius, pi=>pi_8, omega, grav, kappa, rdgas, cp_air, rvgas
103 is,js,ie,je, isd,jsd,ied,jed, &
104 domain_decomp, fill_corners, xdir, ydir, &
105 mp_stop, mp_reduce_sum, mp_reduce_max, mp_gather, mp_bcst
115 use mpp_mod,
only: mpp_error, fatal, mpp_root_pe, mpp_broadcast, mpp_sum
116 use mpp_domains_mod,
only: mpp_update_domains, domain2d
117 use mpp_parameter_mod,
only: agrid_param=>agrid,cgrid_ne_param=>cgrid_ne, &
122 use mpp_mod,
only: mpp_pe, mpp_chksum, stdout
125 use tracer_manager_mod,
only: get_tracer_index
126 use field_manager_mod,
only: model_atmos
171 real(kind=R_GRID),
parameter ::
radius = cnst_radius
172 real(kind=R_GRID),
parameter ::
one = 1.d0
196 real,
allocatable,
dimension(:) ::
pz0,
zz0 211 real ,
allocatable ::
phi0(:,:,:)
212 real ,
allocatable ::
ua0(:,:,:)
213 real ,
allocatable ::
va0(:,:,:)
229 public :: output, output_ncdf
244 subroutine init_winds(UBar, u,v,ua,va,uc,vc, defOnGrid, npx, npy, ng, ndims, nregions, nested, gridstruct, domain, tile)
247 real ,
intent(INOUT) :: UBar
248 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
249 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
250 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed )
251 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1)
252 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed )
253 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed )
254 integer,
intent(IN) :: defOnGrid
255 integer,
intent(IN) :: npx, npy
256 integer,
intent(IN) :: ng
257 integer,
intent(IN) :: ndims
258 integer,
intent(IN) :: nregions
259 logical,
intent(IN) :: nested
260 type(fv_grid_type),
intent(IN),
target :: gridstruct
261 type(domain2d),
intent(INOUT) :: domain
262 integer,
intent(IN) :: tile
264 real(kind=R_GRID) :: p1(2), p2(2), p3(2), p4(2), pt(2)
265 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3)
271 real :: psi_b(isd:ied+1,jsd:jed+1), psi(isd:ied,jsd:jed), psi1, psi2
272 integer :: is2, ie2, js2, je2
274 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
275 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
276 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
277 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
278 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
280 logical,
pointer :: cubed_sphere, latlon
282 logical,
pointer :: have_south_pole, have_north_pole
284 integer,
pointer :: ntiles_g
285 real,
pointer :: acapN, acapS, globalarea
287 grid => gridstruct%grid_64
288 agrid=> gridstruct%agrid_64
290 area => gridstruct%area
291 rarea => gridstruct%rarea
296 ee1 => gridstruct%ee1
297 ee2 => gridstruct%ee2
300 en1 => gridstruct%en1
301 en2 => gridstruct%en2
305 dxa => gridstruct%dxa
306 dya => gridstruct%dya
307 rdxa => gridstruct%rdxa
308 rdya => gridstruct%rdya
309 dxc => gridstruct%dxc
310 dyc => gridstruct%dyc
312 cubed_sphere => gridstruct%cubed_sphere
313 latlon => gridstruct%latlon
315 have_south_pole => gridstruct%have_south_pole
316 have_north_pole => gridstruct%have_north_pole
318 ntiles_g => gridstruct%ntiles_g
319 acapn => gridstruct%acapN
320 acaps => gridstruct%acapS
321 globalarea => gridstruct%globalarea
339 200
format(i4.4,
'x',i4.4,
'x',i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
345 psi(i,j) = (-1.0 * ubar *
radius *( sin(agrid(i,j,2)) *cos(
alpha) - &
346 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
349 call mpp_update_domains( psi, domain )
352 psi_b(i,j) = (-1.0 * ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
353 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
357 if ( (cubed_sphere) .and. (defongrid==0) )
then 361 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
362 if (dist==0) vc(i,j) = 0.
368 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
369 if (dist==0) uc(i,j) = 0.
372 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
373 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
377 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
378 if (dist==0) v(i,j) = 0.
384 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
385 if (dist==0) u(i,j) = 0.
391 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
392 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
394 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
395 if (dist==0) ua(i,j) = 0.
396 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
397 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
399 va(i,j) = (psi2 - psi1) / (dist)
400 if (dist==0) va(i,j) = 0.
404 elseif ( (cubed_sphere) .and. (defongrid==1) )
then 408 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
409 if (dist==0) vc(i,j) = 0.
415 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
416 if (dist==0) uc(i,j) = 0.
419 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
420 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
421 call ctoa(uc,vc,ua,va,dx, dy, dxc,dyc,dxa,dya,npx,npy,ng)
422 call atod(ua,va,u ,v ,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
425 elseif ( (cubed_sphere) .and. (defongrid==2) )
then 429 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
430 if (dist==0) v(i,j) = 0.
436 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
437 if (dist==0) u(i,j) = 0.
441 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng)
442 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
443 elseif ( (cubed_sphere) .and. (defongrid==3) )
then 446 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
447 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
449 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
450 if (dist==0) ua(i,j) = 0.
451 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
452 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
454 va(i,j) = (psi2 - psi1) / (dist)
455 if (dist==0) va(i,j) = 0.
458 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
459 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
460 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested,domain)
461 elseif ( (latlon) .or. (defongrid==4) )
then 465 ua(i,j) = ubar * ( cos(agrid(i,j,2))*cos(
alpha) + &
466 sin(agrid(i,j,2))*cos(agrid(i,j,1))*sin(
alpha) )
467 va(i,j) = -ubar * sin(agrid(i,j,1))*sin(
alpha)
472 if (cubed_sphere)
call rotate_winds(ua(i,j), va(i,j), p1,p2,p3,p4, agrid(i,j,1:2), 2, 1)
474 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
475 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
477 if ( (tile==1) .and.(i==1) ) print*, ua(i,j), -1.0 * (psi2 - psi1) / (dist)
481 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
482 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
483 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
484 elseif ( (latlon) .or. (defongrid==5) )
then 489 p1(:) = grid(i ,j ,1:2)
490 p2(:) = grid(i,j+1 ,1:2)
494 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
495 sin(pt(2))*cos(pt(1))*sin(
alpha) )
496 vtmp = -ubar * sin(pt(1))*sin(
alpha)
503 p1(:) = grid(i ,j ,1:2)
504 p2(:) = grid(i+1,j ,1:2)
508 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
509 sin(pt(2))*cos(pt(1))*sin(
alpha) )
510 vtmp = -ubar * sin(pt(1))*sin(
alpha)
516 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng)
517 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
537 subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
538 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, &
539 dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, adiabatic, &
540 ks, npx_global, ptop, domain_in, tile_in, bd)
543 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
544 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
545 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
546 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
547 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
548 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
550 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
552 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
553 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
554 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
555 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
556 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
558 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
559 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
560 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
561 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
562 real ,
intent(inout) :: delz(bd%isd:,bd%jsd:,1:)
563 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
565 real ,
intent(inout) :: ak(npz+1)
566 real ,
intent(inout) :: bk(npz+1)
568 integer,
intent(IN) :: npx, npy, npz
569 integer,
intent(IN) ::
ng, ncnst, nwat
570 integer,
intent(IN) :: ndims
571 integer,
intent(IN) :: nregions
573 real,
intent(IN) :: dry_mass
574 logical,
intent(IN) :: mountain
575 logical,
intent(IN) :: moist_phys
576 logical,
intent(IN) :: hydrostatic
577 logical,
intent(IN) :: hybrid_z
578 logical,
intent(IN) :: adiabatic
579 integer,
intent(IN) :: ks
584 integer,
intent(IN) :: npx_global
585 integer,
intent(IN),
target :: tile_in
586 real,
intent(INOUT) :: ptop
588 type(domain2d),
intent(IN),
target :: domain_in
590 real :: tmp(1-
ng:npx +
ng,1-
ng:npy +
ng,1:nregions)
591 real :: tmp1(1 :npx ,1 :npy ,1:nregions)
593 real(kind=R_GRID) :: p0(2)
594 real(kind=R_GRID) :: p1(2)
595 real(kind=R_GRID) :: p2(2)
596 real(kind=R_GRID) :: p3(2)
597 real(kind=R_GRID) :: p4(2)
598 real(kind=R_GRID) :: pa(2)
599 real(kind=R_GRID) :: pb(2)
600 real(kind=R_GRID) :: pcen(2)
601 real(kind=R_GRID) :: e1(3), e2(3), e3(3), ex(3), ey(3)
602 real :: dist, r, r1, r2, r0, omg, a, b, c
603 integer :: i,j,k,nreg,z,zz
604 integer :: i0,j0,n0, nt
605 real :: utmp,vtmp,ftmp
608 integer,
parameter :: jm = 5761
615 real :: ddeg, deg, ddp, dp, ph5
620 real :: x1,y1,z1,x2,y2,z2,ang
622 integer :: initwindscase
634 real :: grad(bd%isd:bd%ied ,bd%jsd:bd%jed,2)
635 real :: div0(bd%isd:bd%ied ,bd%jsd:bd%jed )
636 real :: vor0(bd%isd:bd%ied ,bd%jsd:bd%jed )
637 real :: divg(bd%isd:bd%ied ,bd%jsd:bd%jed )
638 real :: vort(bd%isd:bd%ied ,bd%jsd:bd%jed )
639 real :: ztop, rgrav, p00, pturb, zmid, pk0, t00
640 real :: dz1(npz), ppt(npz)
641 real :: ze1(npz+1), pe1(npz+1)
644 character(len=80) :: oflnm, hgtflnm
645 integer :: is2, ie2, js2, je2
647 real :: psi(bd%isd:bd%ied,bd%jsd:bd%jed)
648 real :: psi_b(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
652 real :: eta(npz), eta_0, eta_s, eta_t
653 real :: eta_v(npz), press, anti_rot
654 real :: t_0, t_mean, delta_t, lapse_rate, n2, zeta, s0
655 real :: pt1,pt2,pt3,pt4,pt5,pt6, pt7, pt8, pt9, u1, pt0
656 real :: uu1, uu2, uu3, vv1, vv2, vv3
659 real wbuffer(npy+2,npz)
660 real sbuffer(npx+2,npz)
662 real :: gz(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1), zt, zdist
669 real,
dimension(npz):: pk1, ts1, qs1, uz1, zs1, dudz
671 real(kind=R_GRID):: pp0(2)
676 real :: omega0, k_cell, z0, h, px
677 real :: d1, d2, p1p(2), rt, s
678 real :: wind_alpha, period, h0, rm, zp3(3), dz3(3), k0, lp
682 real,
dimension(npz+1) :: pe0, gz0, ue, ve, we, pte, qe
683 real :: d, cor, exppr, exppz, gamma, ts0, q00, exponent, ztrop, height, zp, rp
684 real :: qtrop, ttrop, zq1, zq2
685 real :: dum, dum1, dum2, dum3, dum4, dum5, dum6, ptmp, uetmp, vetmp
686 real :: pe_u(bd%is:bd%ie,npz+1,bd%js:bd%je+1)
687 real :: pe_v(bd%is:bd%ie+1,npz+1,bd%js:bd%je)
688 real :: ps_u(bd%is:bd%ie,bd%js:bd%je+1)
689 real :: ps_v(bd%is:bd%ie+1,bd%js:bd%je)
694 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
695 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
696 real,
pointer,
dimension(:,:) :: rarea, fc, f0
697 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
698 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
699 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
701 logical,
pointer :: cubed_sphere, latlon
703 type(domain2d),
pointer :: domain
704 integer,
pointer :: tile
706 logical,
pointer :: have_south_pole, have_north_pole
708 integer,
pointer :: ntiles_g
709 real,
pointer :: acapn, acaps, globalarea
720 grid => gridstruct%grid_64
721 agrid=> gridstruct%agrid_64
723 area => gridstruct%area_64
724 rarea => gridstruct%rarea
729 ee1 => gridstruct%ee1
730 ee2 => gridstruct%ee2
733 en1 => gridstruct%en1
734 en2 => gridstruct%en2
738 dxa => gridstruct%dxa
739 dya => gridstruct%dya
740 rdxa => gridstruct%rdxa
741 rdya => gridstruct%rdya
742 dxc => gridstruct%dxc
743 dyc => gridstruct%dyc
745 cubed_sphere => gridstruct%cubed_sphere
746 latlon => gridstruct%latlon
751 have_south_pole => gridstruct%have_south_pole
752 have_north_pole => gridstruct%have_north_pole
754 ntiles_g => gridstruct%ntiles_g
755 acapn => gridstruct%acapN
756 acaps => gridstruct%acapS
757 globalarea => gridstruct%globalarea
759 if (gridstruct%nested)
then 773 f0(:,:) = huge(dummy)
774 fc(:,:) = huge(dummy)
777 fc(i,j) = 2.*omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
778 sin(grid(i,j,2))*cos(
alpha) )
783 f0(i,j) = 2.*omega*( -1.*cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) + &
784 sin(agrid(i,j,2))*cos(
alpha) )
787 call mpp_update_domains( f0, domain )
788 if (cubed_sphere)
call fill_corners(f0, npx, npy, ydir)
790 delp(isd:is-1,jsd:js-1,1:npz)=0.
791 delp(isd:is-1,je+1:jed,1:npz)=0.
792 delp(ie+1:ied,jsd:js-1,1:npz)=0.
793 delp(ie+1:ied,je+1:jed,1:npz)=0.
795 #if defined(SW_DYNAMICS) 805 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
806 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
809 call init_winds(
ubar, u,v,ua,va,uc,vc, 1, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
814 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
815 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
816 if ( (tile==1) .and. (i==1) )
write(*,200) i,j,tile, divg(i,j), uc(i,j,1), uc(i+1,j,1), vc(i,j,1), vc(i,j+1,1)
822 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
823 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
832 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
833 200
format(i4.4,
'x',i4.4,
'x',i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
834 201
format(
' ',a,e21.14,
' ',e21.14)
835 202
format(
' ',a,i4.4,
'x',i4.4,
'x',i4.4)
836 if ( is_master() )
then 837 write(*,*)
' Error Norms of Analytical Divergence field C-Winds initialized' 838 write(*,201)
'Divergence MAX error : ', pmax
839 write(*,201)
'Divergence MIN error : ', pmin
840 write(*,201)
'Divergence L1_norm : ', l1_norm
841 write(*,201)
'Divergence L2_norm : ', l2_norm
842 write(*,201)
'Divergence Linf_norm : ', linf_norm
845 call init_winds(
ubar, u,v,ua,va,uc,vc, 3, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
849 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
850 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
851 if ( (tile==1) .and. (i==1) )
write(*,200) i,j,tile, divg(i,j), uc(i,j,1), uc(i+1,j,1), vc(i,j,1), vc(i,j+1,1)
857 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
858 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
865 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
866 if ( is_master() )
then 867 write(*,*)
' Error Norms of Analytical Divergence field A-Winds initialized' 868 write(*,201)
'Divergence MAX error : ', pmax
869 write(*,201)
'Divergence MIN error : ', pmin
870 write(*,201)
'Divergence L1_norm : ', l1_norm
871 write(*,201)
'Divergence L2_norm : ', l2_norm
872 write(*,201)
'Divergence Linf_norm : ', linf_norm
875 call init_winds(
ubar, u,v,ua,va,uc,vc, 2, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
881 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
882 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
883 if ( (tile==1) .and. ((i==1) .or.(i==npx-1)) )
write(*,200) i,j,tile, divg(i,j), uc(i,j,1), uc(i+1,j,1), vc(i,j,1), vc(i,j+1,1)
889 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
890 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
895 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
896 if ( is_master() )
then 897 write(*,*)
' Error Norms of Analytical Divergence field D-Winds initialized' 898 write(*,201)
'Divergence MAX error : ', pmax
899 write(*,201)
'Divergence MIN error : ', pmin
900 write(*,201)
'Divergence L1_norm : ', l1_norm
901 write(*,201)
'Divergence L2_norm : ', l2_norm
902 write(*,201)
'Divergence Linf_norm : ', linf_norm
916 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
918 if (p /= 0.0) w_p = vtx/p
919 delp(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*0.0) )
920 ua(i,j,1) = w_p*(sin(
lat0)*cos(agrid(i,j,2)) + cos(
lat0)*cos(agrid(i,j,1) -
lon0)*sin(agrid(i,j,2)))
921 va(i,j,1) = w_p*cos(
lat0)*sin(agrid(i,j,1) -
lon0)
922 ua(i,j,1) = ua(i,j,1)*
radius/86400.0
923 va(i,j,1) = va(i,j,1)*
radius/86400.0
929 if (cubed_sphere)
call rotate_winds(ua(i,j,1),va(i,j,1), p1,p2,p3,p4, agrid(i,j,1:2), 2, 1)
933 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
934 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,
ng, gridstruct%nested, domain)
936 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
937 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
938 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
953 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
955 delp(i,j,1) = phis(i,j)
976 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.))
then 994 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
995 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
997 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
998 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0
1000 ( -1.*cos(grid(i+1,j ,1))*cos(grid(i+1,j ,2))*sin(
alpha) + &
1001 sin(grid(i+1,j ,2))*cos(
alpha) ) ** 2.0
1003 ( -1.*cos(grid(i+1,j+1,1))*cos(grid(i+1,j+1,2))*sin(
alpha) + &
1004 sin(grid(i+1,j+1,2))*cos(
alpha) ) ** 2.0
1006 ( -1.*cos(grid(i,j+1,1))*cos(grid(i,j+1,2))*sin(
alpha) + &
1007 sin(grid(i,j+1,2))*cos(
alpha) ) ** 2.0
1008 delp(i,j,1) = (0.25*(pt1+pt2+pt3+pt4) + 3.*pt5) / 4.
1011 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
1012 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
1033 p2(1) = agrid(i,j,1)
1034 p2(2) = agrid(i,j,2)
1037 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
1039 delp(i,j,1) = phis(i,j)
1042 delp(i,j,1) = delp(i,j,1) + grav*2.e3
1053 p1(:) = grid(i ,j ,1:2)
1054 p2(:) = grid(i,j+1 ,1:2)
1058 utmp =
ubar * cos(p3(2))
1065 p1(:) = grid(i, j,1:2)
1066 p2(:) = grid(i+1,j,1:2)
1070 utmp =
ubar * cos(p3(2))
1079 fc(i,j) = 2.*anti_rot*sin(grid(i,j,2))
1084 f0(i,j) = 2.*anti_rot*sin(agrid(i,j,2))
1113 p1(1) = pi*1.5 - ddeg
1117 p2(1) = pi*1.5 + ddeg
1121 #ifndef SINGULAR_VORTEX 1153 p2(1) = agrid(i,j,1)
1154 p2(2) = agrid(i,j,2)
1155 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1157 phis(i,j) = 2000.0*grav*(1.0-(r/r0))
1163 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
1164 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2 - phis(i,j)
1176 a = 0.5*omg*(2.*omega+omg)*(cos(agrid(i,j,2))**2) + &
1177 0.25*rk*rk*(cos(agrid(i,j,2))**(r+r)) * &
1178 ( (r+1)*(cos(agrid(i,j,2))**2) + (2.*r*r-r-2.) - &
1179 2.*(r*r)*cos(agrid(i,j,2))**(-2.) )
1180 b = (2.*(omega+omg)*rk / ((r+1)*(r+2))) * (cos(agrid(i,j,2))**r) * &
1181 ( (r*r+2.*r+2.) - ((r+1.)*cos(agrid(i,j,2)))**2 )
1182 c = 0.25*rk*rk*(cos(agrid(i,j,2))**(2.*r)) * ( &
1183 (r+1) * (cos(agrid(i,j,2))**2.) - (r+2.) )
1184 delp(i,j,1) =
gh0 +
radius*
radius*(a+b*cos(r*agrid(i,j,1))+c*cos(2.*r*agrid(i,j,1)))
1185 delp(i,j,1) = delp(i,j,1) - phis(i,j)
1190 p1(:) = grid(i ,j ,1:2)
1191 p2(:) = grid(i,j+1 ,1:2)
1195 utmp =
radius*omg*cos(p3(2)) + &
1196 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1197 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1203 p1(:) = grid(i, j,1:2)
1204 p2(:) = grid(i+1,j,1:2)
1208 utmp =
radius*omg*cos(p3(2)) + &
1209 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1210 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1215 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,
ng)
1217 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
1235 pt1 =
gh_jet(npy, agrid(i,j,2))
1244 pt6 =
gh_jet(npy, grid(i, j, 2))
1245 pt7 =
gh_jet(npy, grid(i+1,j, 2))
1246 pt8 =
gh_jet(npy, grid(i+1,j+1,2))
1247 pt9 =
gh_jet(npy, grid(i ,j+1,2))
1248 ftmp = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1250 delp(i,j,1) = ftmp + 120.*grav*cos(agrid(i,j,2)) * &
1251 exp( -(3.*(agrid(i,j,1)-pi))**2 ) * exp( -(15.*(agrid(i,j,2)-pi/4.))**2 )
1257 p1(:) = agrid(i,j,1:2)
1260 if ( r < 3.*r0 )
then 1261 delp(i,j,1) = delp(i,j,1) + 1000.*grav*exp(-(r/r0)**2)
1270 p2(:) = grid(i,j+1,1:2)
1271 vv1 =
u_jet(p2(2))*(ee2(2,i,j+1)*cos(p2(1)) - ee2(1,i,j+1)*sin(p2(1)))
1272 p1(:) = grid(i,j,1:2)
1273 vv3 =
u_jet(p1(2))*(ee2(2,i,j)*cos(p1(1)) - ee2(1,i,j)*sin(p1(1)))
1276 vv2 =
u_jet(pa(2))*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1278 v(i,j,1) = 0.25*(vv1 + 2.*vv2 + vv3)
1285 p1(:) = grid(i,j,1:2)
1286 uu1 =
u_jet(p1(2))*(ee1(2,i,j)*cos(p1(1)) - ee1(1,i,j)*sin(p1(1)))
1287 p2(:) = grid(i+1,j,1:2)
1288 uu3 =
u_jet(p2(2))*(ee1(2,i+1,j)*cos(p2(1)) - ee1(1,i+1,j)*sin(p2(1)))
1291 uu2 =
u_jet(pa(2))*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1293 u(i,j,1) = 0.25*(uu1 + 2.*uu2 + uu3)
1300 call get_vorticity(is, ie, js, je, isd, ied, jsd, jed, npz, u, v, q(is:ie,js:je,:,1), dx, dy, rarea)
1303 fc(i,j) = 2.*omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
1304 sin(grid(i,j,2))*cos(
alpha) )
1309 f0(i,j) = 2.*omega*( -1.*cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) + &
1310 sin(agrid(i,j,2))*cos(
alpha) )
1313 call mpp_update_domains( f0, domain )
1314 if (cubed_sphere)
call fill_corners(f0, npx, npy, ydir)
1317 q(i,j,npz,1) = ( q(i,j,npz,1) + f0(i,j) ) / delp(i,j,npz) * 1.e6
1335 p2(1) = agrid(i,j,1)
1336 p2(2) = agrid(i,j,2)
1337 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1339 phis(i,j) = 2000.0*grav*(1.0-(r/r0))
1354 if ( is_master() )
write(*,*)
'Initialzing case-8: soliton twin cycolne...' 1377 p1(:) = grid(i ,j ,1:2)
1378 p2(:) = grid(i,j+1 ,1:2)
1381 utmp =
ubar*exp(-(r/r0)**2)
1389 p1(:) = grid(i, j,1:2)
1390 p2(:) = grid(i+1,j,1:2)
1393 utmp =
ubar*exp(-(r/r0)**2)
1406 p1(:) = grid(i ,j ,1:2)
1407 p2(:) = grid(i,j+1 ,1:2)
1410 utmp =
ubar*exp(-(r/r0)**2)
1418 p1(:) = grid(i, j,1:2)
1419 p2(:) = grid(i+1,j,1:2)
1422 utmp =
ubar*exp(-(r/r0)**2)
1437 ph5 = -0.5*pi + (dble(j-1)-0.5)*ddp
1438 ll_j(j) = -0.5*pi + (dble(j-1)*ddp)
1444 cosp(j) = (sine(j+1)-sine(j)) / dp
1447 cose(j) = 0.5 * (cosp(j-1) + cosp(j))
1450 ddeg = 180./float(jm-1)
1452 deg = -90. + (float(j-1)-0.5)*ddeg
1454 ll_u(j) = -10.*(deg+90.)/90.
1455 elseif (deg <= 60.)
then 1456 ll_u(j) = -10. + deg
1458 ll_u(j) = 50. - (50./30.)* (deg - 60.)
1461 ll_phi(1) = 6000. * grav
1463 ll_phi(j)=ll_phi(j-1) - dp*sine(j) * &
1464 (
radius*2.*omega + ll_u(j)/cose(j))*ll_u(j)
1470 if ( (ll_j(jj) <= agrid(i,j,2)) .and. (agrid(i,j,2) <= ll_j(jj+1)) )
then 1471 delp(i,j,1)=0.5*(ll_phi(jj)+ll_phi(jj+1))
1479 if (agrid(i,j,2)*
todeg <= 0.0)
then 1480 ua(i,j,1) = -10.*(agrid(i,j,2)*
todeg + 90.)/90.
1481 elseif (agrid(i,j,2)*
todeg <= 60.0)
then 1482 ua(i,j,1) = -10. + agrid(i,j,2)*
todeg 1484 ua(i,j,1) = 50. - (50./30.)* (agrid(i,j,2)*
todeg - 60.)
1491 if (cubed_sphere)
call rotate_winds(ua(i,j,1), va(i,j,1), p1,p2,p3,p4, agrid(i,j,1:2), 2, 1)
1495 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
1496 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
1497 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
1498 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
1499 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,
ng, gridstruct%nested, domain)
1510 if ( is_master() )
write(*,*)
'Initialzing case-9: soliton cyclones...' 1532 p1(:) = grid(i ,j ,1:2)
1533 p2(:) = grid(i,j+1 ,1:2)
1536 utmp =
ubar*exp(-(r/r0)**2)
1544 p1(:) = grid(i, j,1:2)
1545 p2(:) = grid(i+1,j,1:2)
1548 utmp =
ubar*exp(-(r/r0)**2)
1561 delp(:,:,z) = delp(:,:,1)
1564 call mpp_update_domains( delp, domain )
1565 call mpp_update_domains( phis, domain )
1568 call init_winds(
ubar, u,v,ua,va,uc,vc, initwindscase, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
1577 ps(i,j) = delp(i,j,1)
1591 if (.not.hydrostatic) w(:,:,:)= 0.0
1596 call hydro_eq(npz, is, ie, js, je, ps, phis, 1.e5, &
1597 delp, ak, bk, pt, delz, area,
ng, .false., hydrostatic, hybrid_z, domain)
1604 p1(2) = pi/6. + (7.5/180.0)*pi
1607 p2(1) = agrid(i,j,1)
1608 p2(2) = agrid(i,j,2)
1609 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1611 phis(i,j) =
gh0*(1.0-(r/r0))
1614 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1615 delp, ak, bk, pt, delz, area,
ng, mountain, hydrostatic, hybrid_z, domain)
1619 call surfdrv(npx, npy, gridstruct%grid_64, gridstruct%agrid_64, &
1620 gridstruct%area_64, dx, dy, dxa, dya, dxc, dyc, &
1621 gridstruct%sin_sg, phis, &
1622 flagstruct%stretch_fac, gridstruct%nested, &
1623 npx_global, domain, flagstruct%grid_number, bd, flagstruct%regional)
1624 call mpp_update_domains( phis, domain )
1626 if ( hybrid_z )
then 1632 if ( is_master() )
write(*,*)
'Using const DZ' 1634 dz1(1) = ztop /
real(npz) 1635 dz1(npz) = 0.5*dz1(1)
1654 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1655 delp, ak, bk, pt, delz, area,
ng, mountain, hydrostatic, hybrid_z, domain)
1660 if (is_master()) print*,
'TEST TRACER enabled for this test case' 1663 ncnst, npz, q, agrid(is:ie,js:je,1), agrid(is:ie,js:je,2), 9., 9.)
1670 p1(1) = 205.*pi/180.
1674 p2(1) = agrid(i,j,1)
1675 p2(2) = agrid(i,j,2)
1677 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.) .and. k > 16)
then 1693 cl = get_tracer_index(model_atmos,
'cl')
1694 cl2 = get_tracer_index(model_atmos,
'cl2')
1695 if (cl > 0 .and. cl2 > 0)
then 1697 q, delp,ncnst,agrid(isd:ied,jsd:jed,1),agrid(isd:ied,jsd:jed,2))
1698 call mpp_update_domains(q,domain)
1710 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
1719 peln(i,1,j) = log(ptop)
1720 pk(i,j,1) = ptop**kappa
1725 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
1726 pk(i,j,k) = exp( kappa*log(pe(i,k,j)) )
1727 peln(i,k,j) = log(pe(i,k,j))
1736 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
1744 eta(k) = 0.5*( (ak(k)+ak(k+1))/1.e5 + bk(k)+bk(k+1) )
1745 eta_v(k) = (eta(k) - eta_0)*pi*0.5
1748 if ( .not. adiabatic )
then 1750 sphum = get_tracer_index(model_atmos,
'sphum')
1761 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j)) - 100000.
1762 q(i,j,k,
sphum) = 0.021*exp(-(agrid(i,j,2)/pcen(2))**4.)*exp(-(ptmp/34000.)**2.)
1791 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j+1,2))**2.0
1794 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1795 vv1 = utmp*(ee2(2,i,j+1)*cos(grid(i,j+1,1)) - ee2(1,i,j+1)*sin(grid(i,j+1,1)))
1797 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1800 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1801 vv3 = utmp*(ee2(2,i,j)*cos(grid(i,j,1)) - ee2(1,i,j)*sin(grid(i,j,1)))
1803 p1(:) = grid(i ,j ,1:2)
1804 p2(:) = grid(i,j+1 ,1:2)
1806 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1809 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1810 vv2 = utmp*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1812 v(i,j,z) = 0.25*(vv1 + 2.*vv2 + vv3)
1817 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1820 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1821 uu1 = utmp*(ee1(2,i,j)*cos(grid(i,j,1)) - ee1(1,i,j)*sin(grid(i,j,1)))
1823 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i+1,j,2))**2.0
1826 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1827 uu3 = utmp*(ee1(2,i+1,j)*cos(grid(i+1,j,1)) - ee1(1,i+1,j)*sin(grid(i+1,j,1)))
1829 p1(:) = grid(i ,j ,1:2)
1830 p2(:) = grid(i+1,j ,1:2)
1832 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1835 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1836 uu2 = utmp*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1838 u(i,j,z) = 0.25*(uu1 + 2.*uu2 + uu3)
1853 eta(z) = 0.5*( (ak(z)+ak(z+1))/1.e5 + bk(z)+bk(z+1) )
1855 t_mean = t_0 * eta(z)**(rdgas*lapse_rate/grav)
1856 if (eta_t > eta(z)) t_mean = t_mean + delta_t*(eta_t - eta(z))**5.0
1858 230
format(i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
1861 press = press + delp(is,js,zz)
1863 if (is_master())
write(*,230) z, eta(z), press/100., t_mean
1867 pt1 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1868 ( -2.0*(sin(agrid(i,j,2))**6.0) *(cos(agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1869 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1870 ( (8.0/5.0)*(cos(agrid(i,j,2))**3.0)*(sin(agrid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1881 pt2 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1882 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1883 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1884 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1886 pt3 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1887 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1888 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1889 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1891 pt4 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1892 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1893 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1894 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1896 pt5 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1897 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1898 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1899 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1901 pt6 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1902 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1903 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1904 ( (8.0/5.0)*(cos(grid(i,j,2))**3.0)*(sin(grid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1905 pt7 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1906 ( -2.0*(sin(grid(i+1,j,2))**6.0) *(cos(grid(i+1,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1907 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1908 ( (8.0/5.0)*(cos(grid(i+1,j,2))**3.0)*(sin(grid(i+1,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1909 pt8 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1910 ( -2.0*(sin(grid(i+1,j+1,2))**6.0) *(cos(grid(i+1,j+1,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1911 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1912 ( (8.0/5.0)*(cos(grid(i+1,j+1,2))**3.0)*(sin(grid(i+1,j+1,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1913 pt9 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1914 ( -2.0*(sin(grid(i,j+1,2))**6.0) *(cos(grid(i,j+1,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1915 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1916 ( (8.0/5.0)*(cos(grid(i,j+1,2))**3.0)*(sin(grid(i,j+1,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1917 pt(i,j,z) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1924 if ( (r/r0)**2 < 40. )
then 1925 pt(i,j,z) = pt(i,j,z) + pt0*exp(-(r/r0)**2)
1932 if (is_master()) print*,
' ' 1939 pt1 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1940 ( -2.0*(sin(agrid(i,j,2))**6.0) *(cos(agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1941 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1942 ( (8.0/5.0)*(cos(agrid(i,j,2))**3.0)*(sin(agrid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1953 pt2 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1954 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1955 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1956 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1958 pt3 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1959 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1960 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1961 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1963 pt4 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1964 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1965 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1966 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1968 pt5 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1969 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1970 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1971 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1973 pt6 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1974 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1975 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1976 ( (8.0/5.0)*(cos(grid(i,j,2))**3.0)*(sin(grid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1977 pt7 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1978 ( -2.0*(sin(grid(i+1,j,2))**6.0) *(cos(grid(i+1,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1979 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1980 ( (8.0/5.0)*(cos(grid(i+1,j,2))**3.0)*(sin(grid(i+1,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1981 pt8 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1982 ( -2.0*(sin(grid(i+1,j+1,2))**6.0) *(cos(grid(i+1,j+1,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1983 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1984 ( (8.0/5.0)*(cos(grid(i+1,j+1,2))**3.0)*(sin(grid(i+1,j+1,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1985 pt9 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1986 ( -2.0*(sin(grid(i,j+1,2))**6.0) *(cos(grid(i,j+1,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1987 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1988 ( (8.0/5.0)*(cos(grid(i,j+1,2))**3.0)*(sin(grid(i,j+1,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1989 phis(i,j) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1996 if ( .not.hydrostatic )
then 2002 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
2008 if (.not. adiabatic)
then 2009 zvir = rvgas/rdgas - 1.
2014 pt(i,j,k) = pt(i,j,k)/(1. + zvir*q(i,j,k,
sphum))
2024 write(stdout(), *)
'PI:', pi
2025 write(stdout(), *)
'PHIS:', mpp_chksum(phis(is:ie,js:je))
2030 is,ie,js,je,isd,ied,jsd,jed,npz,ncnst,ak,bk,ptop, &
2031 pk,peln,pe,pkz,gz,phis,ps,grid,agrid,hydrostatic, &
2032 nwat, adiabatic,
test_case == -13, domain)
2034 write(stdout(), *)
'PHIS:', mpp_chksum(phis(is:ie,js:je))
2060 ze1(k) = ze1(k+1) + ztop/
real(npz)
2064 ze1(1) = ztop + 1.5*ztop/
real(npz)
2077 delz(i,j,k) = ze1(k+1) - ze1(k)
2078 pk(i,j,k) = pk(i,j,k+1) + grav*delz(i,j,k)/(cp_air*t00)*pk0
2079 pe(i,k,j) = pk(i,j,k)**(1./kappa)
2085 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
2090 peln(i,k,j) = log(pe(i,k,j))
2099 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2100 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2113 r0 = 0.5*(ze1(k)+ze1(k+1)) - 3.2e3
2115 r0 = (0.5*(ze1(k)+ze1(k+1)) - 3.0e3) / 2.e3
2120 p2(1) = agrid(i,j,1)
2121 p2(2) = agrid(i,j,2)
2124 dist = sqrt( r**2 + r0**2 ) / 3.2e3
2127 dist = sqrt( r**2 + r0**2 )
2129 if ( dist<=1. )
then 2130 q(i,j,k,1) = pk0 * pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
2131 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
2136 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
2153 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2156 pe1(z) = ak(z) + bk(z)*p00
2161 ze1(z) = ze1(z+1) + ztop/
real(npz)
2165 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2169 ps(i,j) = pe1(npz+1)
2177 peln(i,z,j) = log(pe1(z))
2178 pk(i,j,z) = exp(kappa*peln(i,z,j))
2191 vort(i,j) = 0.5*(1.+cos(pi*r/r0))
2202 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*pi/ztop )
2205 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(kappa*(peln(i,z+1,j)-peln(i,z,j)))
2206 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2208 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2209 q(i,j,z,1) = q(i,j,z,1) + vort(i,j)*zmid
2222 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2225 pe1(z) = ak(z) + bk(z)*p00
2230 ze1(z) = ze1(z+1) + ztop/
real(npz)
2234 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2238 ps(i,j) = pe1(npz+1)
2246 peln(i,z,j) = log(pe1(z))
2247 pk(i,j,z) = exp(kappa*peln(i,z,j))
2260 vort(i,j) = 0.5*(1.+cos(pi*r/r0))
2270 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*pi/ztop )
2273 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(kappa*(peln(i,z+1,j)-peln(i,z,j)))
2274 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2276 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2284 n2 = grav**2 / (cp_air*pt0)
2293 phis(i,j) = grav*2.e3*exp( -(r/1500.e3)**2 )
2295 (sin(agrid(i,j,2))**2-1.) - n2/(grav*grav*kappa)*phis(i,j))
2303 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
2309 p1(:) = grid(i ,j ,1:2)
2310 p2(:) = grid(i,j+1 ,1:2)
2314 utmp =
ubar * cos(p3(2))
2323 p1(:) = grid(i, j,1:2)
2324 p2(:) = grid(i+1,j,1:2)
2328 utmp =
ubar * cos(p3(2))
2354 p1(:) = grid(i ,j ,1:2)
2355 p2(:) = grid(i,j+1 ,1:2)
2359 utmp =
ubar * cos(p3(2))
2366 p1(:) = grid(i, j,1:2)
2367 p2(:) = grid(i+1,j,1:2)
2371 utmp =
ubar * cos(p3(2))
2392 p1(1) = (0.5-0.125) * pi
2399 p2(:) = agrid(i,j,1:2)
2401 if ( r < pi*
radius )
then 2402 p4(:) = p2(:) - p1(:)
2403 if ( abs(p4(1)) > 1.e-12 )
then 2404 zeta = asin( p4(2) / sqrt(p4(1)**2 + p4(2)**2) )
2408 if ( p4(1) <= 0. ) zeta = pi - zeta
2410 v1 = r/uu1 * cos( zeta )
2411 v2 = r/uu2 * sin( zeta )
2412 phis(i,j) = ftop / ( 1. + v1**2 + v2**2 )
2419 if ( hybrid_z )
then 2423 elseif( npz.eq.31 .or. npz.eq.41 .or. npz.eq.51 )
then 2427 if ( is_master() )
write(*,*)
'Using const DZ' 2429 dz1(1) = ztop /
real(npz) 2434 dz1(1) = max( 1.0e3, 3.*dz1(2) )
2440 ze1(k) = ze1(k+1) + dz1(k)
2447 call mpp_error(fatal,
'This test case is only currently setup for hybrid_z')
2453 delz(i,j,k) = ze0(i,j,k+1) - ze0(i,j,k)
2463 s0 = grav*grav / (cp_air*n2)
2467 pe1(k) = p00*( (1.-s0/t00) + s0/t00*exp(-n2*ze1(k)/grav) )**(1./kappa)
2471 if ( is_master() )
write(*,*)
'Lee vortex testcase: model top (mb)=', ptop/100.
2477 bk(k) = (pe1(k) - pe1(1)) / (pe1(npz+1)-pe1(1))
2478 ak(k) = pe1(1)*(1.-bk(k))
2487 pk(i,j,k) = pk0 - (1.-exp(-n2/grav*ze0(i,j,k))) * (grav*grav)/(n2*cp_air*pt0)
2488 pe(i,k,j) = pk(i,j,k) ** (1./kappa)
2489 peln(i,k,j) = log(pe(i,k,j))
2497 peln(i,1,j) = log(pe(i,1,j))
2498 pk(i,j,1) = pe(i,1,j) ** kappa
2499 ps(i,j) = pe(i,npz+1,j)
2506 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2507 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2508 pt(i,j,k) = pkz(i,j,k)*grav*delz(i,j,k) / ( cp_air*(pk(i,j,k)-pk(i,j,k+1)) )
2519 if (.not.hydrostatic) w(:,:,:)= 0.0
2532 dz = 12000./
real(npz)
2534 allocate(
zz0(npz+1))
2535 allocate(
pz0(npz+1))
2543 if (is_master()) print*,
'TRACER ADVECTION TEST CASE' 2544 if (is_master()) print*,
'INITIAL LEVELS' 2550 if (is_master())
write(*,*) k,
pz0(k),
zz0(k)
2563 ptmp = 0.5*(
pz0(k) +
pz0(k+1))
2570 delp(i,j,k) =
pz0(k+1)-
pz0(k)
2575 ptop = 100000.*exp(-12000.*grav/t00/rdgas)
2583 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
2586 call mpp_update_domains( psi, domain )
2589 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
2590 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
2598 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
2599 if (dist==0) vc(i,j,k) = 0.
2605 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
2606 if (dist==0) uc(i,j,k) = 0.
2613 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
2614 if (dist==0) v(i,j,k) = 0.
2620 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
2621 if (dist==0) u(i,j,k) = 0.
2627 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
2628 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
2630 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
2631 if (dist==0) ua(i,j,k) = 0.
2632 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
2633 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
2635 va(i,j,k) = (psi2 - psi1) / (dist)
2636 if (dist==0) va(i,j,k) = 0.
2643 uc(:,:,k) = uc(:,:,1)
2644 vc(:,:,k) = vc(:,:,1)
2645 ua(:,:,k) = ua(:,:,1)
2646 va(:,:,k) = va(:,:,1)
2649 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
2650 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
2658 call mpp_error(fatal,
'Value of tracer_test not implemented ')
2676 if (.not.hydrostatic) w(:,:,:)= 0.0
2680 dz = 12000./
real(npz)
2685 exponent = rdgas*gamma/grav
2686 px = ((t00-9000.*gamma)/t00)**(1./exponent)
2690 height = 12000. - dz*
real(k-1)
2691 if (height >= 9000. )
then 2692 ak(k) = p00*((t00-height*gamma)/t00)**(1./exponent)
2695 ak(k) = (((t00-height*gamma)/t00)**(1./exponent)-1.)/(px - 1.)*px*p00
2696 bk(k) = (((t00-height*gamma)/t00)**(1./exponent)-px)/(1.-px)
2698 if (is_master())
write(*,*) k, ak(k), bk(k), height, ak(k)+bk(k)*p00
2704 p1(1) = 3.*pi/2. ; p1(2) = 0.
2711 p2(:) = agrid(i,j,1:2)
2714 phis(i,j) = grav*0.5*2000.*(1. + cos(pi*r/r0))*cos(pi*r/zetam)**2.
2715 pe(i,npz+1,j) = p00*(1.-gamma/t00*phis(i,j)/grav)**(1./exponent)
2720 ps(i,j) = pe(i,npz+1,j)
2727 pe(i,k,j) = ak(k) + bk(k)*ps(i,j)
2728 gz(i,j,k) = t00/gamma*(1. - (pe(i,k,j)/p00)**exponent)
2740 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
2745 pt(i,j,k) = -grav*t00*p00/(rdgas*gamma + grav)/delp(i,j,k) * &
2746 ( (pe(i,k,j)/p00)**(exponent+1.) - (pe(i,k+1,j)/p00)**(exponent+1.) )
2762 zvir = rvgas/rdgas - 1.
2768 pk(i,j,1) = ptop**kappa
2770 peln(i,1,j) = log(ptop)
2777 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
2778 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
2779 peln(i,k+1,j) = log(pe(i,k+1,j))
2780 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
2788 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2795 pp0(1) = 262.0/180.*pi
2796 pp0(2) = 35.0/180.*pi
2803 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
2810 ze1(k) = ze1(k+1) - delz(is,js,k)
2814 if (is_master())
then 2816 write(6,*)
'Toy supercell winds, piecewise approximation' 2818 write(6,*)
'Toy supercell winds, tanh approximation' 2823 zm = 0.5*(ze1(k)+ze1(k+1))
2828 if ( zm .le. 2.e3 )
then 2829 utmp = 8.*(1.-cos(pi*zm/4.e3))
2830 vtmp = 8.*sin(pi*zm/4.e3)
2831 elseif (zm .le. 6.e3 )
then 2832 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
2842 utmp = 15.0*(1.+tanh(zm/2000. - 1.5))
2843 vtmp = 8.5*tanh(zm/1000.)
2858 if( is_master() )
then 2859 write(6,*) k, utmp, vtmp
2864 p1(:) = grid(i ,j ,1:2)
2865 p2(:) = grid(i,j+1 ,1:2)
2876 p1(:) = grid(i, j,1:2)
2877 p2(:) = grid(i+1,j,1:2)
2888 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
2889 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
2890 .true., hydrostatic, nwat, domain)
2897 zm = 0.5*(ze1(k)+ze1(k+1))
2898 ptmp = ( (zm-zc)/zc ) **2
2899 if ( ptmp < 1. )
then 2903 if ( dist < 1. )
then 2904 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
2913 call mpp_error(fatal,
' test_case 32 not yet implemented')
2939 p0(1) = 60./180. * pi
2954 p2(2) = agrid(i,j,2)
2955 pt1 = exp( dum*(us0*sin(p2(2)))**2 )
2957 pt2 = exp( dum*(us0*sin(p2(2)))**2 )
2959 pt3 = exp( dum*(us0*sin(p2(2)))**2 )
2961 pt4 = exp( dum*(us0*sin(p2(2)))**2 )
2963 pt5 = exp( dum*(us0*sin(p2(2)))**2 )
2965 pt6 = exp( dum*(us0*sin(p2(2)))**2 )
2966 p2(2) = grid(i+1,j,2)
2967 pt7 = exp( dum*(us0*sin(p2(2)))**2 )
2968 p2(2) = grid(i+1,j+1,2)
2969 pt8 = exp( dum*(us0*sin(p2(2)))**2 )
2970 p2(2) = grid(i,j+1,2)
2971 pt9 = exp( dum*(us0*sin(p2(2)))**2 )
2972 ptmp = t00*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
2974 ptmp = t00*exp( dum*(us0*sin(agrid(i,j,2)))**2 )
2990 p2(1:2) = agrid(i,j,1:2)
2992 pt1 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
2995 pt2 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
2998 pt3 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3001 pt4 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3004 pt5 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3005 p2(1:2) = grid(i,j,1:2)
3007 pt6 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3008 p2(1:2) = grid(i+1,j,1:2)
3010 pt7 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3011 p2(1:2) = grid(i+1,j+1,1:2)
3013 pt8 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3014 p2(1:2) = grid(i,j+1,1:2)
3016 pt9 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3017 phis(i,j) = grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
3019 p2(1:2) = agrid(i,j,1:2)
3021 phis(i,j) = grav*h0*cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3038 pt1 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3041 pt2 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3044 pt3 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3047 pt4 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3050 pt5 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3052 pt6 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3054 pt7 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3056 pt8 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3058 pt9 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3059 phis(i,j) = grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
3062 pt1 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3063 phis(i,j) = grav*h0*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3072 ps(i,j) = p00*exp( -0.5*(us0*sin(agrid(i,j,2)))**2/(rdgas*t00)-phis(i,j)/(rdgas*pt(i,j,1)) )
3074 peln(i,1,j) = log(ptop)
3075 pk(i,j,1) = ptop**kappa
3082 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
3083 peln(i,k,j) = log(pe(i,k,j))
3084 pk(i,j,k) = exp( kappa*peln(i,k,j) )
3092 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3093 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
3101 ze1(npz+1) = phis(i,j)/grav
3103 ze1(k) = ze1(k+1) - delz(i,j,k)
3106 w(i,j,k) = 0.5*(ze1(k)+ze1(k+1))
3110 call mpp_update_domains( w, domain )
3115 p1(:) = grid(i ,j, 1:2)
3116 p2(:) = grid(i,j+1, 1:2)
3121 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i-1,j,k)+w(i,j,k)) )
3127 p1(:) = grid(i, j, 1:2)
3128 p2(:) = grid(i+1,j, 1:2)
3132 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i,j-1,k)+w(i,j,k)) )
3141 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3142 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
3143 .true., hydrostatic, nwat, domain)
3155 zvir = rvgas/rdgas - 1.
3168 ze1(k) = ze1(k+1) + ztop/
real(npz)
3177 call var_dz(npz, ztop, ze1)
3180 zs1(k) = 0.5*(ze1(k)+ze1(k+1))
3188 ts1(k) = cp_air*ts1(k)*(1.+zvir*qs1(k))
3194 call balanced_k(npz, is, ie, js, je,
ng, pe1(npz+1), ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
3195 delz, zvir, ptop, ak, bk, agrid)
3198 ps(i,j) = pe(i,npz+1,j)
3205 peln(i,k,j) = log(pe(i,k,j))
3206 pk(i,j,k) = exp( kappa*peln(i,k,j) )
3214 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3215 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3224 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
3231 delz(i,j,k) = ze1(k+1) - ze1(k)
3232 pt(i,j,k) = delz(i,j,k)*grav/(rdgas*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j)))
3241 p1(:) = grid(i ,j ,1:2)
3242 p2(:) = grid(i,j+1 ,1:2)
3246 v(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e2,ex)
3251 p1(:) = grid(i, j,1:2)
3252 p2(:) = grid(i+1,j,1:2)
3256 u(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e1,ex)
3273 zm = 0.5*(ze1(k)+ze1(k+1))
3274 ptmp = ( (zm-zc)/zc ) **2
3275 if ( ptmp < 1. )
then 3280 if ( dist < 1. )
then 3281 pt(i,j,k) = pt(i,j,k) + (pkz(i,j,k)/pk0)*pturb*cos(0.5*pi*dist)**2
3302 zvir = rvgas/rdgas - 1.
3309 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3317 peln(i,1,j) = log(pe(i,1,j))
3318 pk(i,j,1) = exp(kappa*peln(i,1,j))
3322 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3323 peln(i,k,j) = log(pe(i,k,j))
3324 pk(i,j,k) = exp(kappa*peln(i,k,j))
3336 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3338 if ( dist .le. r0 )
then 3350 if (.not.hydrostatic)
then 3354 delz(i,j,k) = rdgas*pt(i,j,k)*(1.+zvir*q(i,j,k,1))/grav*log(pe(i,k,j)/pe(i,k+1,j))
3376 zvir = rvgas/rdgas - 1.
3383 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
3391 peln(i,1,j) = log(pe(i,1,j))
3392 pk(i,j,1) = exp(kappa*peln(i,1,j))
3396 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3397 peln(i,k,j) = log(pe(i,k,j))
3398 pk(i,j,k) = exp(kappa*peln(i,k,j))
3419 p1(:) = grid(i ,j ,1:2)
3420 p2(:) = grid(i,j+1 ,1:2)
3423 utmp =
ubar*exp(-(r/r0)**2)
3431 p1(:) = grid(i, j,1:2)
3432 p2(:) = grid(i+1,j,1:2)
3435 utmp =
ubar*exp(-(r/r0)**2)
3444 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3446 pt(i,j,k) = pt0/p00**kappa
3448 pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
3462 q(i,j,k,:) = agrid(i,j,1)*0.180/pi
3468 ncnst, npz, q, agrid(is:ie,js:je,1), agrid(is:ie,js:je,2), 9., 9.)
3471 if ( .not. hydrostatic )
then 3475 delz(i,j,k) = rdgas*pt(i,j,k)/grav*log(pe(i,k,j)/pe(i,k+1,j))
3492 p0(1) = 180. * pi / 180.
3493 p0(2) = 10. * pi / 180.
3508 p2(:) = agrid(i,j,1:2)
3510 ps(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3515 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3521 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3533 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3541 p2(:) = 0.5*(grid(i,j,1:2)+grid(i,j+1,1:2))
3543 ps_v(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3548 p2(:) = 0.5*(grid(i,j,1:2)+grid(i+1,j,1:2))
3550 ps_u(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3561 pe_v(i,k,j) = ak(k) + ps_v(i,j)*bk(k)
3571 pe_u(i,k,j) = ak(k) + ps_u(i,j)*bk(k)
3580 zvir = rvgas/rdgas - 1.
3583 p0 = (/ pi, pi/18. /)
3590 t00 = ts0*(1.+zvir*q00)
3591 exponent = rdgas*gamma/grav
3595 cor = 2.*omega*sin(p0(2))
3600 p1(:) = grid(i ,j ,1:2)
3601 p2(:) = grid(i,j+1 ,1:2)
3606 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3607 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3608 d = max(1.e-15,sqrt(d1**2+d2**2))
3613 ptmp = 0.5*(pe_v(i,k,j)+pe_v(i,k+1,j))
3614 height = (t00/gamma)*(1.-(ptmp/ps_v(i,j))**exponent)
3615 if (height > ztrop)
then 3618 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3619 - exppr*(r/rp)**exppr*rdgas*(t00-gamma*height) &
3620 /(exppz*height*rdgas*(t00-gamma*height)/(grav*zp**exppz) &
3621 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3633 p1(:) = grid(i, j,1:2)
3634 p2(:) = grid(i+1,j,1:2)
3639 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3640 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3641 d = max(1.e-15,sqrt(d1**2+d2**2))
3646 ptmp = 0.5*(pe_u(i,k,j)+pe_u(i,k+1,j))
3647 height = (t00/gamma)*(1.-(ptmp/ps_u(i,j))**exponent)
3648 if (height > ztrop)
then 3651 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3652 - exppr*(r/rp)**exppr*rdgas*(t00-gamma*height) &
3653 /(exppz*height*rdgas*(t00-gamma*height)/(grav*zp**exppz) &
3654 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3666 ttrop = t00 - gamma*ztrop
3675 ptmp = 0.5*(pe(i,k,j)+pe(i,k+1,j))
3676 height = (t00/gamma)*(1.-(ptmp/ps(i,j))**exponent)
3677 if (height > ztrop)
then 3681 q(i,j,k,1) = q00*exp(-height/zq1)*exp(-(height/zq2)**exppz)
3682 p2(:) = agrid(i,j,1:2)
3684 pt(i,j,k) = (t00-gamma*height)/(1.d0+zvir*q(i,j,k,1))/(1.d0+exppz*rdgas*(t00-gamma*height)*height &
3685 /(grav*zp**exppz*(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz))))
3694 ps(i,j) = pe(i,npz+1,j)
3698 if (.not.hydrostatic)
then 3702 delz(i,j,k) = rdgas*pt(i,j,k)*(1.+zvir*q(i,j,k,1))/grav*log(pe(i,k,j)/pe(i,k+1,j))
3709 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy,
ng)
3711 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3729 call dcmip16_tc (delp, pt, u, v, q, w, delz, &
3730 is, ie, js, je, isd, ied, jsd, jed, npz, ncnst, &
3731 ak, bk, ptop, pk, peln, pe, pkz, gz, phis, &
3732 ps, grid, agrid, hydrostatic, nwat, adiabatic)
3736 call mpp_error(fatal,
" test_case not defined" )
3740 call mpp_update_domains( phis, domain )
3742 ftop =
g_sum(domain, phis(is:ie,js:je), is, ie, js, je,
ng, area, 1)
3743 if(is_master())
write(*,*)
'mean terrain height (m)=', ftop/grav
3747 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3748 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., mountain, &
3749 moist_phys, hydrostatic, nwat, domain, .not.hydrostatic)
3752 #ifdef COLUMN_TRACER 3753 if( ncnst>1 ) q(:,:,:,2:ncnst) = 0.0
3761 p1(:) = grid(i ,j ,1:2)
3762 p2(:) = grid(i,j+1 ,1:2)
3768 if (-(r/r0)**2.0 > -40.0) q(i,j,z,1) = exp(-(r/r0)**2.0)
3804 nullify(cubed_sphere)
3809 nullify(have_south_pole)
3810 nullify(have_north_pole)
3819 subroutine get_vorticity(isc, iec, jsc, jec ,isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
3820 integer isd, ied, jsd, jed, npz
3821 integer isc, iec, jsc, jec
3822 real,
intent(in) :: u(isd:ied, jsd:jed+1, npz), v(isd:ied+1, jsd:jed, npz)
3823 real,
intent(out) :: vort(isc:iec, jsc:jec, npz)
3824 real,
intent(IN) :: dx(isd:ied,jsd:jed+1)
3825 real,
intent(IN) :: dy(isd:ied+1,jsd:jed)
3826 real,
intent(IN) :: rarea(isd:ied,jsd:jed)
3828 real :: utmp(isc:iec, jsc:jec+1), vtmp(isc:iec+1, jsc:jec)
3834 utmp(i,j) = u(i,j,k)*dx(i,j)
3839 vtmp(i,j) = v(i,j,k)*dy(i,j)
3845 vort(i,j,k) = rarea(i,j)*(utmp(i,j)-utmp(i,j+1)-vtmp(i,j)+vtmp(i+1,j))
3852 subroutine checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, &
3853 nq, km, q, lon, lat, nx, ny, rn)
3863 integer,
intent(in):: nq
3864 integer,
intent(in):: km
3865 integer,
intent(in):: i0, i1
3866 integer,
intent(in):: j0, j1
3867 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3868 real,
intent(in):: nx
3869 real,
intent(in):: ny
3870 real,
intent(in),
optional:: rn
3871 real(kind=R_GRID),
intent(in),
dimension(i0:i1,j0:j1):: lon, lat
3872 real,
intent(out):: q(ifirst:ilast,jfirst:jlast,km,nq)
3874 real:: qt(i0:i1,j0:j1)
3882 qtmp = sin(nx*lon(i,j))*sin(ny*lat(i,j))
3883 if ( qtmp < 0. )
then 3891 if (
present(rn) )
then 3899 call random_number(ftmp)
3900 q(i,j,k,iq) = qt(i,j) + rn*ftmp
3912 q(i,j,k,iq) = qt(i,j)
3922 km, q, delp, ncnst, lon, lat)
3927 integer,
intent(in):: km
3928 integer,
intent(in):: i0, i1
3929 integer,
intent(in):: j0, j1
3930 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3931 integer,
intent(in):: ncnst
3932 real(kind=R_GRID),
intent(in),
dimension(ifirst:ilast,jfirst:jlast):: lon, lat
3933 real,
intent(inout):: q(ifirst:ilast,jfirst:jlast,km,ncnst)
3934 real,
intent(in):: delp(ifirst:ilast,jfirst:jlast,km)
3936 real:: D, k1, r, ll, sinthc, costhc, mm
3942 real,
parameter :: qcly = 4.e-6
3943 real,
parameter :: lc = 5.*pi/3.
3944 real,
parameter :: thc = pi/9.
3945 real,
parameter :: k2 = 1.
3950 cl = get_tracer_index(model_atmos,
'Cl')
3951 cl2 = get_tracer_index(model_atmos,
'Cl2')
3955 k1 = max(0., sin(lat(i,j))*sinthc + cos(lat(i,j))*costhc*cos(lon(i,j) - lc))
3957 d = sqrt(r*r + 2.*r*qcly)
3959 q(i,j,1,cl2) = 0.5*(qcly - q(i,j,1,cl))
3966 q(i,j,k,cl) = q(i,j,1,cl)
3967 q(i,j,k,cl2) = q(i,j,1,cl2)
3974 if (is_master())
then 3979 qcly0 =
qcly0 + (q(i,j,k,cl) + 2.*q(i,j,k,cl2))*delp(i,j,k)
3980 mm = mm + delp(i,j,k)
3985 if (is_master()) print*,
' qcly0 = ',
qcly0 3994 real,
intent(in):: ubar
3995 real,
intent(in):: r0
3996 real,
intent(in):: p1(2)
3997 real,
intent(inout):: u(isd:ied, jsd:jed+1)
3998 real,
intent(inout):: v(isd:ied+1,jsd:jed)
3999 real(kind=R_GRID),
intent(IN) :: grid(isd:ied+1,jsd:jed+1,2)
4001 real(kind=R_GRID):: p2(2), p3(2), p4(2)
4002 real(kind=R_GRID):: e1(3), e2(3), ex(3), ey(3)
4003 real:: vr, r, d2, cos_p, x1, y1
4012 p2(1) = p2(1) - p1(1)
4013 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
4021 x1 = cos(p2(2))*sin(p2(1))
4022 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
4023 d2 = max(1.e-25, sqrt(x1**2 + y1**2))
4026 p3(1) = grid(i,j, 1) - p1(1)
4027 p3(2) = grid(i,j, 2)
4028 p4(1) = grid(i+1,j,1) - p1(1)
4029 p4(2) = grid(i+1,j,2)
4041 p2(1) = p2(1) - p1(1)
4042 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
4049 x1 = cos(p2(2))*sin(p2(1))
4050 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
4051 d2 = max(1.e-25, sqrt(x1**2 + y1**2))
4054 p3(1) = grid(i,j, 1) - p1(1)
4055 p3(2) = grid(i,j, 2)
4056 p4(1) = grid(i,j+1,1) - p1(1)
4057 p4(2) = grid(i,j+1,2)
4067 real function gh_jet(npy, lat_in)
4068 integer,
intent(in):: npy
4069 real,
intent(in):: lat_in
4070 real lat, lon, dp, uu
4077 dp = pi /
real(jm-1)
4087 lat = -pi/2. + (
real(j-1)-0.5)*dp
4089 ft = 2.*omega*sin(lat)
4114 real function u_jet(lat)
4116 real umax, en, ph0, ph1
4121 en = exp( -4./(ph1-ph0)**2 )
4123 if ( lat>ph0 .and. lat<ph1 )
then 4124 u_jet = (umax/en)*exp( 1./( (lat-ph0)*(lat-ph1) ) )
4131 real,
intent(OUT) :: B(isd:ied,jsd:jed)
4132 real,
intent(IN) :: agrid(isd:ied,jsd:jed,2)
4140 if (sin(agrid(i,j,2)) > 0.)
then 4141 myc = sin(agrid(i,j,1))
4142 yy = (cos(agrid(i,j,2))/sin(agrid(i,j,2)))**2
4143 myb =
gh0*yy*exp(1.-yy)
4161 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4162 real ,
intent(IN) :: time_since_start
4168 tday = time_since_start/86400.0
4169 if (tday >= 20.)
then 4170 aoft(2) = 0.5*(1.-cos(0.25*pi*(tday-20)))
4171 if (tday == 24)
aoft(2) = 1.0
4172 elseif (tday <= 4.)
then 4173 aoft(2) = 0.5*(1.-cos(0.25*pi*tday))
4174 elseif (tday <= 16.)
then 4177 aoft(2) = 0.5*(1.+cos(0.25*pi*(tday-16.)))
4182 phis(i,j) = amean*
case9_b(i,j)
4195 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4212 subroutine case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain)
4214 real,
intent(INOUT) :: delp(isd:ied,jsd:jed,npz)
4215 real,
intent(INOUT) :: uc(isd:ied+1,jsd:jed,npz)
4216 real,
intent(INOUT) :: vc(isd:ied,jsd:jed+1,npz)
4217 real,
intent(INOUT) :: u(isd:ied,jsd:jed+1,npz)
4218 real,
intent(INOUT) :: v(isd:ied+1,jsd:jed,npz)
4219 real,
intent(INOUT) :: ua(isd:ied,jsd:jed,npz)
4220 real,
intent(INOUT) :: va(isd:ied,jsd:jed,npz)
4221 real,
intent(INOUT) :: pe(is-1:ie+1, npz+1,js-1:je+1)
4222 real,
intent(IN) :: time, dt
4223 real,
intent(INOUT) :: ptop
4224 integer,
intent(IN) :: npx, npy, npz
4226 type(domain2d),
intent(INOUT) :: domain
4233 real :: s, l, dt2, v0, phase
4234 real :: ull, vll, lonp
4235 real :: p0(2), elon(3), elat(3)
4237 real :: psi(isd:ied,jsd:jed)
4238 real :: psi_b(isd:ied+1,jsd:jed+1)
4239 real :: dist, psi1, psi2
4244 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3), pt(2), p1(2), p2(2), p3(2), rperiod, timefac, t00
4248 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
4249 real,
pointer,
dimension(:,:) :: dx, dxa, dy, dya, dxc, dyc
4251 agrid => gridstruct%agrid_64
4252 grid => gridstruct%grid_64
4255 dxa => gridstruct%dxa
4256 dxc => gridstruct%dxc
4258 dya => gridstruct%dya
4259 dyc => gridstruct%dyc
4261 period =
real( 12*24*3600 ) 4266 phase = pi*time/period
4279 omega0 = 23000.*pi/period
4282 ptop = 100000.*exp(-12000.*grav/t00/rdgas)
4287 s = min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*pi)))
4288 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(agrid(i,j,1)-period*(time+dt2))*cos(agrid(i,j,2))* &
4289 cos(period*(time+dt2))*sin(s*0.5*pi)
4297 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4308 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
4311 call mpp_update_domains( psi, domain )
4314 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
4315 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
4324 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
4325 if (dist==0) vc(i,j,k) = 0.
4331 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
4332 if (dist==0) uc(i,j,k) = 0.
4339 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
4340 if (dist==0) v(i,j,k) = 0.
4346 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
4347 if (dist==0) u(i,j,k) = 0.
4353 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
4354 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
4356 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
4357 if (dist==0) ua(i,j,k) = 0.
4358 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
4359 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
4361 va(i,j,k) = (psi2 - psi1) / (dist)
4362 if (dist==0) va(i,j,k) = 0.
4368 omega0 = 23000.*pi/period
4373 s = min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*pi)))
4374 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(agrid(i,j,1)-period*(time+dt2))*cos(agrid(i,j,2))* &
4375 cos(period*(time+dt2))*sin(s*0.5*pi)
4383 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4393 p1(:) = grid(i ,j ,1:2)
4394 p2(:) = grid(i,j+1 ,1:2)
4398 l = p3(1) - 2.*pi*time/period
4399 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(pi*time/period) + 2.*pi*
radius/period*cos(p3(2))
4400 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(pi*time/period)
4406 p1(:) = grid(i, j,1:2)
4407 p2(:) = grid(i+1,j,1:2)
4411 l = p3(1) - 2.*pi*time/period
4412 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(pi*time/period) + 2.*pi*
radius/period*cos(p3(2))
4413 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(pi*time/period)
4436 call dtoa( u(:,:,1), v(:,:,1),ua(:,:,1),va(:,:,1),dx,dy,dxa,dya,dxc,dyc,npx,npy,
ng)
4437 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
4438 call atoc(ua(:,:,1),va(:,:,1),uc(:,:,1),vc(:,:,1),dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
4443 ua(i,j,k) = ua(i,j,1)
4448 va(i,j,k) = va(i,j,1)
4456 vc(i,j,k) = vc(i,j,1)
4461 uc(i,j,k) = uc(i,j,1)
4474 pe(i,k,j) = pe(i,k,j) + dt*omega0*grav*pe(i,k,j)/rdgas/300./k_cell* &
4475 (-2.*sin(k_cell*agrid(i,j,2))*sin(agrid(i,j,2)) + k_cell*cos(agrid(i,j,2))*cos(k_cell*agrid(i,j,2)))* &
4476 sin(pi*
zz0(k)/12000.)*cos(phase)
4484 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4495 utmp =
ubar*cos(agrid(i,j,2))
4496 vtmp = -
radius * omega0 * pi / k_cell / 12000. * &
4497 cos(agrid(i,j,2)) * sin(k_cell * agrid(i,j,2)) * &
4498 sin(pi*
zz0(k)/12000.)*cos(phase)
4507 uc(:,:,k) = uc(:,:,1)
4508 vc(:,:,k) = vc(:,:,1)
4509 ua(:,:,k) = ua(:,:,1)
4510 va(:,:,k) = va(:,:,1)
4513 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
4514 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
4531 subroutine get_stats(dt, dtout, nt, maxnt, ndays, u,v,pt,delp,q,phis, ps, &
4532 uc,vc, ua,va, npx, npy, npz, ncnst, ndims, nregions, &
4533 gridstruct, stats_lun, consv_lun, monitorFreq, tile, &
4535 integer,
intent(IN) :: nt, maxnt
4536 real ,
intent(IN) :: dt, dtout, ndays
4537 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
4538 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
4539 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
4540 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
4541 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
4542 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4543 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
4544 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
4545 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
4546 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
4547 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
4548 integer,
intent(IN) :: npx, npy, npz, ncnst, tile
4549 integer,
intent(IN) :: ndims
4550 integer,
intent(IN) :: nregions
4551 integer,
intent(IN) :: stats_lun
4552 integer,
intent(IN) :: consv_lun
4553 integer,
intent(IN) :: monitorfreq
4555 type(domain2d),
intent(INOUT) :: domain
4556 logical,
intent(IN) :: nested
4561 real :: pmin, pmin1, uamin1, vamin1
4562 real :: pmax, pmax1, uamax1, vamax1
4563 real(kind=4) :: arr_r4(5)
4564 real :: tmass0, tvort0, tener0, tke0
4565 real :: tmass, tvort, tener, tke
4566 real :: temp(is:ie,js:je)
4567 integer :: i0, j0, k0, n0
4568 integer :: i, j, k, n, iq
4570 real :: psmo, vtx, p, w_p, p0
4571 real :: x1,y1,z1,x2,y2,z2,ang
4573 real :: p1(2), p2(2), p3(2), r, r0, dist, heading
4575 real :: uc0(isd:ied+1,jsd:jed ,npz)
4576 real :: vc0(isd:ied ,jsd:jed+1,npz)
4581 real,
save,
allocatable,
dimension(:,:,:) :: u0, v0
4582 real :: up(isd:ied ,jsd:jed+1,npz)
4583 real :: vp(isd:ied+1,jsd:jed ,npz)
4585 real,
dimension(:,:,:),
pointer :: grid, agrid
4586 real,
dimension(:,:),
pointer :: area, f0, dx, dy, dxa, dya, dxc, dyc
4588 grid => gridstruct%grid
4589 agrid=> gridstruct%agrid
4591 area => gridstruct%area
4596 dxa => gridstruct%dxa
4597 dya => gridstruct%dya
4598 dxc => gridstruct%dxc
4599 dyc => gridstruct%dyc
4602 if (nt == 0 .and. is_master()) print*,
'INITIALIZING GET_STATS' 4605 myday = ndays*((float(nt)/float(maxnt)))
4607 #if defined(SW_DYNAMICS) 4616 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
4618 if (p /= 0.0) w_p = vtx/p
4620 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
4630 dist = 2.0*pi*
radius* ((float(nt)/float(maxnt)))
4631 heading = 3.0*pi/2.0 -
alpha 4636 p2(1) = agrid(i,j,1)
4637 p2(2) = agrid(i,j,2)
4640 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
4642 phi0(i,j,1) = phis(i,j)
4649 call pmxn(delp(:,:,1), npx, npy, nregions, tile, gridstruct, pmin1, pmax1, i0, j0, n0)
4653 call get_scalar_stats( delp(:,:,1),
phi0(:,:,1), npx, npy, ndims, nregions, &
4654 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4661 arr_r4(5) = linf_norm
4674 200
format(i6.6,a,i6.6,a,e21.14)
4675 201
format(
' ',a,e21.14,
' ',e21.14)
4676 202
format(
' ',a,i4.4,
'x',i4.4,
'x',i4.4)
4678 if ( (is_master()) .and. mod(nt,monitorfreq)==0 )
then 4679 write(*,200) nt,
' step of ', maxnt,
' DAY ', myday
4680 write(*,201)
'Height MAX : ', pmax1
4681 write(*,201)
'Height MIN : ', pmin1
4682 write(*,202)
'HGT MAX location : ', i0, j0, n0
4684 write(*,201)
'Height L1_norm : ', l1_norm
4685 write(*,201)
'Height L2_norm : ', l2_norm
4686 write(*,201)
'Height Linf_norm : ', linf_norm
4691 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy,
ng)
4692 call pmxn(ua(:,:,1), npx, npy, nregions, tile, gridstruct, pmin1, pmax1, i0, j0, n0)
4694 call get_vector_stats( ua(:,:,1),
ua0(:,:,1), va(:,:,1),
va0(:,:,1), npx, npy, ndims, nregions, &
4695 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4701 arr_r4(5) = linf_norm
4703 if ( (is_master()) .and. mod(nt,monitorfreq)==0)
then 4704 write(*,201)
'UV MAX : ', pmax1
4705 write(*,201)
'UV MIN : ', pmin1
4706 write(*,202)
'UV MAX location : ', i0, j0, n0
4708 write(*,201)
'UV L1_norm : ', l1_norm
4709 write(*,201)
'UV L2_norm : ', l2_norm
4710 write(*,201)
'UV Linf_norm : ', linf_norm
4715 200
format(i6.6,a,i6.6,a,e10.4)
4716 201
format(
' ',a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4)
4717 202
format(
' ',a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4,
' ',e10.4)
4718 203
format(
' ',a,i3.3,a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4)
4720 if(is_master())
write(*,200) nt,
' step of ', maxnt,
' DAY ', myday
4723 psmo =
globalsum(ps(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4724 if(is_master())
write(*,*)
' Total surface pressure =', 0.01*psmo
4725 call pmxn(ps, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4726 if (is_master())
then 4727 write(*,201)
'PS MAX|MIN : ', 0.01*pmax, 0.01*pmin, i0, j0, n0
4738 call pmxn(pt(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4739 pmin1 = min(pmin, pmin1)
4740 pmax1 = max(pmax, pmax1)
4741 if (pmax1 == pmax) k0 = k
4743 if (is_master())
then 4744 write(*,201)
'PT MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4747 #if defined(DEBUG_TEST_CASES) 4748 if(is_master())
write(*,*)
' ' 4756 call pmxn(pt(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4757 pmin1 = min(pmin, pmin1)
4758 pmax1 = max(pmax, pmax1)
4759 if (is_master())
then 4760 write(*,202)
'PT MAX|MIN : ', pmax1, pmin1, i0, j0, k, n0, 0.5*( (ak(k)+ak(k+1))/1.e5 + bk(k)+bk(k+1) )
4763 if(is_master())
write(*,*)
' ' 4774 call pmxn(delp(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4775 pmin1 = min(pmin, pmin1)
4776 pmax1 = max(pmax, pmax1)
4777 if (pmax1 == pmax) k0 = k
4779 if (is_master())
then 4780 write(*,201)
'Delp MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4791 call dtoa(u(isd,jsd,k), v(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), dx,dy,dxa,dya,dxc,dyc,npx, npy,
ng)
4792 call pmxn(ua(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4793 uamin1 = min(pmin, uamin1)
4794 uamax1 = max(pmax, uamax1)
4795 if (uamax1 == pmax) k0 = k
4797 if (is_master())
then 4798 write(*,201)
'U MAX|MIN : ', uamax1, uamin1, i0, j0, k0, n0
4808 call pmxn(va(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4809 vamin1 = min(pmin, vamin1)
4810 vamax1 = max(pmax, vamax1)
4811 if (vamax1 == pmax) k0 = k
4813 if (is_master())
then 4814 write(*,201)
'V MAX|MIN : ', vamax1, vamin1, i0, j0, k0, n0
4825 call pmxn(q(isd,jsd,k,1), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4826 pmin1 = min(pmin, pmin1)
4827 pmax1 = max(pmax, pmax1)
4828 if (pmax1 == pmax) k0 = k
4830 if (is_master())
then 4831 write(*,201)
'Q MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4843 call pmxn(q(isd,jsd,k,iq), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4844 pmin1 = min(pmin, pmin1)
4845 pmax1 = max(pmax, pmax1)
4846 if (pmax1 == pmax) k0 = k
4848 if (is_master())
then 4849 write(*,203)
'TR',iq-1,
' MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4857 call get_vector_stats( ua(:,:,22),
ua0(:,:,22), va(:,:,22),
va0(:,:,22), npx, npy, ndims, nregions, &
4858 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4859 if (is_master())
then 4860 write(*,201)
'UV(850) L1_norm : ', l1_norm
4861 write(*,201)
'UV(850) L2_norm : ', l2_norm
4862 write(*,201)
'UV(850) Linf_norm : ', linf_norm
4870 #if defined(SW_DYNAMICS) 4878 temp(:,:) = delp(is:ie,js:je,k)
4879 tmass0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4880 tmass = tmass + tmass0
4884 allocate(u0(isd:ied,jsd:jed+1,npz))
4885 allocate(v0(isd:ied+1,jsd:jed,npz))
4894 call dtoa(up(isd,jsd,k), vp(isd,jsd,k), ua, va, dx,dy, dxa, dya, dxc, dyc, npx, npy,
ng)
4895 call atoc(ua(isd,jsd,k),va(isd,jsd,k),uc0(isd,jsd,k),vc0(isd,jsd,k),dx,dy,dxa,dya,npx,npy,
ng,nested, domain, nocomm=.true.)
4899 temp(i,j) = ( uc0(i,j,k)*uc0(i,j,k) + uc0(i+1,j,k)*uc0(i+1,j,k) + &
4900 vc0(i,j,k)*vc0(i,j,k) + vc0(i,j+1,k)*vc0(i,j+1,k) )
4903 tke0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4909 temp(i,j) = 0.5 * (delp(i,j,k)/grav) * temp(i,j)
4910 temp(i,j) = temp(i,j) + &
4911 grav*((delp(i,j,k)/grav + phis(i,j))*(delp(i,j,k)/grav + phis(i,j))) - &
4915 tener0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4916 tener = tener + tener0
4922 temp(i,j) = f0(i,j) + (1./area(i,j)) * ( (v(i+1,j,k)*dy(i+1,j) - v(i,j,k)*dy(i,j)) - &
4923 (u(i,j+1,k)*dx(i,j+1) - u(i,j,k)*dx(i,j)) )
4924 temp(i,j) = ( grav*(temp(i,j)*temp(i,j))/delp(i,j,k) )
4927 tvort0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4928 tvort = tvort + tvort0
4944 #if defined(SW_DYNAMICS) 4947 myrec = myday*86400.0/dtout + 1
4949 if (is_master())
write(consv_lun,rec=myrec) arr_r4(1:4)
4950 #if defined(SW_DYNAMICS) 4951 if ( (is_master()) .and. mod(nt,monitorfreq)==0)
then 4953 if ( (is_master()) )
then 4955 write(*,201)
'MASS TOTAL : ', tmass
4958 write(*,201)
'Kinetic Energy KE : ', tke
4959 write(*,201)
'ENERGY TOTAL : ', tener
4961 write(*,201)
'ENSTR TOTAL : ', tvort
4986 real ,
intent(IN) :: p1(2), p2(2)
4987 real ,
intent(IN) :: dist
4988 real ,
intent(IN) :: heading
4989 real ,
intent(OUT) :: p3(2)
4995 p3(2) = asin( (cos(heading)*cos(p1(2))*sin(pha)) + (sin(p1(2))*cos(pha)) )
4996 dp = atan2( sin(heading)*sin(pha)*cos(p1(2)) , cos(pha) - sin(p1(2))*sin(p3(2)) )
4997 p3(1) = mod( (p1(1)-pi)-dp+pi , 2.*pi )
5014 vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
5015 integer,
intent(IN) :: npx, npy
5016 integer,
intent(IN) :: ndims
5017 integer,
intent(IN) :: nregions, tile
5018 real ,
intent(IN) :: var(isd:ied,jsd:jed)
5019 real ,
intent(IN) :: varT(isd:ied,jsd:jed)
5020 real ,
intent(OUT) :: vmin
5021 real ,
intent(OUT) :: vmax
5022 real ,
intent(OUT) :: L1_norm
5023 real ,
intent(OUT) :: L2_norm
5024 real ,
intent(OUT) :: Linf_norm
5026 type(fv_grid_type),
target :: gridstruct
5035 real :: varSUM, varSUM2, varMAX
5037 real :: vminT, vmaxT, vmeanT, vvarT
5038 integer :: i0, j0, n0
5040 real,
dimension(:,:,:),
pointer :: grid, agrid
5041 real,
dimension(:,:),
pointer :: area
5043 grid => gridstruct%grid
5044 agrid=> gridstruct%agrid
5046 area => gridstruct%area
5065 vmean =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5066 vmeant =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5067 vmean = vmean / (4.0*pi)
5068 vmeant = vmeant / (4.0*pi)
5070 call pmxn(var, npx, npy, nregions, tile, gridstruct, vmin , vmax , i0, j0, n0)
5071 call pmxn(vart, npx, npy, nregions, tile, gridstruct, vmint, vmaxt, i0, j0, n0)
5072 call pmxn(var-vart, npx, npy, nregions, tile, gridstruct, pdiffmn, pdiffmx, i0, j0, n0)
5074 vmax = (vmax - vmaxt) / (vmaxt-vmint)
5075 vmin = (vmin - vmint) / (vmaxt-vmint)
5077 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5078 varsum2 =
globalsum(vart(is:ie,js:je)**2., npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5079 l1_norm =
globalsum(abs(var(is:ie,js:je)-vart(is:ie,js:je)), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5080 l2_norm =
globalsum((var(is:ie,js:je)-vart(is:ie,js:je))**2., npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5081 l1_norm = l1_norm/varsum
5082 l2_norm = sqrt(l2_norm)/sqrt(varsum2)
5084 call pmxn(abs(vart), npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5086 call pmxn(abs(var-vart), npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5087 linf_norm = vmax/varmax
5102 npx, npy, ndims, nregions, &
5103 vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
5104 integer,
intent(IN) :: npx, npy
5105 integer,
intent(IN) :: ndims
5106 integer,
intent(IN) :: nregions, tile
5107 real ,
intent(IN) :: varU(isd:ied,jsd:jed)
5108 real ,
intent(IN) :: varUT(isd:ied,jsd:jed)
5109 real ,
intent(IN) :: varV(isd:ied,jsd:jed)
5110 real ,
intent(IN) :: varVT(isd:ied,jsd:jed)
5111 real ,
intent(OUT) :: vmin
5112 real ,
intent(OUT) :: vmax
5113 real ,
intent(OUT) :: L1_norm
5114 real ,
intent(OUT) :: L2_norm
5115 real ,
intent(OUT) :: Linf_norm
5117 real :: var(isd:ied,jsd:jed)
5118 real :: varT(isd:ied,jsd:jed)
5126 real :: varSUM, varSUM2, varMAX
5128 real :: vminT, vmaxT, vmeanT, vvarT
5130 integer :: i0, j0, n0
5132 type(fv_grid_type),
target :: gridstruct
5134 real,
dimension(:,:,:),
pointer :: grid, agrid
5135 real,
dimension(:,:),
pointer :: area
5137 grid => gridstruct%grid
5138 agrid=> gridstruct%agrid
5140 area => gridstruct%area
5161 var(i,j) = sqrt( (varu(i,j)-varut(i,j))**2. + &
5162 (varv(i,j)-varvt(i,j))**2. )
5163 vart(i,j) = sqrt( varut(i,j)*varut(i,j) + &
5164 varvt(i,j)*varvt(i,j) )
5167 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5168 l1_norm =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5169 l1_norm = l1_norm/varsum
5171 call pmxn(vart, npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5173 call pmxn(var, npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5174 linf_norm = vmax/varmax
5178 var(i,j) = ( (varu(i,j)-varut(i,j))**2. + &
5179 (varv(i,j)-varvt(i,j))**2. )
5180 vart(i,j) = ( varut(i,j)*varut(i,j) + &
5181 varvt(i,j)*varvt(i,j) )
5184 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5185 l2_norm =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5186 l2_norm = sqrt(l2_norm)/sqrt(varsum)
5200 real,
intent(IN) :: ndt
5201 integer,
intent(IN) :: n_split
5202 integer,
intent(IN) :: npx, npy, npz, tile
5203 logical,
OPTIONAL,
intent(IN) :: noprint
5204 real ,
intent(IN) :: uc(isd:ied+1,jsd:jed ,npz)
5205 real ,
intent(IN) :: vc(isd:ied ,jsd:jed+1,npz)
5207 real :: ideal_c=0.06
5208 real :: tolerance= 1.e-3
5209 real :: dt_inc, dt_orig
5210 real :: meancy, mincy, maxcy, meancx, mincx, maxcx
5219 real,
dimension(:,:),
pointer :: dxc, dyc
5221 dxc => gridstruct%dxc
5222 dyc => gridstruct%dyc
5224 dt = ndt/
real(n_split)
5226 300
format(i4.4,
' ',i4.4,
' ',i4.4,
' ',i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
5232 do while(.not. ideal)
5244 mincx = min(mincx, abs( (dt/dxc(i,j))*uc(i,j,k) ))
5245 maxcx = max(maxcx, abs( (dt/dxc(i,j))*uc(i,j,k) ))
5246 meancx = meancx + abs( (dt/dxc(i,j))*uc(i,j,k) )
5248 if (abs( (dt/dxc(i,j))*uc(i,j,k) ) > 1.0)
then 5250 write(*,300) i,j,k,tile, abs( (dt/dxc(i,j))*uc(i,j,k) ), dt, dxc(i,j), uc(i,j,k), counter
5258 mincy = min(mincy, abs( (dt/dyc(i,j))*vc(i,j,k) ))
5259 maxcy = max(maxcy, abs( (dt/dyc(i,j))*vc(i,j,k) ))
5260 meancy = meancy + abs( (dt/dyc(i,j))*vc(i,j,k) )
5262 if (abs( (dt/dyc(i,j))*vc(i,j,k) ) > 1.0)
then 5264 write(*,300) i,j,k,tile, abs( (dt/dyc(i,j))*vc(i,j,k) ), dt, dyc(i,j), vc(i,j,k), counter
5272 call mp_reduce_max(maxcx)
5273 call mp_reduce_max(maxcy)
5276 call mp_reduce_max(mincx)
5277 call mp_reduce_max(mincy)
5280 call mp_reduce_sum(meancx)
5281 call mp_reduce_sum(meancy)
5282 meancx = meancx/(6.0*dble(npx)*dble(npy-1))
5283 meancy = meancy/(6.0*dble(npx-1)*dble(npy))
5295 if ( (.not.
present(noprint)) .and. (is_master()) )
then 5297 print*,
'--------------------------------------------' 5298 print*,
'Y-dir Courant number MIN : ', mincy
5299 print*,
'Y-dir Courant number MAX : ', maxcy
5301 print*,
'X-dir Courant number MIN : ', mincx
5302 print*,
'X-dir Courant number MAX : ', maxcx
5304 print*,
'X-dir Courant number MEAN : ', meancx
5305 print*,
'Y-dir Courant number MEAN : ', meancy
5307 print*,
'NDT: ', ndt
5308 print*,
'n_split: ', n_split
5311 print*,
'--------------------------------------------' 5325 subroutine pmxn(p, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
5326 integer,
intent(IN) :: npx
5327 integer,
intent(IN) :: npy
5328 integer,
intent(IN) :: nregions, tile
5329 real ,
intent(IN) :: p(isd:ied,jsd:jed)
5330 type(fv_grid_type),
intent(IN),
target :: gridstruct
5331 real ,
intent(OUT) :: pmin
5332 real ,
intent(OUT) :: pmax
5333 integer,
intent(OUT) :: i0
5334 integer,
intent(OUT) :: j0
5335 integer,
intent(OUT) :: n0
5341 real,
pointer,
dimension(:,:,:) :: agrid, grid
5342 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
5343 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
5344 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
5345 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5347 logical,
pointer :: cubed_sphere, latlon
5349 logical,
pointer :: have_south_pole, have_north_pole
5351 integer,
pointer :: ntiles_g
5352 real,
pointer :: acapN, acapS, globalarea
5354 grid => gridstruct%grid
5355 agrid=> gridstruct%agrid
5357 area => gridstruct%area
5358 rarea => gridstruct%rarea
5363 ee1 => gridstruct%ee1
5364 ee2 => gridstruct%ee2
5367 en1 => gridstruct%en1
5368 en2 => gridstruct%en2
5372 dxa => gridstruct%dxa
5373 dya => gridstruct%dya
5374 rdxa => gridstruct%rdxa
5375 rdya => gridstruct%rdya
5376 dxc => gridstruct%dxc
5377 dyc => gridstruct%dyc
5379 cubed_sphere => gridstruct%cubed_sphere
5380 latlon => gridstruct%latlon
5382 have_south_pole => gridstruct%have_south_pole
5383 have_north_pole => gridstruct%have_north_pole
5385 ntiles_g => gridstruct%ntiles_g
5386 acapn => gridstruct%acapN
5387 acaps => gridstruct%acapS
5388 globalarea => gridstruct%globalarea
5399 if (temp > pmax)
then 5403 elseif (temp < pmin)
then 5410 call mp_reduce_max(temp)
5411 if (temp /= pmax)
then 5417 call mp_reduce_max(i0)
5418 call mp_reduce_max(j0)
5419 call mp_reduce_max(n0)
5422 call mp_reduce_max(pmin)
5438 subroutine output_ncdf(dt, nt, maxnt, nout, u,v,pt,delp,q,phis,ps, uc,vc, ua,va, &
5439 omga, npx, npy, npz, ng, ncnst, ndims, nregions, ncid, &
5440 npx_p1_id, npy_p1_id, npx_id, npy_id, npz_id, ntiles_id, ncnst_id, nt_id, &
5441 phis_id, delp_id, ps_id, pt_id, pv_id, om_id, u_id, v_id, q_id, tracers_ids, &
5442 lats_id, lons_id, gridstruct, flagstruct)
5443 real,
intent(IN) :: dt
5444 integer,
intent(IN) :: nt, maxnt
5445 integer,
intent(INOUT) :: nout
5447 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
5448 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
5449 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
5450 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
5451 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
5453 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
5454 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
5456 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
5457 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
5458 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
5459 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
5460 real ,
intent(INOUT) :: omga(isd:ied ,jsd:jed ,npz)
5462 integer,
intent(IN) :: npx, npy, npz
5463 integer,
intent(IN) :: ng, ncnst
5464 integer,
intent(IN) :: ndims
5465 integer,
intent(IN) :: nregions
5466 integer,
intent(IN) :: ncid
5467 integer,
intent(IN) :: npx_p1_id, npy_p1_id, npx_id, npy_id, npz_id, ncnst_id
5468 integer,
intent(IN) :: ntiles_id, nt_id
5469 integer,
intent(IN) :: phis_id, delp_id, ps_id, pt_id, pv_id, u_id, v_id, q_id
5470 integer,
intent(IN) :: om_id
5471 integer,
intent(IN) :: tracers_ids(ncnst-1)
5472 integer,
intent(IN) :: lats_id, lons_id
5474 type(fv_grid_type),
target :: gridstruct
5475 type(fv_flags_type),
intent(IN) :: flagstruct
5477 real,
allocatable :: tmp(:,:,:)
5478 real,
allocatable :: tmpA(:,:,:)
5479 #if defined(SW_DYNAMICS) 5480 real,
allocatable :: ut(:,:,:)
5481 real,
allocatable :: vt(:,:,:)
5483 real,
allocatable :: ut(:,:,:,:)
5484 real,
allocatable :: vt(:,:,:,:)
5485 real,
allocatable :: tmpA_3d(:,:,:,:)
5487 real,
allocatable :: vort(:,:)
5494 real :: utmp, vtmp, r, r0, dist, heading
5495 integer :: i,j,k,n,iq,nreg
5498 real :: x1,y1,z1,x2,y2,z2,ang
5500 real,
pointer,
dimension(:,:,:) :: agrid, grid
5501 real,
pointer,
dimension(:,:) :: area, rarea
5502 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5504 grid => gridstruct%grid
5505 agrid => gridstruct%agrid
5507 area => gridstruct%area
5508 rarea => gridstruct%rarea
5512 dxa => gridstruct%dxa
5513 dya => gridstruct%dya
5514 rdxa => gridstruct%rdxa
5515 rdya => gridstruct%rdya
5516 dxc => gridstruct%dxc
5517 dyc => gridstruct%dyc
5519 allocate( tmp(npx ,npy ,nregions) )
5520 allocate( tmpa(npx-1,npy-1,nregions) )
5521 #if defined(SW_DYNAMICS) 5522 allocate( ut(npx-1,npy-1,nregions) )
5523 allocate( vt(npx-1,npy-1,nregions) )
5525 allocate( ut(npx-1,npy-1,npz,nregions) )
5526 allocate( vt(npx-1,npy-1,npz,nregions) )
5527 allocate( tmpa_3d(npx-1,npy-1,npz,nregions) )
5529 allocate( vort(isd:ied,jsd:jed) )
5534 tmp(is:ie+1,js:je+1,tile) = grid(is:ie+1,js:je+1,2)
5535 call wrtvar_ncdf(ncid, lats_id, nout, is,ie+1, js,je+1, npx+1, npy+1, 1, nregions, tmp(1:npx,1:npy,1:nregions), 3)
5536 tmp(is:ie+1,js:je+1,tile) = grid(is:ie+1,js:je+1,1)
5537 call wrtvar_ncdf(ncid, lons_id, nout, is,ie+1, js,je+1, npx+1, npy+1, 1, nregions, tmp(1:npx,1:npy,1:nregions), 3)
5540 #if defined(SW_DYNAMICS) 5542 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)/grav
5551 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
5552 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0) / grav
5566 dist = 2.0*pi*
radius* ((float(nt)/float(maxnt)))
5567 heading = 5.0*pi/2.0 -
alpha 5571 p2(1) = agrid(i,j,1)
5572 p2(2) = agrid(i,j,2)
5575 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
5577 phi0(i,j,1) = phis(i,j)
5589 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
5591 if (p /= 0.0) w_p = vtx/p
5592 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
5597 tmpa(is:ie,js:je,tile) =
phi0(is:ie,js:je,1)
5598 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5599 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)
5601 call wrtvar_ncdf(ncid, ps_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5607 vort(i,j) = f0(i,j) + (1./area(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
5608 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
5609 vort(i,j) = grav*vort(i,j)/delp(i,j,1)
5612 tmpa(is:ie,js:je,tile) = vort(is:ie,js:je)
5613 call wrtvar_ncdf(ncid, pv_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5616 call cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, 1, 1, gridstruct%grid_type, gridstruct%nested, flagstruct%c2l_ord, bd)
5619 ut(i,j,tile) = ua(i,j,1)
5620 vt(i,j,tile) = va(i,j,1)
5624 call wrtvar_ncdf(ncid, u_id, nout, is,ie, js,je, npx, npy, npz, nregions, ut(1:npx-1,1:npy-1,1:nregions), 3)
5625 call wrtvar_ncdf(ncid, v_id, nout, is,ie, js,je, npx, npy, npz, nregions, vt(1:npx-1,1:npy-1,1:nregions), 3)
5627 if ((
test_case >= 2) .and. (nt==0) )
then 5628 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/grav
5629 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5634 tmpa_3d(is:ie,js:je,1:npz,tile) = q(is:ie,js:je,1:npz,1)
5635 call wrtvar_ncdf(ncid, q_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5639 tmpa_3d(is:ie,js:je,1:npz,tile) = q(is:ie,js:je,1:npz,iq)
5640 call wrtvar_ncdf(ncid, tracers_ids(iq-1), nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5644 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/grav
5645 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, 1, nregions, tmpa, 3)
5648 tmpa(is:ie,js:je,tile) = ps(is:ie,js:je)
5649 call wrtvar_ncdf(ncid, ps_id, nout, is,ie, js,je, npx, npy, 1, nregions, tmpa, 3)
5651 tmpa_3d(is:ie,js:je,k,tile) = delp(is:ie,js:je,k)/grav
5653 call wrtvar_ncdf(ncid, delp_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5657 tmpa_3d(is:ie,js:je,k,tile) = pt(is:ie,js:je,k)
5659 call wrtvar_ncdf(ncid, pt_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5662 call cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, npz, gridstruct%grid_type, gridstruct%nested, flagstruct%c2l_ord)
5666 ut(i,j,k,tile) = ua(i,j,k)
5667 vt(i,j,k,tile) = va(i,j,k)
5671 call wrtvar_ncdf(ncid, u_id, nout, is,ie, js,je, npx, npy, npz, nregions, ut(1:npx-1,1:npy-1,1:npz,1:nregions), 4)
5672 call wrtvar_ncdf(ncid, v_id, nout, is,ie, js,je, npx, npy, npz, nregions, vt(1:npx-1,1:npy-1,1:npz,1:nregions), 4)
5679 tmpa_3d(i,j,k,tile) = rarea(i,j) * ( (v(i+1,j,k)*dy(i+1,j) - v(i,j,k)*dy(i,j)) - &
5680 (u(i,j+1,k)*dx(i,j+1) - u(i,j,k)*dx(i,j)) )
5684 call wrtvar_ncdf(ncid, pv_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5690 tmpa_3d(i,j,k,tile) = omga(i,j,k)
5694 call wrtvar_ncdf(ncid, om_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5700 #if defined(SW_DYNAMICS) 5706 deallocate( tmpa_3d )
5724 end subroutine output_ncdf
5735 subroutine output(dt, nt, maxnt, nout, u,v,pt,delp,q,phis,ps, uc,vc, ua,va, &
5736 npx, npy, npz, ng, ncnst, ndims, nregions, phis_lun, phi_lun, &
5737 pt_lun, pv_lun, uv_lun, gridstruct)
5739 real,
intent(IN) :: dt
5740 integer,
intent(IN) :: nt, maxnt
5741 integer,
intent(INOUT) :: nout
5743 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
5744 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
5745 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
5746 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
5747 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
5749 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
5750 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
5752 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
5753 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
5754 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
5755 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
5757 integer,
intent(IN) :: npx, npy, npz
5758 integer,
intent(IN) :: ng, ncnst
5759 integer,
intent(IN) :: ndims
5760 integer,
intent(IN) :: nregions
5761 integer,
intent(IN) :: phis_lun, phi_lun, pt_lun, pv_lun, uv_lun
5763 type(fv_grid_type),
target :: gridstruct
5765 real :: tmp(1-ng:npx +ng,1-ng:npy +ng,1:nregions)
5766 real :: tmpA(1-ng:npx-1+ng,1-ng:npy-1+ng,1:nregions)
5772 real :: ut(1:npx,1:npy,1:nregions)
5773 real :: vt(1:npx,1:npy,1:nregions)
5774 real :: utmp, vtmp, r, r0, dist, heading
5775 integer :: i,j,k,n,nreg
5776 real :: vort(isd:ied,jsd:jed)
5779 real :: x1,y1,z1,x2,y2,z2,ang
5781 real,
pointer,
dimension(:,:,:) :: agrid, grid
5782 real,
pointer,
dimension(:,:) :: area, rarea
5783 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5785 grid => gridstruct%grid
5786 agrid => gridstruct%agrid
5788 area => gridstruct%area
5792 dxa => gridstruct%dxa
5793 dya => gridstruct%dya
5794 rdxa => gridstruct%rdxa
5795 rdya => gridstruct%rdya
5796 dxc => gridstruct%dxc
5797 dyc => gridstruct%dyc
5799 cubed_sphere => gridstruct%cubed_sphere
5803 #if defined(SW_DYNAMICS) 5805 call atob_s(delp(:,:,1)/grav, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5806 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)/grav
5815 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
5816 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0) / grav
5830 dist = 2.0*pi*
radius* ((float(nt)/float(maxnt)))
5831 heading = 5.0*pi/2.0 -
alpha 5835 p2(1) = agrid(i,j,1)
5836 p2(2) = agrid(i,j,2)
5839 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
5841 phi0(i,j,1) = phis(i,j)
5853 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
5855 if (p /= 0.0) w_p = vtx/p
5856 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
5861 call atob_s(
phi0(:,:,1), tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5862 tmpa(is:ie,js:je,tile) =
phi0(is:ie,js:je,1)
5863 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5864 call atob_s(delp(:,:,1), tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5865 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)
5868 call wrt2d(phi_lun, nout, is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5874 vort(i,j) = f0(i,j) + (1./area(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
5875 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
5876 vort(i,j) = grav*vort(i,j)/delp(i,j,1)
5879 call atob_s(vort, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5880 call wrt2d(pv_lun, nout, is,ie+1, js,je+1, npx+1, npy+1, nregions, tmp(1:npx,1:npy,1:nregions))
5883 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy, ng)
5885 if (cubed_sphere)
then 5894 if (cubed_sphere)
call rotate_winds(utmp,vtmp, p1,p2,p3,p4, agrid(i,j,1:2), 2, 2)
5901 call wrt2d(uv_lun, 2*(nout-1) + 1, is,ie, js,je, npx, npy, nregions, ut(1:npx-1,1:npy-1,1:nregions))
5902 call wrt2d(uv_lun, 2*(nout-1) + 2, is,ie, js,je, npx, npy, nregions, vt(1:npx-1,1:npy-1,1:nregions))
5904 if ((
test_case >= 2) .and. (nt==0) )
then 5905 call atob_s(phis/grav, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5907 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/grav
5908 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5914 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/grav
5915 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5925 tmpa(is:ie,js:je,tile) = ps(is:ie,js:je)
5926 call wrt2d(phi_lun, (nout-1)*(npz+1) + 1, is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5928 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,k)/grav
5929 call wrt2d(phi_lun, (nout-1)*(npz+1) + 1 + k, is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5934 tmpa(is:ie,js:je,tile) = pt(is:ie,js:je,k)
5935 call wrt2d(pt_lun, (nout-1)*npz + (k-1) + 1, is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5940 call dtoa(u(isd,jsd,k), v(isd,jsd,k), ua(isd,jsd,k), va(isd,jsd,k), dx,dy,dxa,dya,dxc,dyc,npx, npy, ng)
5942 if (cubed_sphere)
then 5951 if (cubed_sphere)
call rotate_winds(utmp,vtmp, p1,p2,p3,p4, agrid(i,j,1:2), 2, 2)
5957 call wrt2d(uv_lun, 2*((nout-1)*npz + (k-1)) + 1, is,ie, js,je, npx, npy, nregions, ut(1:npx-1,1:npy-1,1:nregions))
5958 call wrt2d(uv_lun, 2*((nout-1)*npz + (k-1)) + 2, is,ie, js,je, npx, npy, nregions, vt(1:npx-1,1:npy-1,1:nregions))
5976 nullify(cubed_sphere)
5978 end subroutine output
5987 subroutine wrtvar_ncdf(ncid, varid, nrec, i1,i2, j1,j2, npx, npy, npz, ntiles, p, ndims)
5988 #include <netcdf.inc> 5989 integer,
intent(IN) :: ncid, varid
5990 integer,
intent(IN) :: nrec
5991 integer,
intent(IN) :: i1,i2,j1,j2
5992 integer,
intent(IN) :: npx
5993 integer,
intent(IN) :: npy
5994 integer,
intent(IN) :: npz
5995 integer,
intent(IN) :: ntiles
5996 real ,
intent(IN) :: p(npx-1,npy-1,npz,ntiles)
5997 integer,
intent(IN) :: ndims
6000 real(kind=4),
allocatable :: p_R4(:,:,:,:)
6002 integer :: istart(ndims+1), icount(ndims+1)
6004 allocate( p_r4(npx-1,npy-1,npz,ntiles) )
6007 p_r4(i1:i2,j1:j2,1:npz,tile) = p(i1:i2,j1:j2,1:npz,tile)
6008 call mp_gather(p_r4, i1,i2, j1,j2, npx-1, npy-1, npz, ntiles)
6011 istart(ndims+1) = nrec
6015 if (ndims == 3) icount(3) = ntiles
6016 if (ndims == 4) icount(4) = ntiles
6019 if (is_master())
then 6020 error = nf_put_vara_real(ncid, varid, istart, icount, p_r4)
6025 end subroutine wrtvar_ncdf
6034 subroutine wrt2d(iout, nrec, i1,i2, j1,j2, npx, npy, nregions, p)
6035 integer,
intent(IN) :: iout
6036 integer,
intent(IN) :: nrec
6037 integer,
intent(IN) :: i1,i2,j1,j2
6038 integer,
intent(IN) :: npx
6039 integer,
intent(IN) :: npy
6040 integer,
intent(IN) :: nregions
6041 real ,
intent(IN) :: p(npx-1,npy-1,nregions)
6043 real(kind=4) :: p_R4(npx-1,npy-1,nregions)
6049 p_r4(i,j,n) = p(i,j,n)
6054 call mp_gather(p_r4, i1,i2, j1,j2, npx-1, npy-1, nregions)
6056 if (is_master())
then 6057 write(iout,rec=nrec) p_r4(1:npx-1,1:npy-1,1:nregions)
6060 end subroutine wrt2d
6069 subroutine init_double_periodic(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
6070 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, &
6071 mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd)
6075 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6076 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6077 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
6078 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6079 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6080 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
6082 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
6084 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
6085 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
6086 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
6087 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
6088 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
6089 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6090 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6091 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6092 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6093 real ,
intent(inout) :: delz(bd%isd:,bd%jsd:,1:)
6094 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
6096 real ,
intent(inout) :: ak(npz+1)
6097 real ,
intent(inout) :: bk(npz+1)
6099 integer,
intent(IN) :: npx, npy, npz
6100 integer,
intent(IN) ::
ng, ncnst, nwat
6101 integer,
intent(IN) :: ndims
6102 integer,
intent(IN) :: nregions
6104 real,
intent(IN) :: dry_mass
6105 logical,
intent(IN) :: mountain
6106 logical,
intent(IN) :: moist_phys
6107 logical,
intent(IN) :: hydrostatic, hybrid_z
6108 integer,
intent(INOUT) :: ks
6109 integer,
intent(INOUT),
target :: tile_in
6110 real,
intent(INOUT) :: ptop
6112 type(domain2d),
intent(IN),
target :: domain_in
6117 real,
dimension(bd%is:bd%ie):: pm, qs
6118 real,
dimension(1:npz):: pk1, ts1, qs1
6120 real :: dist, r0, f0_const, prf, rgrav
6121 real :: ptmp, ze, zc, zm, utmp, vtmp
6122 real :: t00, p00, xmax, xc, xx, yy, pk0, pturb, ztop
6126 integer :: i, j, k, m, icenter, jcenter
6128 real,
pointer,
dimension(:,:,:) :: agrid, grid
6129 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
6130 real,
pointer,
dimension(:,:) :: rarea, fc, f0
6131 real,
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
6132 real,
pointer,
dimension(:,:,:,:) :: ew, es
6133 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
6135 logical,
pointer :: cubed_sphere, latlon
6137 type(domain2d),
pointer :: domain
6138 integer,
pointer :: tile
6140 logical,
pointer :: have_south_pole, have_north_pole
6142 integer,
pointer :: ntiles_g
6143 real,
pointer :: acapn, acaps, globalarea
6145 real(kind=R_GRID),
pointer :: dx_const, dy_const
6147 integer :: is, ie, js, je
6148 integer :: isd, ied, jsd, jed
6159 agrid => gridstruct%agrid
6160 grid => gridstruct%grid
6162 area => gridstruct%area_64
6166 dxa => gridstruct%dxa
6167 dya => gridstruct%dya
6168 rdxa => gridstruct%rdxa
6169 rdya => gridstruct%rdya
6170 dxc => gridstruct%dxc
6171 dyc => gridstruct%dyc
6177 dx_const => flagstruct%dx_const
6178 dy_const => flagstruct%dy_const
6183 have_south_pole => gridstruct%have_south_pole
6184 have_north_pole => gridstruct%have_north_pole
6186 ntiles_g => gridstruct%ntiles_g
6187 acapn => gridstruct%acapN
6188 acaps => gridstruct%acapS
6189 globalarea => gridstruct%globalarea
6191 f0_const = 2.*omega*sin(flagstruct%deglat/180.*pi)
6212 if (j>0 .and. j<5)
then 6214 if (i>0 .and. i<5)
then 6220 call mpp_update_domains( delp, domain )
6227 r0 = 5.*sqrt(dx_const**2 + dy_const**2)
6232 dist=(i-icenter)*dx_const*(i-icenter)*dx_const &
6233 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6234 dist=min(r0,sqrt(dist))
6235 phis(i,j)=1500.*(1. - (dist/r0))
6256 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
6257 delp, ak, bk, pt, delz, area,
ng, .false., hydrostatic, hybrid_z, domain)
6261 r0 = 100.*sqrt(dx_const**2 + dy_const**2)
6267 dist = (i-icenter)*dx_const*(i-icenter)*dx_const &
6268 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6269 dist = min(r0, sqrt(dist))
6271 prf = ak(k) + ps(i,j)*bk(k)
6272 if ( prf > 100.e2 )
then 6273 pt(i,j,k) = pt(i,j,k) + 0.01*(1. - (dist/r0)) * prf/ps(i,j)
6279 if ( hydrostatic )
then 6280 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6281 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6282 moist_phys, .true., nwat , domain)
6285 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6286 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6287 moist_phys, hydrostatic, nwat, domain, .true. )
6294 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6297 call qsmith((ie-is+1)*(je-js+1), npz, &
6298 ie-is+1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6300 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6303 q(i,j,k,1) = max(2.e-6, 0.8*pm(i)/ps(i,j)*qs(i) )
6319 if ( .not. hydrostatic ) w(:,:,:) = 0.
6331 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6340 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6346 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6347 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6348 moist_phys, .false., nwat, domain)
6351 r0 = 5.*max(dx_const, dy_const)
6360 zm = ze - 0.5*delz(i,j,k)
6361 ze = ze - delz(i,j,k)
6362 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + &
6365 if ( dist <= r0 )
then 6366 pt(i,j,k) = pt(i,j,k) + 5.*(1.-dist/r0)
6389 ze1(k) = ze1(k+1) + ztop/
real(npz)
6403 delz(i,j,k) = ze1(k+1) - ze1(k)
6404 pk(i,j,k) = pk(i,j,k+1) + grav*delz(i,j,k)/(cp_air*t00)*pk0
6405 pe(i,k,j) = pk(i,j,k)**(1./kappa)
6411 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
6416 peln(i,k,j) = log(pe(i,k,j))
6425 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6426 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6437 zm = (0.5*(ze1(k)+ze1(k+1))-3.e3) / 2.e3
6441 xx = (dx_const * (0.5+
real(i-1)) - xc) / 4.e3
6442 yy = (dy_const * (0.5+
real(j-1)) - xc) / 4.e3
6443 dist = sqrt( xx**2 + yy**2 + zm**2 )
6444 if ( dist<=1. )
then 6445 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
6448 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
6457 zvir = rvgas/rdgas - 1.
6463 pk(i,j,1) = ptop**kappa
6465 peln(i,1,j) = log(ptop)
6472 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6473 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6474 peln(i,k+1,j) = log(pe(i,k+1,j))
6475 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
6483 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6496 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6503 ze1(k) = ze1(k+1) - delz(is,js,k)
6507 zm = 0.5*(ze1(k)+ze1(k+1))
6508 utmp = us0*tanh(zm/3.e3)
6516 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6517 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6518 .true., hydrostatic, nwat, domain)
6524 icenter = (npx-1)/3 + 1
6525 jcenter = (npy-1)/2 + 1
6527 zm = 0.5*(ze1(k)+ze1(k+1))
6528 ptmp = ( (zm-zc)/zc ) **2
6529 if ( ptmp < 1. )
then 6532 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6533 if ( dist < 1. )
then 6534 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6546 zvir = rvgas/rdgas - 1.
6552 pk(i,j,1) = ptop**kappa
6554 peln(i,1,j) = log(ptop)
6561 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6562 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6563 peln(i,k+1,j) = log(pe(i,k+1,j))
6564 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
6572 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6584 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6591 ze1(k) = ze1(k+1) - delz(is,js,k)
6597 zm = 0.5*(ze1(k)+ze1(k+1))
6598 if ( zm .le. 2.e3 )
then 6599 utmp = 8.*(1.-cos(pi*zm/4.e3))
6600 vtmp = 8.*sin(pi*zm/4.e3)
6601 elseif (zm .le. 6.e3 )
then 6602 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
6611 u(i,j,k) = utmp - 8.
6617 v(i,j,k) = vtmp - 4.
6623 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6624 pe, peln, pk, pkz, kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6625 .true., hydrostatic, nwat, domain)
6631 icenter = (npx-1)/2 + 1
6632 jcenter = (npy-1)/2 + 1
6634 zm = 0.5*(ze1(k)+ze1(k+1))
6635 ptmp = ( (zm-zc)/zc ) **2
6636 if ( ptmp < 1. )
then 6639 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6640 if ( dist < 1. )
then 6641 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6663 if (.not.hybrid_z)
call mpp_error(fatal,
'hybrid_z must be .TRUE.')
6668 call mpp_error(fatal,
'npz must be == 101 ')
6681 peln(i,npz+1,j) = log(p00)
6688 peln(i,k,j) = peln(i,k+1,j) + grav*delz(i,j,k)/(rdgas*t00)
6689 pe(i,k,j) = exp(peln(i,k,j))
6690 pk(i,j,k) = pe(i,k,j)**kappa
6699 if ( is_master() )
write(*,*)
'LES testcase: computed model top (mb)=', ptop/100.
6704 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6705 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6713 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6716 call qsmith((ie-is+1)*(je-js+1), npz, &
6717 ie-is+1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6719 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6722 if ( pm(i) > 100.e2 )
then 6723 q(i,j,k,1) = 0.9*qs(i)
6740 zm = 0.5*(ze0(i,j,k)+ze0(i,j,k+1))
6741 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + (zm-zc)**2
6743 if ( dist <= r0 )
then 6744 pt(i,j,k) = pt(i,j,k) + 2.0*(1.-dist/r0)
6782 nullify(have_south_pole)
6783 nullify(have_north_pole)
6793 integer,
intent(in):: km
6794 real,
intent(in):: p00
6795 real,
intent(inout),
dimension(km+1):: pe
6796 real,
intent(in),
dimension(km+1):: ze
6799 real,
intent(out),
dimension(km):: pt, qz
6801 integer,
parameter:: nx = 5
6802 real,
parameter:: qst = 1.0e-6
6803 real,
parameter:: qv0 = 1.4e-2
6804 real,
parameter:: ztr = 12.e3
6805 real,
parameter:: ttr = 213.
6806 real,
parameter:: ptr = 343.
6807 real,
parameter:: pt0 = 300.
6808 real,
dimension(km):: zs, rh, temp, dp, dp0
6809 real,
dimension(km+1):: peln, pk
6810 real:: qs, zvir, fac_z, pk0, temp1, pm
6813 zvir = rvgas/rdgas - 1.
6815 if ( (is_master()) )
then 6816 write(*,*)
'Computing sounding for HIWPP super-cell test using p00=', p00
6823 zs(k) = 0.5*(ze(k)+ze(k+1))
6825 if ( zs(k) .gt. ztr )
then 6827 pt(k) = ptr*exp(grav*(zs(k)-ztr)/(cp_air*ttr))
6830 fac_z = (zs(k)/ztr)**1.25
6831 pt(k) = pt0 + (ptr-pt0)* fac_z
6832 rh(k) = 1. - 0.75 * fac_z
6834 qz(k) = qv0 - (qv0-qst)*fac_z
6836 if ( is_master() )
write(*,*) zs(k), pt(k), qz(k)
6841 #ifdef USE_MOIST_P00 6848 peln(km+1) = log(p00)
6853 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k)*(1.+zvir*qz(k)))
6854 peln(k) = log(pk(k)) / kappa
6855 pe(k) = exp(peln(k))
6858 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
6859 temp(k) = pt(k)*pm**kappa
6861 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6862 qz(k) = min( qv0, rh(k)*qs )
6863 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
6870 peln(km+1) = log(p00)
6874 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k))
6875 peln(k) = log(pk(k)) / kappa
6876 pe(k) = exp(peln(k))
6879 dp0(k) = pe(k+1) - pe(k)
6880 pm = dp0(k)/(peln(k+1)-peln(k))
6881 temp(k) = pt(k)*pm**kappa
6883 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6884 qz(k) = min( qv0, rh(k)*qs )
6890 dp(k) = dp0(k)*(1. + qz(k))
6891 pe(k+1) = pe(k) + dp(k)
6894 pk(km+1) = pe(km+1)**kappa
6895 peln(km+1) = log(pe(km+1))
6899 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k)*(1.+zvir*qz(k)))
6900 peln(k) = log(pk(k)) / kappa
6901 pe(k) = exp(peln(k))
6904 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
6905 temp(k) = pt(k)*pm**kappa
6907 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6908 qz(k) = min( qv0, rh(k)*qs )
6909 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
6914 if ( is_master() )
then 6915 write(*,*)
'Super_K: computed ptop (mb)=', 0.01*pe(1),
' PS=', 0.01*pe(km+1)
6916 call prt_m1(
'1D Sounding T0', temp, 1, km, 1, 1, 0, 1, 1.)
6921 subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
6922 delz, zvir, ptop, ak, bk, agrid)
6923 integer,
intent(in):: is, ie, js, je, ng, km
6924 real,
intent(in),
dimension(km ):: ts1, qs1, uz1, dudz
6925 real,
intent(in),
dimension(km+1):: ze1
6926 real,
intent(in):: zvir, ps0
6927 real,
intent(inout):: ptop
6928 real(kind=R_GRID),
intent(in):: agrid(is-ng:ie+ng,js-ng:je+ng,2)
6929 real,
intent(inout),
dimension(km+1):: ak, bk
6930 real,
intent(inout),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delz
6931 real,
intent(out),
dimension(is:ie,js:je,km+1):: pk
6933 real,
intent(inout),
dimension(is-1:ie+1,km+1,js-1:je+1):: pe
6935 integer,
parameter:: nt=5
6936 integer,
parameter:: nlat=1001
6937 real,
dimension(nlat,km):: pt2, pky, dzc
6938 real,
dimension(nlat,km+1):: pk2, pe2, peln2, pte
6939 real,
dimension(km+1):: pe1
6940 real:: lat(nlat), latc(nlat-1)
6941 real:: fac_y, dlat, dz0, pk0, tmp1, tmp2, tmp3, pint
6942 integer::i,j,k,n, jj, k1
6946 dz0 = ze1(km) - ze1(km+1)
6949 dlat = 0.5*pi/
real(nlat-1)
6951 lat(j) = dlat*
real(j-1)
6953 dzc(j,k) = ze1(k) - ze1(k+1)
6957 latc(j) = 0.5*(lat(j)+lat(j+1))
6966 if ( is_master() )
then 6968 call prt_m1(
'Super_K PT0', pt2, 1, nlat, 1, km, 0, 1, tmp1)
6975 call ppme(pt2, pte, dzc, nlat, km)
6978 tmp1 = 0.5*(pte(j-1,k ) + pte(j,k ))
6979 tmp3 = 0.5*(pte(j-1,k+1) + pte(j,k+1))
6980 pt2(j,k) = pt2(j-1,k) + dlat/(2.*grav)*sin(2.*latc(j-1))*uz1(k)* &
6981 ( uz1(k)*(tmp1-tmp3)/dzc(j,k) - (pt2(j-1,k)+pt2(j,k))*dudz(k) )
6984 if ( is_master() )
then 6985 call prt_m1(
'Super_K PT', pt2, 1, nlat, 1, km, 0, 1, pk0/cp_air)
6991 pk2(1,km+1) = ps0**kappa
6993 pk2(j,km+1) = pk2(j-1,km+1) - dlat*uz1(km)*uz1(km)*sin(2.*latc(j-1)) &
6994 / (pt2(j-1,km) + pt2(j,km))
6999 pk2(j,k) = pk2(j,k+1) - grav*dzc(j,k)/pt2(j,k)
7005 peln2(j,k) = log(pk2(j,k)) / kappa
7006 pe2(j,k) = exp(peln2(j,k))
7012 pky(j,k) = (pk2(j,k+1)-pk2(j,k))/(kappa*(peln2(j,k+1)-peln2(j,k)))
7013 pt2(j,k) = pt2(j,k)*pky(j,k)/(cp_air*(1.+zvir*qs1(k)))
7021 if ( is_master() )
then 7022 write(*,*)
'SuperK ptop at EQ=', 0.01*pe1(1),
'new ptop=', 0.01*ptop
7023 call prt_m1(
'Super_K pe', pe2, 1, nlat, 1, km+1, 0, 1, 0.01)
7024 call prt_m1(
'Super_K Temp', pt2, 1, nlat, 1, km, 0, 1, 1.)
7031 if (abs(agrid(i,j,2))>=lat(jj) .and. abs(agrid(i,j,2))<=lat(jj+1) )
then 7033 fac_y = (abs(agrid(i,j,2))-lat(jj)) / dlat
7035 pt(i, j,k) = pt2(jj, k) + fac_y*(pt2(jj+1, k)-pt2(jj,k))
7038 pe(i,k,j) = pe2(jj,k) + fac_y*(pe2(jj+1,k)-pe2(jj,k))
7061 bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint)
7062 ak(k) = pe1(k) - bk(k) * pe1(km+1)
7063 if ( is_master() )
write(*,*) k, ak(k), bk(k)
7076 subroutine superk_u(km, zz, um, dudz)
7077 integer,
intent(in):: km
7078 real,
intent(in):: zz(km)
7079 real,
intent(out):: um(km), dudz(km)
7081 real,
parameter:: zs = 5.e3
7082 real,
parameter:: us = 30.
7089 if ( zz(k) .gt. zs+1.e3 )
then 7092 elseif ( abs(zz(k)-zs) .le. 1.e3 )
then 7093 um(k) = us*(-4./5. + 3.*zz(k)/zs - 5./4.*(zz(k)/zs)**2)
7094 dudz(k) = us/zs*(3. - 5./2.*zz(k)/zs)
7103 um(k) = us*tanh( zz(k)/zs ) - uc
7104 dudz(k) = (us/zs)/cosh(zz(k)/zs)**2
7112 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7113 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7114 hydrostatic, nwat, adiabatic, do_pert, domain)
7116 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7117 real,
intent(IN) :: ptop
7118 real,
intent(IN),
dimension(npz+1) :: ak, bk
7119 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7120 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w, delz
7121 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7122 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7123 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7124 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7125 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7126 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7127 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7128 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7129 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7130 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7131 logical,
intent(IN) :: hydrostatic,adiabatic,do_pert
7132 type(domain2d),
intent(INOUT) :: domain
7134 real,
parameter :: p0 = 1.e5
7135 real,
parameter :: u0 = 35.
7136 real,
parameter :: b = 2.
7137 real,
parameter :: KK = 3.
7138 real,
parameter :: Te = 310.
7139 real,
parameter :: Tp = 240.
7140 real,
parameter :: T0 = 0.5*(te + tp)
7141 real,
parameter :: up = 1.
7142 real,
parameter :: zp = 1.5e4
7143 real(kind=R_GRID),
parameter :: lamp = pi/9.
7144 real(kind=R_GRID),
parameter :: phip = 2.*lamp
7145 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7146 real,
parameter :: Rp =
radius/10.
7147 real,
parameter :: lapse = 5.e-3
7148 real,
parameter :: dT = 4.8e5
7149 real,
parameter :: phiW = 2.*pi/9.
7150 real,
parameter :: pW = 34000.
7151 real,
parameter :: q0 = .018
7152 real,
parameter :: qt = 1.e-12
7153 real,
parameter :: ptrop = 1.e4
7155 real,
parameter :: zconv = 1.e-6
7156 real,
parameter :: rdgrav = rdgas/grav
7157 real,
parameter :: zvir = rvgas/rdgas - 1.
7158 real,
parameter :: rrdgrav = grav/rdgas
7160 integer :: i,j,k,iter, sphum, cl, cl2, n
7161 real :: p,z,z0,ziter,piter,titer,uu,vv,pl,pt_u,pt_v
7162 real(kind=R_GRID),
dimension(2) :: pa
7163 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7164 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2
7165 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7166 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2
7167 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7186 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7197 peln(i,1,j) = log(ptop)
7198 pk(i,j,1) = ptop**kappa
7202 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7205 pk(i,j,k) = exp(kappa*log(pe(i,k,j)))
7206 peln(i,k,j) = log(pe(i,k,j))
7214 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
7236 z = ziter + (piter - p)*rdgrav*titer/piter
7242 if (abs(z - ziter) < zconv)
exit 7253 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7264 peln_u(i,j) = log(p0)
7279 p = ak(k) + ps_u(i,j)*bk(k)
7288 z = ziter + (piter - p)*rdgrav*titer/piter
7289 if (abs(z - ziter) < zconv)
exit 7292 pt_u = rrdgrav * ( z - gz_u(i,j) ) / (peln_u(i,j) - pl)
7299 u(i,j,k) = u1(i,j)*uu
7312 peln_v(i,j) = log(p0)
7327 p = ak(k) + ps_v(i,j)*bk(k)
7336 z = ziter + (piter - p)*rdgrav*titer/piter
7337 if (abs(z - ziter) < zconv)
exit 7340 pt_v = rrdgrav * ( z - gz_v(i,j) ) / (peln_v(i,j) - pl)
7346 v(i,j,k) = v1(i,j)*uu
7364 if (.not. adiabatic)
then 7365 sphum = get_tracer_index(model_atmos,
'sphum')
7369 p = delp(i,j,k)/(peln(i,k+1,j) - peln(i,k,j))
7370 q(i,j,k,sphum) =
dcmip16_bc_sphum(p,ps(i,j),agrid(i,j,2),agrid(i,j,1))
7372 pt(i,j,k) = pt(i,j,k) / ( 1. + zvir*q(i,j,k,sphum))
7378 cl = get_tracer_index(model_atmos,
'cl')
7379 cl2 = get_tracer_index(model_atmos,
'cl2')
7380 if (cl > 0 .and. cl2 > 0)
then 7382 q, delp,nq,agrid(isd,jsd,1),agrid(isd,jsd,2))
7383 call mpp_update_domains(q,domain)
7387 if (.not. hydrostatic)
then 7392 delz(i,j,k) = gz(i,j,k) - gz(i,j,k+1)
7403 real,
intent(IN) :: z
7404 real(kind=R_GRID),
intent(IN) :: lat
7405 real :: it, t1, t2, tr, zsc
7407 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7408 zsc = z*grav/(b*rdgas*t0)
7409 tr = ( 1. - 2.*zsc**2.) * exp(-zsc**2. )
7411 t1 = (1./t0)*exp(lapse*z/t0) + (t0 - tp)/(t0*tp) * tr
7412 t2 = 0.5* ( kk + 2.) * (te - tp)/(te*tp) * tr
7420 real,
intent(IN) :: z
7421 real(kind=R_GRID),
intent(IN) :: lat
7422 real :: it, ti1, ti2, tir
7424 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7425 tir = z*exp(-(z*grav/(b*rdgas*t0))*(z*grav/(b*rdgas*t0)) )
7427 ti1 = 1./lapse* (exp(lapse*z/t0) - 1.) + tir*(t0-tp)/(t0*tp)
7428 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7436 real,
intent(IN) :: z, t
7437 real(kind=R_GRID),
intent(IN) :: lat
7438 real :: tir, ti2, uu, ur
7440 tir = z*exp(-(z*grav/(b*rdgas*t0))*(z*grav/(b*rdgas*t0)) )
7441 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7443 uu = grav*kk/
radius * ti2 * ( cos(lat)**(int(kk)-1) - cos(lat)**(int(kk)+1) ) * t
7444 ur = - omega *
radius * cos(lat) + sqrt( (omega*
radius*cos(lat))**2 +
radius*cos(lat)*uu)
7452 real,
intent(IN) :: z
7453 real(kind=R_GRID),
intent(IN) :: lat, lon
7455 real(kind=R_GRID) :: dst, pphere(2)
7458 zz = max(1. - 3.*zrat*zrat + 2.*zrat*zrat*zrat, 0.)
7460 pphere = (/ lon, lat /)
7469 real,
intent(IN) :: p, ps
7470 real(kind=R_GRID),
intent(IN) :: lat, lon
7485 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7486 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7487 hydrostatic, nwat, adiabatic)
7489 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7490 real,
intent(IN) :: ptop
7491 real,
intent(IN),
dimension(npz+1) :: ak, bk
7492 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7493 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w, delz
7494 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7495 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7496 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7497 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7498 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7499 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7500 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7501 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7502 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7503 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7504 logical,
intent(IN) :: hydrostatic,adiabatic
7506 real,
parameter :: zt = 15000
7507 real,
parameter :: q0 = 0.021
7508 real,
parameter :: qt = 1.e-11
7509 real,
parameter :: T0 = 302.15
7510 real,
parameter :: Tv0 = 302.15*(1.+0.608*q0)
7511 real,
parameter :: Ts = 302.15
7512 real,
parameter :: zq1 = 3000.
7513 real,
parameter :: zq2 = 8000.
7514 real,
parameter :: lapse = 7.e-3
7515 real,
parameter :: Tvt = tv0 - lapse*zt
7516 real,
parameter :: pb = 101500.
7517 real,
parameter :: ptt = pb*(tvt/tv0)**(grav/rdgas/lapse)
7518 real(kind=R_GRID),
parameter :: lamp = pi
7519 real(kind=R_GRID),
parameter :: phip = pi/18.
7520 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7521 real,
parameter :: dp = 1115.
7522 real,
parameter :: rp = 282000.
7523 real,
parameter :: zp = 7000.
7524 real,
parameter :: fc = 2.*omega*sin(phip)
7526 real,
parameter :: zconv = 1.e-6
7527 real,
parameter :: rdgrav = rdgas/grav
7528 real,
parameter :: rrdgrav = grav/rdgas
7530 integer :: i,j,k,iter, sphum, cl, cl2, n
7531 real :: p,z,z0,ziter,piter,titer,uu,vv,pl, r
7532 real(kind=R_GRID),
dimension(2) :: pa
7533 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7534 real,
dimension(is:ie,js:je) :: rc
7535 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2, rc_u
7536 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7537 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2, rc_v
7538 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7556 ps(i,j) = pb - dp*exp( -sqrt((rc(i,j)/rp)**3) )
7564 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7575 peln(i,1,j) = log(ptop)
7576 pk(i,j,1) = ptop**kappa
7580 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7583 pk(i,j,k) = exp(kappa*log(pe(i,k,j)))
7584 peln(i,k,j) = log(pe(i,k,j))
7592 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
7614 z = ziter + (piter - p)*rdgrav*titer/piter
7620 if (abs(z - ziter) < zconv)
exit 7631 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7649 p_u(i,j) = pb - dp*exp( -sqrt((rc_u(i,j)/rp)**3) )
7650 peln_u(i,j) = log(p_u(i,j))
7651 ps_u(i,j) = p_u(i,j)
7658 p = ak(k) + ps_u(i,j)*bk(k)
7667 z = ziter + (piter - p)*rdgrav*titer/piter
7668 if (abs(z - ziter) < zconv)
exit 7672 u(i,j,k) = u1(i,j)*uu + u2(i,j)*vv
7692 p_v(i,j) = pb - dp*exp( - sqrt((rc_v(i,j)/rp)**3) )
7693 peln_v(i,j) = log(p_v(i,j))
7694 ps_v(i,j) = p_v(i,j)
7701 p = ak(k) + ps_v(i,j)*bk(k)
7710 z = ziter + (piter - p)*rdgrav*titer/piter
7711 if (abs(z - ziter) < zconv)
exit 7715 v(i,j,k) = v1(i,j)*uu + v2(i,j)*vv
7733 if (.not. adiabatic)
then 7734 sphum = get_tracer_index(model_atmos,
'sphum')
7738 z = 0.5*(gz(i,j,k) + gz(i,j,k+1))
7746 if (.not. hydrostatic)
then 7751 delz(i,j,k) = gz(i,j,k) - gz(i,j,k+1)
7762 real,
intent(IN) :: z, r
7763 real :: tv, term1, term2
7771 term1 = grav*zp*zp* ( 1. - pb/dp * exp( sqrt(r/rp)**3 + (z/zp)**2 ) )
7772 term2 = 2*rdgas*tv*z
7780 real,
intent(IN) :: z, r
7783 dcmip16_tc_pressure = pb*exp(grav/(rdgas*lapse) * log( (tv0-lapse*z)/tv0) ) -dp* exp(-sqrt((r/rp)**3) - (z/zp)**2) * &
7784 exp( grav/(rdgas*lapse) * log( (tv0-lapse*z)/tv0) )
7793 real,
intent(IN) :: z, r
7794 real(kind=R_GRID),
intent(IN) :: lon, lat
7795 real,
intent(OUT) :: uu, vv
7796 real :: rfac, Tvrd, vt, fr5, d1, d2, d
7797 real(kind=R_GRID) :: dst, pphere(2)
7805 rfac = sqrt(r/rp)**3
7808 tvrd = (tv0 - lapse*z)*rdgas
7810 vt = -fr5 + sqrt( fr5**2 - (1.5 * rfac * tvrd) / &
7811 ( 1. + 2*tvrd*z/(grav*zp**2) - pb/dp*exp( rfac + (z/zp)**2) ) )
7813 d1 = sin(phip)*cos(lat) - cos(phip)*sin(lat)*cos(lon - lamp)
7814 d2 = cos(phip)*sin(lon - lamp)
7815 d = max(1.e-25,sqrt(d1*d1 + d2*d2))
7824 real,
intent(IN) :: z
7835 subroutine init_latlon(u,v,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
7836 gridstruct, npx, npy, npz, ng, ncnst, ndims, nregions, dry_mass, &
7837 mountain, moist_phys, hybrid_z, delz, ze0, domain_in, tile_in)
7839 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
7840 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
7841 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
7842 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
7843 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
7845 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
7847 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
7848 real ,
intent(INOUT) :: pe(is-1:ie+1,npz+1,js-1:je+1)
7849 real ,
intent(INOUT) :: pk(is:ie ,js:je ,npz+1)
7850 real ,
intent(INOUT) :: peln(is :ie ,npz+1 ,js:je)
7851 real ,
intent(INOUT) :: pkz(is:ie ,js:je ,npz )
7852 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
7853 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
7854 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
7855 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
7856 real ,
intent(inout) :: delz(isd:,jsd:,1:)
7857 real ,
intent(inout) :: ze0(is:,js:,1:)
7859 real ,
intent(IN) :: ak(npz+1)
7860 real ,
intent(IN) :: bk(npz+1)
7862 integer,
intent(IN) :: npx, npy, npz
7863 integer,
intent(IN) ::
ng, ncnst
7864 integer,
intent(IN) :: ndims
7865 integer,
intent(IN) :: nregions
7866 integer,
target,
intent(IN):: tile_in
7868 real,
intent(IN) :: dry_mass
7869 logical,
intent(IN) :: mountain
7870 logical,
intent(IN) :: moist_phys
7871 logical,
intent(IN) :: hybrid_z
7874 type(domain2d),
intent(IN),
target :: domain_in
7876 real,
pointer,
dimension(:,:,:) :: agrid, grid
7877 real,
pointer,
dimension(:,:) :: area, rarea, fc, f0
7878 real,
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
7879 real,
pointer,
dimension(:,:,:,:) :: ew, es
7880 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
7882 logical,
pointer :: cubed_sphere, latlon
7884 type(domain2d),
pointer :: domain
7885 integer,
pointer :: tile
7887 logical,
pointer :: have_south_pole, have_north_pole
7889 integer,
pointer :: ntiles_g
7890 real,
pointer :: acapn, acaps, globalarea
7892 real(kind=R_GRID) :: p1(2), p2(2)
7896 agrid => gridstruct%agrid
7897 grid => gridstruct%grid
7899 area => gridstruct%area
7903 dxa => gridstruct%dxa
7904 dya => gridstruct%dya
7905 rdxa => gridstruct%rdxa
7906 rdya => gridstruct%rdya
7907 dxc => gridstruct%dxc
7908 dyc => gridstruct%dyc
7913 ntiles_g => gridstruct%ntiles_g
7914 acapn => gridstruct%acapN
7915 acaps => gridstruct%acapS
7916 globalarea => gridstruct%globalarea
7921 have_south_pole => gridstruct%have_south_pole
7922 have_north_pole => gridstruct%have_north_pole
7926 fc(i,j) = 2.*omega*( -cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) &
7927 +sin(grid(i,j,2))*cos(
alpha) )
7932 f0(i,j) = 2.*omega*( -cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) &
7933 +sin(agrid(i,j,2))*cos(
alpha) )
7948 p2(1) = agrid(i,j,1)
7949 p2(2) = agrid(i,j,2)
7952 delp(i,j,1) = phis(i,j) + 0.5*(1.0+cos(pi*r/r0))
7954 delp(i,j,1) = phis(i,j)
8005 nullify(have_south_pole)
8006 nullify(have_north_pole)
8019 real,
intent(INOUT) :: UBar
8020 real,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
8021 real,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
8022 real,
intent(INOUT) :: uc(isd:ied+1,jsd:jed )
8023 real,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1)
8024 real,
intent(INOUT) :: ua(isd:ied ,jsd:jed )
8025 real,
intent(INOUT) :: va(isd:ied ,jsd:jed )
8026 integer,
intent(IN) :: defOnGrid
8027 type(fv_grid_type),
intent(IN),
target :: gridstruct
8029 real :: p1(2),p2(2),p3(2),p4(2), pt(2)
8030 real :: e1(3), e2(3), ex(3), ey(3)
8036 real :: psi_b(isd:ied+1,jsd:jed+1), psi(isd:ied,jsd:jed), psi1, psi2
8038 real,
dimension(:,:,:),
pointer :: grid, agrid
8039 real,
dimension(:,:),
pointer :: area, dx, dy, dxc, dyc
8041 grid => gridstruct%grid
8042 agrid=> gridstruct%agrid
8044 area => gridstruct%area
8047 dxc => gridstruct%dxc
8048 dyc => gridstruct%dyc
8054 psi(i,j) = (-1.0 * ubar *
radius *( sin(agrid(i,j,2)) *cos(
alpha) - &
8055 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
8060 psi_b(i,j) = (-1.0 * ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
8061 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
8065 if ( defongrid == 1 )
then 8069 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
8070 if (dist==0) vc(i,j) = 0.
8076 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
8077 if (dist==0) uc(i,j) = 0.
8085 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
8086 if (dist==0) v(i,j) = 0.
8092 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
8093 if (dist==0) u(i,j) = 0.
8100 subroutine d2a2c(im,jm,km, ifirst,ilast, jfirst,jlast, ng, nested, &
8101 u,v, ua,va, uc,vc, gridstruct, domain)
8104 integer,
intent(IN) :: im,jm,km
8105 integer,
intent(IN) :: ifirst,ilast
8106 integer,
intent(IN) :: jfirst,jlast
8107 integer,
intent(IN) :: ng
8108 logical,
intent(IN) :: nested
8109 type(fv_grid_type),
intent(IN),
target :: gridstruct
8110 type(domain2d),
intent(INOUT) :: domain
8125 real ,
intent(inout) :: u(isd:ied,jsd:jed+1)
8126 real ,
intent(inout) :: v(isd:ied+1,jsd:jed)
8127 real ,
intent(inout) :: ua(isd:ied,jsd:jed)
8128 real ,
intent(inout) :: va(isd:ied,jsd:jed)
8129 real ,
intent(inout) :: uc(isd:ied+1,jsd:jed)
8130 real ,
intent(inout) :: vc(isd:ied,jsd:jed+1)
8135 real :: sinlon(im,jm)
8136 real :: coslon(im,jm)
8137 real :: sinl5(im,jm)
8138 real :: cosl5(im,jm)
8140 real :: tmp1(jsd:jed+1)
8141 real :: tmp2(jsd:jed)
8142 real :: tmp3(jsd:jed)
8144 real mag,mag1,mag2, ang,ang1,ang2
8146 integer i, j, k, im2
8160 real,
pointer,
dimension(:,:,:) :: agrid, grid
8161 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
8162 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
8163 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
8164 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
8166 logical,
pointer :: cubed_sphere, latlon
8168 logical,
pointer :: have_south_pole, have_north_pole
8170 integer,
pointer :: ntiles_g
8171 real,
pointer :: acapN, acapS, globalarea
8173 grid => gridstruct%grid
8174 agrid=> gridstruct%agrid
8176 area => gridstruct%area
8177 rarea => gridstruct%rarea
8182 ee1 => gridstruct%ee1
8183 ee2 => gridstruct%ee2
8186 en1 => gridstruct%en1
8187 en2 => gridstruct%en2
8191 dxa => gridstruct%dxa
8192 dya => gridstruct%dya
8193 rdxa => gridstruct%rdxa
8194 rdya => gridstruct%rdya
8195 dxc => gridstruct%dxc
8196 dyc => gridstruct%dyc
8198 cubed_sphere => gridstruct%cubed_sphere
8199 latlon => gridstruct%latlon
8201 have_south_pole => gridstruct%have_south_pole
8202 have_north_pole => gridstruct%have_north_pole
8204 ntiles_g => gridstruct%ntiles_g
8205 acapn => gridstruct%acapN
8206 acaps => gridstruct%acapS
8207 globalarea => gridstruct%globalarea
8209 if (cubed_sphere)
then 8211 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,im,jm,ng)
8212 if (.not. nested)
call fill_corners(ua, va, im, jm, vector=.true., agrid=.true.)
8213 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,im,jm,ng, nested, domain, nocomm=.true.)
8214 if (.not. nested)
call fill_corners(uc, vc, im, jm, vector=.true., cgrid=.true.)
8226 js2gcp1 = jfirst-ng-1
8232 jn2gsp1 = jlast+ng-1
8234 if (have_south_pole)
then 8242 if (have_north_pole)
then 8252 if ( ng == 1 .AND. ng > 1 )
THEN 8255 js2gc1 = jfirst-ng+1
8256 if (have_south_pole) js2gc1 = 2
8261 if ((have_south_pole) .or. (have_north_pole))
then 8263 call vpol5(u(1:im,:), v(1:im,:), im, jm, &
8264 coslon, sinlon, cosl5, sinl5, ng, ng, jfirst, jlast )
8265 call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, v(:,:))
8268 call dtoa(u, v, ua, va, dx,dy,dxa,dya,dxc,dyc,im, jm, ng)
8269 if (.not. nested)
call fill_corners(ua, va, im, jm, vector=.true., agrid=.true.)
8271 if ( have_south_pole )
then 8276 us = us + (ua(i+im2,2)-ua(i,2))*sinlon(i,2) &
8277 + (va(i,2)-va(i+im2,2))*coslon(i,2)
8278 vs = vs + (ua(i+im2,2)-ua(i,2))*coslon(i,2) &
8279 + (va(i+im2,2)-va(i,2))*sinlon(i,2)
8285 ua(i,1) = -us*sinlon(i,1) - vs*coslon(i,1)
8286 va(i,1) = us*coslon(i,1) - vs*sinlon(i,1)
8287 ua(i+im2,1) = -ua(i,1)
8288 va(i+im2,1) = -va(i,1)
8291 ua(im+1,1) = ua(1 ,1)
8292 va(im+1,1) = va(1 ,1)
8295 if ( have_north_pole )
then 8301 un = un + (ua(i+im2,j)-ua(i,j))*sinlon(i,j) &
8302 + (va(i+im2,j)-va(i,j))*coslon(i,j)
8303 vn = vn + (ua(i,j)-ua(i+im2,j))*coslon(i,j) &
8304 + (va(i+im2,j)-va(i,j))*sinlon(i,j)
8310 ua(i,jm) = -un*sinlon(i,jm) + vn*coslon(i,jm)
8311 va(i,jm) = -un*coslon(i,jm) - vn*sinlon(i,jm)
8312 ua(i+im2,jm) = -ua(i,jm)
8313 va(i+im2,jm) = -va(i,jm)
8315 ua(0 ,jm) = ua(im,jm)
8316 ua(im+1,jm) = ua(1 ,jm)
8317 va(im+1,jm) = va(1 ,jm)
8320 if (latlon)
call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, ua(:,:))
8321 if (latlon)
call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, va(:,:))
8324 call atoc(ua, va, uc, vc, dx,dy,dxa,dya,im, jm, ng, nested, domain, nocomm=.true.)
8328 if (.not. nested)
call fill_corners(uc, vc, im, jm, vector=.true., cgrid=.true.)
8332 end subroutine d2a2c 8334 subroutine atob_s(qin, qout, npx, npy, dxa, dya, nested, cubed_sphere, altInterp)
8335 integer,
intent(IN) :: npx, npy
8336 real ,
intent(IN) :: qin(isd:ied ,jsd:jed )
8337 real ,
intent(OUT) :: qout(isd:ied+1,jsd:jed+1)
8338 integer,
OPTIONAL,
intent(IN) :: altInterp
8339 logical,
intent(IN) :: nested, cubed_sphere
8340 real,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8344 real :: tmp1j(jsd:jed+1)
8345 real :: tmp2j(jsd:jed+1)
8346 real :: tmp3j(jsd:jed+1)
8347 real :: tmp1i(isd:ied+1)
8348 real :: tmp2i(isd:ied+1)
8349 real :: tmp3i(isd:ied+1)
8350 real :: tmpq(isd:ied ,jsd:jed )
8351 real :: tmpq1(isd:ied+1,jsd:jed+1)
8352 real :: tmpq2(isd:ied+1,jsd:jed+1)
8354 if (
present(altinterp))
then 8356 tmpq(:,:) = qin(:,:)
8358 if (.not. nested)
call fill_corners(tmpq , npx, npy, fill=xdir, agrid=.true.)
8364 if (.not. nested)
call fill_corners(tmpq , npx, npy, fill=ydir, agrid=.true.)
8367 tmp1j(jsd:jed) = 0.0
8368 tmp2j(jsd:jed) = tmpq(i,jsd:jed)
8369 tmp3j(jsd:jed) = dya(i,jsd:jed)
8371 tmpq2(i,jsd:jed) = tmp1j(jsd:jed)
8376 tmp1j(:) = tmpq1(i,:)
8377 tmp2j(:) = tmpq1(i,:)
8380 tmpq1(i,:) = tmp1j(:)
8385 tmp1i(:) = tmpq2(:,j)
8386 tmp2i(:) = tmpq2(:,j)
8389 tmpq2(:,j) = tmp1i(:)
8395 qout(i,j) = 0.5 * (tmpq1(i,j) + tmpq2(i,j))
8400 if (cubed_sphere .and. .not. nested)
then 8403 if ( (is==i) .and. (js==j) )
then 8404 qout(i,j) = (1./3.) * (qin(i,j) + qin(i-1,j) + qin(i,j-1))
8409 if ( (ie+1==i) .and. (js==j) )
then 8410 qout(i,j) = (1./3.) * (qin(i-1,j) + qin(i-1,j-1) + qin(i,j))
8415 if ( (is==i) .and. (je+1==j) )
then 8416 qout(i,j) = (1./3.) * (qin(i,j-1) + qin(i-1,j-1) + qin(i,j))
8421 if ( (ie+1==i) .and. (je+1==j) )
then 8422 qout(i,j) = (1./3.) * (qin(i-1,j-1) + qin(i,j-1) + qin(i-1,j))
8430 qout(i,j) = 0.25 * (qin(i-1,j) + qin(i-1,j-1) + &
8431 qin(i ,j) + qin(i ,j-1))
8435 if (.not. nested)
then 8438 if ( (is==i) .and. (js==j) )
then 8439 qout(i,j) = (1./3.) * (qin(i,j) + qin(i-1,j) + qin(i,j-1))
8444 if ( (ie+1==i) .and. (js==j) )
then 8445 qout(i,j) = (1./3.) * (qin(i-1,j) + qin(i-1,j-1) + qin(i,j))
8450 if ( (is==i) .and. (je+1==j) )
then 8451 qout(i,j) = (1./3.) * (qin(i,j-1) + qin(i-1,j-1) + qin(i,j))
8456 if ( (ie+1==i) .and. (je+1==j) )
then 8457 qout(i,j) = (1./3.) * (qin(i-1,j-1) + qin(i,j-1) + qin(i-1,j))
8470 subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, nested, domain)
8471 integer,
intent(IN) :: npx, npy, ng
8472 real ,
intent(IN) :: uin(isd:ied ,jsd:jed )
8473 real ,
intent(IN) :: vin(isd:ied ,jsd:jed )
8474 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed+1)
8475 real ,
intent(OUT) :: vout(isd:ied+1,jsd:jed )
8476 logical,
intent(IN) :: nested
8477 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8478 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dxc
8479 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dyc
8480 type(domain2d),
intent(INOUT) :: domain
8483 real :: tmp1i(isd:ied+1)
8484 real :: tmp2i(isd:ied)
8485 real :: tmp3i(isd:ied)
8486 real :: tmp1j(jsd:jed+1)
8487 real :: tmp2j(jsd:jed)
8488 real :: tmp3j(jsd:jed)
8492 tmp2i(:) = vin(:,j)*dxa(:,j)
8495 vout(:,j) = tmp1i(:)/dxc(:,j)
8499 tmp2j(:) = uin(i,:)*dya(i,:)
8502 uout(i,:) = tmp1j(:)/dyc(i,:)
8505 if (.not. nested)
call fill_corners(uout, vout, npx, npy, vector=.true., dgrid=.true.)
8513 subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng)
8514 integer,
intent(IN) :: npx, npy, ng
8515 real ,
intent(IN) :: uin(isd:ied ,jsd:jed+1)
8516 real ,
intent(IN) :: vin(isd:ied+1,jsd:jed )
8517 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed )
8518 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed )
8519 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dx, dyc
8520 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dy, dxc
8521 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8525 real :: tmp1i(isd:ied+1)
8526 real :: tmp2i(isd:ied+1)
8527 real :: tmp3i(isd:ied+1)
8528 real :: tmp1j(jsd:jed+1)
8529 real :: tmp2j(jsd:jed+1)
8530 real :: tmp3j(jsd:jed+1)
8537 uout(i,j) = 0.5*(uin(i,j)*dx(i,j)+uin(i,j+1)*dx(i,j+1))/dxa(i,j)
8538 vout(i,j) = 0.5*(vin(i,j)*dy(i,j)+vin(i+1,j)*dy(i+1,j))/dya(i,j)
8544 tmp2j(:) = uin(i,:)*dyc(i,:)
8547 uout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dya(i,jsd:jed)
8551 tmp2i(:) = vin(:,j)*dxc(:,j)
8554 vout(isd:ied,j) = tmp1i(isd+1:ied+1)/dxa(isd:ied,j)
8565 subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, nested, domain, noComm)
8566 integer,
intent(IN) :: npx, npy, ng
8567 real ,
intent(IN) :: uin(isd:ied ,jsd:jed )
8568 real ,
intent(IN) :: vin(isd:ied ,jsd:jed )
8569 real ,
intent(OUT) :: uout(isd:ied+1,jsd:jed )
8570 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed+1)
8571 logical,
intent(IN) :: nested
8572 logical,
OPTIONAL,
intent(IN) :: noComm
8573 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dx
8574 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dy
8575 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8576 type(domain2d),
intent(INOUT) :: domain
8581 real :: tmp1i(isd:ied+1)
8582 real :: tmp2i(isd:ied)
8583 real :: tmp3i(isd:ied)
8584 real :: tmp1j(jsd:jed+1)
8585 real :: tmp2j(jsd:jed)
8586 real :: tmp3j(jsd:jed)
8588 #if !defined(ALT_INTERP) 8593 uout(i,j) = ( uin(i,j)*dxa(i,j) + uin(i-1,j)*dxa(i-1,j) ) &
8594 / ( dxa(i,j) + dxa(i-1,j) )
8599 vout(i,j) = ( vin(i,j)*dya(i,j) + vin(i,j-1)*dya(i,j-1) ) &
8600 / ( dya(i,j) + dya(i,j-1) )
8612 vout(i,:) = tmp1j(:)
8619 tmp2i(:) = uin(:,j)*dya(:,j)
8622 uout(:,j) = tmp1i(:)/dy(:,j)
8626 tmp2j(:) = vin(i,:)*dxa(i,:)
8629 vout(i,:) = tmp1j(:)/dx(i,:)
8632 if (cubed_sphere .and. .not. nested)
then 8633 csfac = cos(30.0*pi/180.0)
8635 if ( (is==1) .and. (js==1) )
then 8638 uout(i,j)=uout(i,j)*csfac
8639 uout(i,j-1)=uout(i,j-1)*csfac
8640 vout(i,j)=vout(i,j)*csfac
8641 vout(i-1,j)=vout(i-1,j)*csfac
8643 if ( (is==1) .and. (je==npy-1) )
then 8646 uout(i,j)=uout(i,j)*csfac
8647 uout(i,j+1)=uout(i,j+1)*csfac
8648 vout(i,j+1)=vout(i,j+1)*csfac
8649 vout(i-1,j+1)=vout(i-1,j+1)*csfac
8651 if ( (ie==npx-1) .and. (je==npy-1) )
then 8654 uout(i+1,j)=uout(i+1,j)*csfac
8655 uout(i+1,j+1)=uout(i+1,j+1)*csfac
8656 vout(i,j+1)=vout(i,j+1)*csfac
8657 vout(i+1,j+1)=vout(i+1,j+1)*csfac
8659 if ( (ie==npx-1) .and. (js==1) )
then 8662 uout(i+1,j)=uout(i+1,j)*csfac
8663 uout(i+1,j-1)=uout(i+1,j-1)*csfac
8664 vout(i,j)=vout(i,j)*csfac
8665 vout(i+1,j)=vout(i+1,j)*csfac
8671 if (
present(nocomm))
then 8672 if (.not. nocomm)
call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
8674 call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
8676 if (.not. nested)
call fill_corners(uout, vout, npx, npy, vector=.true., cgrid=.true.)
8685 subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng)
8686 integer,
intent(IN) :: npx, npy, ng
8687 real ,
intent(IN) :: uin(isd:ied+1,jsd:jed )
8688 real ,
intent(IN) :: vin(isd:ied ,jsd:jed+1)
8689 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed )
8690 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed )
8691 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dxc, dy
8692 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dyc, dx
8693 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8697 real :: tmp1i(isd:ied+1)
8698 real :: tmp2i(isd:ied+1)
8699 real :: tmp3i(isd:ied+1)
8700 real :: tmp1j(jsd:jed+1)
8701 real :: tmp2j(jsd:jed+1)
8702 real :: tmp3j(jsd:jed+1)
8716 tmp2j(:) = vin(i,:)*dx(i,:)
8719 vout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dxa(i,jsd:jed)
8723 tmp2i(:) = uin(:,j)*dy(:,j)
8726 uout(isd:ied,j) = tmp1i(isd+1:ied+1)/dya(isd:ied,j)
8736 subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
8737 integer,
intent(IN) :: ndims
8738 real ,
intent(INOUT) :: myU
8739 real ,
intent(INOUT) :: myV
8740 real(kind=R_GRID) ,
intent(IN) :: p1(ndims)
8741 real(kind=R_GRID) ,
intent(IN) :: p2(ndims)
8742 real(kind=R_GRID) ,
intent(IN) :: p3(ndims)
8743 real(kind=R_GRID) ,
intent(IN) :: p4(ndims)
8744 real(kind=R_GRID) ,
intent(IN) :: t1(ndims)
8745 integer,
intent(IN) :: dir
8747 real(kind=R_GRID) :: ee1(3), ee2(3), ee3(3), elon(3), elat(3)
8749 real :: g11, g12, g21, g22
8755 elon(1) = -sin(t1(1) - pi)
8756 elon(2) = cos(t1(1) - pi)
8758 elat(1) = -sin(t1(2))*cos(t1(1) - pi)
8759 elat(2) = -sin(t1(2))*sin(t1(1) - pi)
8760 elat(3) = cos(t1(2))
8768 newu = myu*g11 + myv*g12
8769 newv = myu*g21 + myv*g22
8771 newu = ( myu*g22 - myv*g12)/(g11*g22 - g21*g12)
8772 newv = (-myu*g21 + myv*g11)/(g11*g22 - g21*g12)
8780 use mpp_parameter_mod,
only: dgrid_ne
8781 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
8782 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
8783 integer,
intent(IN) :: npx, npy
8784 type(domain2d),
intent(INOUT) :: domain
8786 call mpp_update_domains( u, v, domain, gridtype=dgrid_ne, complete=.true.)
8798 use mpp_parameter_mod,
only: dgrid_ne
8799 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
8800 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
8801 integer,
intent(IN) :: npx, npy, npz
8802 type(domain2d),
intent(INOUT) :: domain
8805 call mpp_update_domains( u, v, domain, gridtype=dgrid_ne, complete=.true.)
8814 real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
result (gsum)
8815 integer,
intent(IN) :: npx, npy
8816 integer,
intent(IN) :: ifirst, ilast
8817 integer,
intent(IN) :: jfirst, jlast
8818 integer,
intent(IN) :: isd, ied
8819 integer,
intent(IN) :: jsd, jed, tile
8820 real ,
intent(IN) :: p(ifirst:ilast,jfirst:jlast)
8826 real,
allocatable :: p_r8(:,:,:)
8828 real,
pointer,
dimension(:,:,:) :: agrid, grid
8829 real,
pointer,
dimension(:,:) :: area, rarea, fc, f0
8830 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
8832 logical,
pointer :: cubed_sphere, latlon
8834 logical,
pointer :: have_south_pole, have_north_pole
8836 integer,
pointer :: ntiles_g
8837 real,
pointer :: acapn, acaps, globalarea
8839 grid => gridstruct%grid
8840 agrid=> gridstruct%agrid
8842 area => gridstruct%area
8843 rarea => gridstruct%rarea
8850 dxa => gridstruct%dxa
8851 dya => gridstruct%dya
8852 rdxa => gridstruct%rdxa
8853 rdya => gridstruct%rdya
8854 dxc => gridstruct%dxc
8855 dyc => gridstruct%dyc
8857 cubed_sphere => gridstruct%cubed_sphere
8858 latlon => gridstruct%latlon
8860 have_south_pole => gridstruct%have_south_pole
8861 have_north_pole => gridstruct%have_north_pole
8863 ntiles_g => gridstruct%ntiles_g
8864 acapn => gridstruct%acapN
8865 acaps => gridstruct%acapS
8866 globalarea => gridstruct%globalarea
8868 allocate(p_r8(npx-1,npy-1,ntiles_g))
8875 gsum = gsum + p(1,1)*acaps
8876 gsum = gsum + p(1,npy-1)*acapn
8879 gsum = gsum + p(i,j)*cos(agrid(i,j,2))
8887 p_r8(i,j,n) = p(i,j)*area(i,j)
8891 call mp_gather(p_r8, ifirst,ilast, jfirst,jlast, npx-1, npy-1, ntiles_g)
8892 if (is_master())
then 8896 gsum = gsum + p_r8(i,j,n)
8900 gsum = gsum/globalarea
8902 call mpp_broadcast(gsum, mpp_root_pe())
8912 real(kind=R_GRID),
intent(in):: p1(2), p2(2), p3(2)
8913 real(kind=R_GRID),
intent(out):: uvect(3)
8916 real(kind=R_GRID) :: xyz1(3), xyz2(3), xyz3(3)
8923 uvect(n) = xyz3(n)-xyz1(n)
8936 integer,
intent(in):: np
8937 real(kind=R_GRID),
intent(inout):: e(3,np)
8943 pdot = sqrt(e(1,n)**2+e(2,n)**2+e(3,n)**2)
8945 e(k,n) = e(k,n) / pdot
8955 subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
8956 kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
8959 integer,
intent(in):: im, jm, km, nq
8960 integer,
intent(in):: ifirst, ilast
8961 integer,
intent(in):: jfirst, jlast
8962 integer,
intent(in):: kfirst, klast
8963 integer,
intent(in):: ng_e
8964 integer,
intent(in):: ng_w
8965 integer,
intent(in):: ng_s
8966 integer,
intent(in):: ng_n
8967 real,
intent(inout):: q_ghst(ifirst-ng_w:ilast+ng_e,jfirst-ng_s:jlast+ng_n,kfirst:klast,nq)
8968 real,
optional,
intent(in):: q(ifirst:ilast,jfirst:jlast,kfirst:klast,nq)
8982 if (
present(q))
then 8983 q_ghst(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq) = &
8984 q(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq)
8990 do j=jfirst-ng_s,jlast+ng_n
8992 q_ghst(ifirst-i,j,k,n) = q_ghst(ilast-i+1,j,k,n)
8995 q_ghst(ilast+i,j,k,n) = q_ghst(ifirst+i-1,j,k,n)
9012 integer,
intent(in):: ifirst,ilast
9013 real,
intent(out) :: qout(ifirst:)
9014 real,
intent(in) :: qin(ifirst:)
9015 real,
intent(in) :: dx(ifirst:)
9016 integer,
intent(in):: order
9019 real :: dm(ifirst:ilast),qmax,qmin
9020 real :: r3, da1, da2, a6da, a6, al, ar
9021 real :: qLa, qLb1, qLb2
9030 qout(i) = 0.5 * (qin(i-1) + qin(i))
9032 elseif (order==2)
then 9035 qout(i) = (dx(i-1)*qin(i-1) + dx(i)*qin(i))/(dx(i-1)+dx(i))
9037 elseif (order==3)
then 9040 do i=ifirst+1,ilast-1
9041 dm(i) = 0.25*(qin(i+1) - qin(i-1))
9046 do i=ifirst+1,ilast-1
9047 qmax = max(qin(i-1),qin(i),qin(i+1)) - qin(i)
9048 qmin = qin(i) - min(qin(i-1),qin(i),qin(i+1))
9049 dm(i) = sign(min(abs(dm(i)),qmin,qmax),dm(i))
9052 do i=ifirst+1,ilast-1
9053 qout(i) = 0.5*(qin(i-1)+qin(i)) + r3*(dm(i-1) - dm(i))
9060 qout(ifirst+1) = 0.5 * (qin(ifirst) + qin(ifirst+1))
9061 qout(ilast) = 0.5 * (qin(ilast-1) + qin(ilast))
9063 elseif (order==4)
then 9066 do i=ifirst+1,ilast-1
9067 dm(i) = ( (2.*dx(i-1) + dx(i) ) / &
9068 ( dx(i+1) + dx(i) ) ) * ( qin(i+1) - qin(i) ) + &
9069 ( (dx(i) + 2.*dx(i+1)) / &
9070 (dx(i-1) + dx(i) ) ) * ( qin(i) - qin(i-1) )
9071 dm(i) = ( dx(i) / ( dx(i-1) + dx(i) + dx(i+1) ) ) * dm(i)
9072 if ( (qin(i+1)-qin(i))*(qin(i)-qin(i-1)) > 0.)
then 9073 dm(i) = sign( min( abs(dm(i)), 2.*abs(qin(i)-qin(i-1)), 2.*abs(qin(i+1)-qin(i)) ) , dm(i) )
9079 do i=ifirst+2,ilast-1
9080 qla = ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) - &
9081 ( (dx(i+1) + dx(i)) / (2.*dx(i) + dx(i-1)) )
9082 qla = ( (2.*dx(i) * dx(i-1)) / (dx(i-1) + dx(i)) ) * qla * &
9084 qlb1 = dx(i-1) * ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) * &
9086 qlb2 = dx(i) * ( (dx(i) + dx(i+1)) / (dx(i-1) + 2.*dx(i)) ) * &
9089 qout(i) = 1. / ( dx(i-2) + dx(i-1) + dx(i) + dx(i+1) )
9090 qout(i) = qout(i) * ( qla - qlb1 + qlb2 )
9091 qout(i) = qin(i-1) + ( dx(i-1) / ( dx(i-1) + dx(i) ) ) * (qin(i) - qin(i-1)) + qout(i)
9094 elseif (order==5)
then 9097 do i=ifirst+1,ilast-1
9098 x = float(i-(ifirst+1))*float(ilast-ifirst+1-1)/float(ilast-ifirst-1)
9099 qout(i) = qin(ifirst+nint(x)) + (x - nint(x)) * (qin(ifirst+nint(x+1)) - qin(ifirst+nint(x)))
9118 subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, &
9119 ng_d, ng_s, jfirst, jlast)
9125 integer,
intent(in):: ng_s, ng_d
9126 real,
intent(in):: coslon(im,jm), sinlon(im,jm)
9127 real,
intent(in):: cosl5(im,jm),sinl5(im,jm)
9128 real,
intent(in):: u(im,jfirst-ng_d:jlast+ng_s)
9131 real,
intent(inout):: v(im,jfirst-ng_d:jlast+ng_d)
9136 real uanp(im), uasp(im), vanp(im), vasp(im)
9137 real un, vn, us, vs, r2im
9140 r2im = 0.5d0/dble(im)
9145 if ( jfirst-ng_d <= 1 )
then 9147 uasp(i) = u(i, 2) + u(i,3)
9151 vasp(i) = v(i, 2) + v(i+1,2)
9153 vasp(im) = v(im,2) + v(1,2)
9159 us = us + (uasp(i+imh)-uasp(i))*sinlon(i,1) &
9160 + (vasp(i)-vasp(i+imh))*coslon(i,1)
9161 vs = vs + (uasp(i+imh)-uasp(i))*coslon(i,1) &
9162 + (vasp(i+imh)-vasp(i))*sinlon(i,1)
9170 v(i, 1) = us*cosl5(i,1) - vs*sinl5(i,1)
9171 v(i+imh,1) = -v(i,1)
9176 if ( jlast+ng_d >= jm )
then 9179 uanp(i) = u(i,jm-1) + u(i,jm)
9183 vanp(i) = v(i,jm-1) + v(i+1,jm-1)
9185 vanp(im) = v(im,jm-1) + v(1,jm-1)
9192 un = un + (uanp(i+imh)-uanp(i))*sinlon(i,jm) &
9193 + (vanp(i+imh)-vanp(i))*coslon(i,jm)
9194 vn = vn + (uanp(i)-uanp(i+imh))*coslon(i,jm) &
9195 + (vanp(i+imh)-vanp(i))*sinlon(i,jm)
9203 v(i, jm) = -un*cosl5(i,jm) - vn*sinl5(i,jm)
9204 v(i+imh,jm) = -v(i,jm)
9209 end subroutine vpol5 9211 subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
9213 character(len=*),
intent(in):: qname
9214 integer,
intent(in):: is, ie, js, je
9215 integer,
intent(in):: n_g, km
9216 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
9217 real,
intent(in):: fac
9228 if( q(i,j,k) < qmin )
then 9230 elseif( q(i,j,k) > qmax )
then 9237 write(*,*) qname,
' max = ', qmax*fac,
' min = ', qmin*fac
9241 subroutine var_dz(km, ztop, ze)
9242 integer,
intent(in):: km
9243 real,
intent(in):: ztop
9244 real,
intent(out),
dimension(km+1):: ze
9246 real,
dimension(km):: dz, s_fac
9257 s_fac(k) = 1.05 * s_fac(k+1)
9259 s_fac(4) = 1.1*s_fac(5)
9260 s_fac(3) = 1.2*s_fac(4)
9261 s_fac(2) = 1.3*s_fac(3)
9262 s_fac(1) = 1.5*s_fac(2)
9266 sum1 = sum1 + s_fac(k)
9272 dz(k) = s_fac(k) * dz0
9277 ze(k) = ze(k+1) + dz(k)
9282 dz(k) = dz(k) * (ztop/ze(1))
9286 ze(k) = ze(k+1) + dz(k)
9290 call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
9292 if ( is_master() )
then 9293 write(*,*)
'var_dz: model top (km)=', ztop*0.001
9295 dz(k) = ze(k) - ze(k+1)
9296 write(*,*) k, 0.5*(ze(k)+ze(k+1)),
'dz=', dz(k)
9302 subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
9303 integer,
intent(in):: is, ie, js, je, km
9304 integer,
intent(in):: ntimes, i, j
9305 real,
intent(inout):: ze(is:ie,js:je,km+1)
9307 real,
parameter:: df = 0.25
9310 integer k, n, k1, k2
9314 dz(k) = ze(i,j,k+1) - ze(i,j,k)
9323 flux(k) = df*(dz(k) - dz(k-1))
9327 dz(k) = dz(k) - flux(k) + flux(k+1)
9332 ze(i,j,k) = ze(i,j,k+1) - dz(k)
real, dimension(:,:), allocatable case9_b
subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng)
subroutine, public cart_to_latlon(np, q, xs, ys)
subroutine, public mid_pt_sphere(p1, p2, pm)
subroutine dcmip16_tc(delp, pt, u, v, q, w, delz, is, ie, js, je, isd, ied, jsd, jed, npz, nq, ak, bk, ptop, pk, peln, pe, pkz, gz, phis, ps, grid, agrid, hydrostatic, nwat, adiabatic)
logical, public bubble_do
subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
subroutine get_case9_b(B, agrid)
subroutine interp_left_edge_1d(qout, qin, dx, ifirst, ilast, order)
subroutine, public init_double_periodic(u, v, w, pt, delp, q, phis, ps, pe, peln, pk, pkz, uc, vc, ua, va, ak, bk, gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd)
subroutine, public init_case(u, v, w, pt, delp, q, phis, ps, pe, peln, pk, pkz, uc, vc, ua, va, ak, bk, gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, adiabatic, ks, npx_global, ptop, domain_in, tile_in, bd)
The module 'fv_mp_mod' is a single program multiple data (SPMD) parallel decompostion/communication m...
integer, parameter interporder
subroutine get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
subroutine get_vector_stats(varU, varUT, varV, varVT, npx, npy, ndims, nregions, vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
The type 'fv_grid_type' is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
real(kind=r_grid), parameter radius
subroutine, public hybrid_z_dz(km, dz, ztop, s_rate)
real tvort_orig
enstrophy (integral of total vorticity)
subroutine, public get_stats(dt, dtout, nt, maxnt, ndays, u, v, pt, delp, q, phis, ps, uc, vc, ua, va, npx, npy, npz, ncnst, ndims, nregions, gridstruct, stats_lun, consv_lun, monitorFreq, tile, domain, nested)
subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
real, dimension(:), allocatable lats_table
subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, nested, domain, noComm)
real, parameter pi_shift
3.0*pi/4.
subroutine mp_update_dwinds_3d(u, v, npx, npy, npz, domain)
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 set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
subroutine superk_u(km, zz, um, dudz)
real function dcmip16_bc_uwind(z, T, lat)
integer, public test_case
subroutine normalize_vect(np, e)
subroutine dcmip16_tc_uwind_pert(z, r, lon, lat, uu, vv)
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'.
subroutine var_dz(km, ztop, ze)
real function, public inner_prod(v1, v2)
subroutine get_scalar_stats(var, varT, npx, npy, ndims, nregions, vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine, public get_unit_vect2(e1, e2, uc)
integer, parameter, public r_grid
subroutine, public gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
The module 'fv_sg' performs FV sub-grid mixing.
integer, parameter initwindscase2
integer, public tracer_test
subroutine terminator_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, km, q, delp, ncnst, lon, lat)
real function gh_jet(npy, lat_in)
subroutine atob_s(qin, qout, npx, npy, dxa, dya, nested, cubed_sphere, altInterp)
subroutine, public make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
subroutine, public hydro_eq(km, is, ie, js, je, ps, hs, drym, delp, ak, bk, pt, delz, area, ng, mountain, hydrostatic, hybrid_z, domain)
The subroutine 'hydro_eq' computes a hydrostatically balanced and isothermal basic state from input h...
integer, parameter initwindscase6
subroutine rankine_vortex(ubar, r0, p1, u, v, grid)
subroutine d2a2c(im, jm, km, ifirst, ilast, jfirst, jlast, ng, nested, u, v, ua, va, uc, vc, gridstruct, domain)
real, parameter, public ptop_min
integer, parameter, public f_p
subroutine, public ppme(p, qe, delp, im, km)
subroutine, public compute_dz_l32(km, ztop, dz)
subroutine get_pt_on_great_circle(p1, p2, dist, heading, p3)
real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
real function dcmip16_tc_temperature(z, r)
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
integer, public wind_field
integer, parameter initwindscase9
real tmass_orig
total mass
The module 'fv_arrays' contains the 'fv_atmos_type' and associated datatypes.
real function, public great_circle_dist(q1, q2, radius)
real, public soliton_umax
integer, parameter initwindscase0
subroutine, public case9_forcing2(phis)
The module 'fv_eta' contains routine to set up the reference (Eulerian) pressure coordinate.
real function dcmip16_bc_sphum(p, ps, lat, lon)
real, dimension(:), allocatable, public pz0
subroutine pmxn(p, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
subroutine init_winds(UBar, u, v, ua, va, uc, vc, defOnGrid, npx, npy, ng, ndims, nregions, nested, gridstruct, domain, tile)
real, public soliton_size
integer, public nsolitons
subroutine, public case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain)
real function dcmip16_tc_sphum(z)
subroutine, public project_sphere_v(np, f, e)
real, dimension(:,:,:), allocatable ua0
Validating U-Wind.
integer, parameter, public ng
subroutine, public surfdrv(npx, npy, grid, agrid, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, phis, stretch_fac, nested, npx_global, domain, grid_number, bd, regional)
subroutine, public case9_forcing1(phis, time_since_start)
@ The module 'fv_diagnostics' contains routines to compute diagnosic fields.
subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, delz, zvir, ptop, ak, bk, agrid)
The module 'fv_grid_utils' contains routines for setting up and computing grid-related quantities...
subroutine get_unit_vector(p1, p2, p3, uvect)
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, nested, domain)
real, dimension(:), allocatable, public zz0
subroutine, public init_latlon(u, v, pt, delp, q, phis, ps, pe, peln, pk, pkz, uc, vc, ua, va, ak, bk, gridstruct, npx, npy, npz, ng, ncnst, ndims, nregions, dry_mass, mountain, moist_phys, hybrid_z, delz, ze0, domain_in, tile_in)
subroutine, public checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, nq, km, q, lon, lat, nx, ny, rn)
real function dcmip16_bc_uwind_pert(z, lat, lon)
subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
subroutine mp_update_dwinds_2d(u, v, npx, npy, domain)
real, dimension(:,:,:), allocatable va0
Validating V-Windfms_io_exit, get_tile_string, &.
real function dcmip16_bc_temperature(z, lat)
real, dimension(:,:,:), allocatable phi0
Validating Field.
subroutine init_latlon_winds(UBar, u, v, ua, va, uc, vc, defOnGrid, gridstruct)
real function dcmip16_tc_pressure(z, r)
real(kind=r_grid), parameter one
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
real function dcmip16_bc_pressure(z, lat)
subroutine dcmip16_bc(delp, pt, u, v, q, w, delz, is, ie, js, je, isd, ied, jsd, jed, npz, nq, ak, bk, ptop, pk, peln, pe, pkz, gz, phis, ps, grid, agrid, hydrostatic, nwat, adiabatic, do_pert, domain)
subroutine superk_sounding(km, pe, p00, ze, pt, qz)
subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, ng_d, ng_s, jfirst, jlast)
integer, parameter initwindscase1
real, dimension(:), allocatable gh_table
subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng)
subroutine, public p_var(km, ifirst, ilast, jfirst, jlast, ptop, ptop_min, delp, delz, pt, ps, pe, peln, pk, pkz, cappa, q, ng, nq, area, dry_mass, adjust_dry_mass, mountain, moist_phys, hydrostatic, nwat, domain, make_nh)
the subroutine 'p_var' computes auxiliary pressure variables for a hydrostatic state.
integer, parameter initwindscase5
subroutine, public compute_dz_l101(km, ztop, dz)
subroutine, public latlon2xyz(p, e, id)
The subroutine 'latlon2xyz' maps (lon, lat) to (x,y,z)
subroutine, public check_courant_numbers(uc, vc, ndt, n_split, gridstruct, npx, npy, npz, tile, noPrint)