26 is,js,ie,je, isd,jsd,ied,jed, &
27 domain_decomp, fill_corners, xdir, ydir, &
28 mp_stop, mp_reduce_sum, mp_reduce_max, mp_gather, mp_bcst
94 real(kind=R_GRID),
parameter ::
radius = cnst_radius
95 real(kind=R_GRID),
parameter ::
one = 1.d0
119 real,
allocatable,
dimension(:) ::
pz0,
zz0 134 real ,
allocatable ::
phi0(:,:,:)
135 real ,
allocatable ::
ua0(:,:,:)
136 real ,
allocatable ::
va0(:,:,:)
152 public :: output, output_ncdf
170 subroutine init_winds(UBar, u,v,ua,va,uc,vc, defOnGrid, npx, npy, ng, ndims, nregions, nested, gridstruct, domain, tile)
173 real ,
intent(INOUT) :: UBar
174 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
175 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
176 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed )
177 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1)
178 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed )
179 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed )
180 integer,
intent(IN) :: defOnGrid
181 integer,
intent(IN) :: npx, npy
182 integer,
intent(IN) :: ng
183 integer,
intent(IN) :: ndims
184 integer,
intent(IN) :: nregions
185 logical,
intent(IN) :: nested
186 type(fv_grid_type),
intent(IN),
target :: gridstruct
187 type(domain2d),
intent(INOUT) :: domain
188 integer,
intent(IN) :: tile
190 real(kind=R_GRID) :: p1(2), p2(2), p3(2), p4(2), pt(2)
191 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3)
197 real :: psi_b(isd:ied+1,jsd:jed+1), psi(isd:ied,jsd:jed), psi1, psi2
198 integer :: is2, ie2, js2, je2
200 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
201 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
202 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
203 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
204 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
206 logical,
pointer :: cubed_sphere, latlon
208 logical,
pointer :: have_south_pole, have_north_pole
210 integer,
pointer :: ntiles_g
211 real,
pointer :: acapN, acapS, globalarea
213 grid => gridstruct%grid_64
214 agrid=> gridstruct%agrid_64
216 area => gridstruct%area
217 rarea => gridstruct%rarea
222 ee1 => gridstruct%ee1
223 ee2 => gridstruct%ee2
226 en1 => gridstruct%en1
227 en2 => gridstruct%en2
231 dxa => gridstruct%dxa
232 dya => gridstruct%dya
233 rdxa => gridstruct%rdxa
234 rdya => gridstruct%rdya
235 dxc => gridstruct%dxc
236 dyc => gridstruct%dyc
238 cubed_sphere => gridstruct%cubed_sphere
239 latlon => gridstruct%latlon
241 have_south_pole => gridstruct%have_south_pole
242 have_north_pole => gridstruct%have_north_pole
244 ntiles_g => gridstruct%ntiles_g
245 acapn => gridstruct%acapN
246 acaps => gridstruct%acapS
247 globalarea => gridstruct%globalarea
265 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)
271 psi(i,j) = (-1.0 * ubar *
radius *( sin(agrid(i,j,2)) *cos(
alpha) - &
272 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
278 psi_b(i,j) = (-1.0 * ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
279 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
283 if ( (cubed_sphere) .and. (defongrid==0) )
then 287 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
288 if (dist==0) vc(i,j) = 0.
294 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
295 if (dist==0) uc(i,j) = 0.
299 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
303 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
304 if (dist==0) v(i,j) = 0.
310 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
311 if (dist==0) u(i,j) = 0.
317 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
318 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
320 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
321 if (dist==0) ua(i,j) = 0.
322 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
323 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
325 va(i,j) = (psi2 - psi1) / (dist)
326 if (dist==0) va(i,j) = 0.
330 elseif ( (cubed_sphere) .and. (defongrid==1) )
then 334 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
335 if (dist==0) vc(i,j) = 0.
341 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
342 if (dist==0) uc(i,j) = 0.
346 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
347 call ctoa(uc,vc,ua,va,dx, dy, dxc,dyc,dxa,dya,npx,npy,ng)
348 call atod(ua,va,u ,v ,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
351 elseif ( (cubed_sphere) .and. (defongrid==2) )
then 355 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
356 if (dist==0) v(i,j) = 0.
362 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
363 if (dist==0) u(i,j) = 0.
367 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng)
368 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
369 elseif ( (cubed_sphere) .and. (defongrid==3) )
then 372 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
373 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
375 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
376 if (dist==0) ua(i,j) = 0.
377 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
378 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
380 va(i,j) = (psi2 - psi1) / (dist)
381 if (dist==0) va(i,j) = 0.
385 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
386 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested,domain)
387 elseif ( (latlon) .or. (defongrid==4) )
then 391 ua(i,j) = ubar * ( cos(agrid(i,j,2))*cos(
alpha) + &
392 sin(agrid(i,j,2))*cos(agrid(i,j,1))*sin(
alpha) )
393 va(i,j) = -ubar * sin(agrid(i,j,1))*sin(
alpha)
398 if (cubed_sphere)
call rotate_winds(ua(i,j), va(i,j), p1,p2,p3,p4, agrid(i,j,1:2), 2, 1)
400 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
401 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
403 if ( (tile==1) .and.(i==1) ) print*, ua(i,j), -1.0 * (psi2 - psi1) / (dist)
408 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, nested, domain)
409 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
410 elseif ( (latlon) .or. (defongrid==5) )
then 415 p1(:) = grid(i ,j ,1:2)
416 p2(:) = grid(i,j+1 ,1:2)
420 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
421 sin(pt(2))*cos(pt(1))*sin(
alpha) )
422 vtmp = -ubar * sin(pt(1))*sin(
alpha)
429 p1(:) = grid(i ,j ,1:2)
430 p2(:) = grid(i+1,j ,1:2)
434 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
435 sin(pt(2))*cos(pt(1))*sin(
alpha) )
436 vtmp = -ubar * sin(pt(1))*sin(
alpha)
442 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng)
443 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, nested, domain)
464 subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
465 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, &
466 dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, adiabatic, &
467 ks, npx_global, ptop, domain_in, tile_in, bd)
470 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
471 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
472 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
473 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
474 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
475 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
477 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
479 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
480 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
481 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
482 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
483 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
485 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
486 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
487 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
488 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
489 real ,
intent(inout) :: delz(bd%isd:,bd%jsd:,1:)
490 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
492 real ,
intent(inout) :: ak(npz+1)
493 real ,
intent(inout) :: bk(npz+1)
495 integer,
intent(IN) :: npx, npy, npz
496 integer,
intent(IN) ::
ng, ncnst, nwat
497 integer,
intent(IN) :: ndims
498 integer,
intent(IN) :: nregions
500 real,
intent(IN) :: dry_mass
501 logical,
intent(IN) :: mountain
502 logical,
intent(IN) :: moist_phys
503 logical,
intent(IN) :: hydrostatic
504 logical,
intent(IN) :: hybrid_z
505 logical,
intent(IN) :: adiabatic
506 integer,
intent(IN) :: ks
511 integer,
intent(IN) :: npx_global
512 integer,
intent(IN),
target :: tile_in
513 real,
intent(INOUT) :: ptop
515 type(
domain2d),
intent(IN),
target :: domain_in
517 real :: tmp(1-
ng:npx +
ng,1-
ng:npy +
ng,1:nregions)
518 real :: tmp1(1 :npx ,1 :npy ,1:nregions)
520 real(kind=R_GRID) :: p0(2)
521 real(kind=R_GRID) :: p1(2)
522 real(kind=R_GRID) :: p2(2)
523 real(kind=R_GRID) :: p3(2)
524 real(kind=R_GRID) :: p4(2)
525 real(kind=R_GRID) :: pa(2)
526 real(kind=R_GRID) :: pb(2)
527 real(kind=R_GRID) :: pcen(2)
528 real(kind=R_GRID) :: e1(3), e2(3), e3(3), ex(3), ey(3)
529 real :: dist, r, r1, r2, r0, omg, a, b, c
530 integer :: i,j,k,nreg,z,zz
531 integer :: i0,j0,n0, nt
532 real :: utmp,vtmp,ftmp
535 integer,
parameter :: jm = 5761
542 real :: ddeg, deg, ddp, dp, ph5
547 real :: x1,y1,z1,x2,y2,z2,ang
549 integer :: initwindscase
561 real :: grad(bd%isd:bd%ied ,bd%jsd:bd%jed,2)
562 real :: div0(bd%isd:bd%ied ,bd%jsd:bd%jed )
563 real :: vor0(bd%isd:bd%ied ,bd%jsd:bd%jed )
564 real :: divg(bd%isd:bd%ied ,bd%jsd:bd%jed )
565 real :: vort(bd%isd:bd%ied ,bd%jsd:bd%jed )
566 real :: ztop, rgrav, p00, pturb, zmid, pk0, t00
567 real :: dz1(npz), ppt(npz)
568 real :: ze1(npz+1), pe1(npz+1)
571 character(len=80) :: oflnm, hgtflnm
572 integer :: is2, ie2, js2, je2
574 real :: psi(bd%isd:bd%ied,bd%jsd:bd%jed)
575 real :: psi_b(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
579 real :: eta(npz), eta_0, eta_s, eta_t
580 real :: eta_v(npz), press, anti_rot
581 real :: t_0, t_mean, delta_t, lapse_rate, n2, zeta, s0
582 real :: pt1,pt2,pt3,pt4,pt5,pt6, pt7, pt8, pt9,
u1, pt0
583 real :: uu1, uu2, uu3, vv1, vv2, vv3
586 real wbuffer(npy+2,npz)
587 real sbuffer(npx+2,npz)
589 real :: gz(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1), zt, zdist
596 real,
dimension(npz):: pk1, ts1, qs1, uz1, zs1, dudz
598 real(kind=R_GRID):: pp0(2)
603 real :: omega0, k_cell, z0, h, px
604 real :: d1, d2, p1p(2), rt, s
605 real :: wind_alpha, period, h0, rm, zp3(3), dz3(3), k0, lp
609 real,
dimension(npz+1) :: pe0, gz0, ue, ve, we, pte, qe
610 real :: d, cor, exppr, exppz, gamma, ts0, q00, exponent, ztrop, height, zp, rp
611 real :: qtrop, ttrop, zq1, zq2
612 real :: dum, dum1, dum2, dum3, dum4, dum5, dum6, ptmp, uetmp, vetmp
613 real :: pe_u(bd%is:bd%ie,npz+1,bd%js:bd%je+1)
614 real :: pe_v(bd%is:bd%ie+1,npz+1,bd%js:bd%je)
615 real :: ps_u(bd%is:bd%ie,bd%js:bd%je+1)
616 real :: ps_v(bd%is:bd%ie+1,bd%js:bd%je)
621 real(kind=R_GRID),
pointer,
dimension(:,:,:) ::
agrid, grid
622 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
623 real,
pointer,
dimension(:,:) :: rarea, fc,
f0 624 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
625 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
626 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
628 logical,
pointer :: cubed_sphere, latlon
631 integer,
pointer :: tile
633 logical,
pointer :: have_south_pole, have_north_pole
635 integer,
pointer :: ntiles_g
636 real,
pointer :: acapn, acaps, globalarea
647 grid => gridstruct%grid_64
648 agrid=> gridstruct%agrid_64
650 area => gridstruct%area_64
651 rarea => gridstruct%rarea
656 ee1 => gridstruct%ee1
657 ee2 => gridstruct%ee2
660 en1 => gridstruct%en1
661 en2 => gridstruct%en2
665 dxa => gridstruct%dxa
666 dya => gridstruct%dya
667 rdxa => gridstruct%rdxa
668 rdya => gridstruct%rdya
669 dxc => gridstruct%dxc
670 dyc => gridstruct%dyc
672 cubed_sphere => gridstruct%cubed_sphere
673 latlon => gridstruct%latlon
678 have_south_pole => gridstruct%have_south_pole
679 have_north_pole => gridstruct%have_north_pole
681 ntiles_g => gridstruct%ntiles_g
682 acapn => gridstruct%acapN
683 acaps => gridstruct%acapS
684 globalarea => gridstruct%globalarea
686 if (gridstruct%nested)
then 700 f0(:,:) = huge(dummy)
701 fc(:,:) = huge(dummy)
704 fc(i,j) = 2.*
omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
705 sin(grid(i,j,2))*cos(
alpha) )
715 if (cubed_sphere)
call fill_corners(
f0, npx, npy, ydir)
717 delp(isd:is-1,jsd:js-1,1:npz)=0.
718 delp(isd:is-1,je+1:jed,1:npz)=0.
719 delp(ie+1:ied,jsd:js-1,1:npz)=0.
720 delp(ie+1:ied,je+1:jed,1:npz)=0.
722 #if defined(SW_DYNAMICS) 736 call init_winds(
ubar, u,v,ua,va,uc,vc, 1, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
741 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
742 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
743 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)
749 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
750 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
759 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
760 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)
761 201
format(
' ',a,e21.14,
' ',e21.14)
762 202
format(
' ',a,i4.4,
'x',i4.4,
'x',i4.4)
763 if ( is_master() )
then 764 write(*,*)
' Error Norms of Analytical Divergence field C-Winds initialized' 765 write(*,201)
'Divergence MAX error : ', pmax
766 write(*,201)
'Divergence MIN error : ', pmin
767 write(*,201)
'Divergence L1_norm : ', l1_norm
768 write(*,201)
'Divergence L2_norm : ', l2_norm
769 write(*,201)
'Divergence Linf_norm : ', linf_norm
772 call init_winds(
ubar, u,v,ua,va,uc,vc, 3, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
776 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
777 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
778 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)
784 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
785 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
792 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
793 if ( is_master() )
then 794 write(*,*)
' Error Norms of Analytical Divergence field A-Winds initialized' 795 write(*,201)
'Divergence MAX error : ', pmax
796 write(*,201)
'Divergence MIN error : ', pmin
797 write(*,201)
'Divergence L1_norm : ', l1_norm
798 write(*,201)
'Divergence L2_norm : ', l2_norm
799 write(*,201)
'Divergence Linf_norm : ', linf_norm
802 call init_winds(
ubar, u,v,ua,va,uc,vc, 2, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
808 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
809 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
810 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)
816 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
817 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
822 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
823 if ( is_master() )
then 824 write(*,*)
' Error Norms of Analytical Divergence field D-Winds initialized' 825 write(*,201)
'Divergence MAX error : ', pmax
826 write(*,201)
'Divergence MIN error : ', pmin
827 write(*,201)
'Divergence L1_norm : ', l1_norm
828 write(*,201)
'Divergence L2_norm : ', l2_norm
829 write(*,201)
'Divergence Linf_norm : ', linf_norm
843 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
845 if (p /= 0.0) w_p = vtx/p
846 delp(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*0.0) )
849 ua(i,j,1) = ua(i,j,1)*
radius/86400.0
850 va(i,j,1) = va(i,j,1)*
radius/86400.0
856 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)
861 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,
ng, gridstruct%nested, domain)
863 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
865 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
880 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(
pi*r/r0))
882 delp(i,j,1) = phis(i,j)
903 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.))
then 924 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
925 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0
927 ( -1.*cos(grid(i+1,j ,1))*cos(grid(i+1,j ,2))*sin(
alpha) + &
928 sin(grid(i+1,j ,2))*cos(
alpha) ) ** 2.0
930 ( -1.*cos(grid(i+1,j+1,1))*cos(grid(i+1,j+1,2))*sin(
alpha) + &
931 sin(grid(i+1,j+1,2))*cos(
alpha) ) ** 2.0
933 ( -1.*cos(grid(i,j+1,1))*cos(grid(i,j+1,2))*sin(
alpha) + &
934 sin(grid(i,j+1,2))*cos(
alpha) ) ** 2.0
935 delp(i,j,1) = (0.25*(pt1+pt2+pt3+pt4) + 3.*pt5) / 4.
964 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(
pi*r/r0))
966 delp(i,j,1) = phis(i,j)
969 delp(i,j,1) = delp(i,j,1) +
grav*2.e3
980 p1(:) = grid(i ,j ,1:2)
981 p2(:) = grid(i,j+1 ,1:2)
985 utmp =
ubar * cos(p3(2))
992 p1(:) = grid(i, j,1:2)
993 p2(:) = grid(i+1,j,1:2)
997 utmp =
ubar * cos(p3(2))
1006 fc(i,j) = 2.*anti_rot*sin(grid(i,j,2))
1011 f0(i,j) = 2.*anti_rot*sin(
agrid(i,j,2))
1040 p1(1) =
pi*1.5 - ddeg
1044 p2(1) =
pi*1.5 + ddeg
1048 #ifndef SINGULAR_VORTEX 1080 p2(1) =
agrid(i,j,1)
1081 p2(2) =
agrid(i,j,2)
1082 r =
min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1084 phis(i,j) = 2000.0*
grav*(1.0-(r/r0))
1091 sin(
agrid(i ,j ,2))*cos(
alpha) ) ** 2 - phis(i,j)
1103 a = 0.5*omg*(2.*
omega+omg)*(cos(
agrid(i,j,2))**2) + &
1104 0.25*rk*rk*(cos(
agrid(i,j,2))**(r+r)) * &
1105 ( (r+1)*(cos(
agrid(i,j,2))**2) + (2.*r*r-r-2.) - &
1106 2.*(r*r)*cos(
agrid(i,j,2))**(-2.) )
1107 b = (2.*(
omega+omg)*rk / ((r+1)*(r+2))) * (cos(
agrid(i,j,2))**r) * &
1108 ( (r*r+2.*r+2.) - ((r+1.)*cos(
agrid(i,j,2)))**2 )
1109 c = 0.25*rk*rk*(cos(
agrid(i,j,2))**(2.*r)) * ( &
1110 (r+1) * (cos(
agrid(i,j,2))**2.) - (r+2.) )
1112 delp(i,j,1) = delp(i,j,1) - phis(i,j)
1117 p1(:) = grid(i ,j ,1:2)
1118 p2(:) = grid(i,j+1 ,1:2)
1122 utmp =
radius*omg*cos(p3(2)) + &
1123 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1124 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1130 p1(:) = grid(i, j,1:2)
1131 p2(:) = grid(i+1,j,1:2)
1135 utmp =
radius*omg*cos(p3(2)) + &
1136 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1137 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1142 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,
ng)
1144 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
1171 pt6 =
gh_jet(npy, grid(i, j, 2))
1172 pt7 =
gh_jet(npy, grid(i+1,j, 2))
1173 pt8 =
gh_jet(npy, grid(i+1,j+1,2))
1174 pt9 =
gh_jet(npy, grid(i ,j+1,2))
1175 ftmp = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1177 delp(i,j,1) = ftmp + 120.*
grav*cos(
agrid(i,j,2)) * &
1178 exp( -(3.*(
agrid(i,j,1)-
pi))**2 ) * exp( -(15.*(
agrid(i,j,2)-
pi/4.))**2 )
1184 p1(:) =
agrid(i,j,1:2)
1187 if ( r < 3.*r0 )
then 1188 delp(i,j,1) = delp(i,j,1) + 1000.*
grav*exp(-(r/r0)**2)
1197 p2(:) = grid(i,j+1,1:2)
1198 vv1 =
u_jet(p2(2))*(ee2(2,i,j+1)*cos(p2(1)) - ee2(1,i,j+1)*sin(p2(1)))
1199 p1(:) = grid(i,j,1:2)
1200 vv3 =
u_jet(p1(2))*(ee2(2,i,j)*cos(p1(1)) - ee2(1,i,j)*sin(p1(1)))
1203 vv2 =
u_jet(pa(2))*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1205 v(i,j,1) = 0.25*(vv1 + 2.*vv2 + vv3)
1212 p1(:) = grid(i,j,1:2)
1213 uu1 =
u_jet(p1(2))*(ee1(2,i,j)*cos(p1(1)) - ee1(1,i,j)*sin(p1(1)))
1214 p2(:) = grid(i+1,j,1:2)
1215 uu3 =
u_jet(p2(2))*(ee1(2,i+1,j)*cos(p2(1)) - ee1(1,i+1,j)*sin(p2(1)))
1218 uu2 =
u_jet(pa(2))*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1220 u(i,j,1) = 0.25*(uu1 + 2.*uu2 + uu3)
1227 call get_vorticity(is, ie, js, je, isd, ied, jsd, jed, npz, u, v, q(is:ie,js:je,:,1), dx, dy, rarea)
1230 fc(i,j) = 2.*
omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
1231 sin(grid(i,j,2))*cos(
alpha) )
1241 if (cubed_sphere)
call fill_corners(
f0, npx, npy, ydir)
1244 q(i,j,npz,1) = ( q(i,j,npz,1) +
f0(i,j) ) / delp(i,j,npz) * 1.e6
1262 p2(1) =
agrid(i,j,1)
1263 p2(2) =
agrid(i,j,2)
1264 r =
min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1266 phis(i,j) = 2000.0*
grav*(1.0-(r/r0))
1281 if ( is_master() )
write(*,*)
'Initialzing case-8: soliton twin cycolne...' 1304 p1(:) = grid(i ,j ,1:2)
1305 p2(:) = grid(i,j+1 ,1:2)
1308 utmp =
ubar*exp(-(r/r0)**2)
1316 p1(:) = grid(i, j,1:2)
1317 p2(:) = grid(i+1,j,1:2)
1320 utmp =
ubar*exp(-(r/r0)**2)
1333 p1(:) = grid(i ,j ,1:2)
1334 p2(:) = grid(i,j+1 ,1:2)
1337 utmp =
ubar*exp(-(r/r0)**2)
1345 p1(:) = grid(i, j,1:2)
1346 p2(:) = grid(i+1,j,1:2)
1349 utmp =
ubar*exp(-(r/r0)**2)
1364 ph5 = -0.5*
pi + (dble(j-1)-0.5)*ddp
1365 ll_j(j) = -0.5*
pi + (dble(j-1)*ddp)
1371 cosp(j) = (sine(j+1)-sine(j)) / dp
1374 cose(j) = 0.5 * (cosp(j-1) + cosp(j))
1377 ddeg = 180./float(jm-1)
1379 deg = -90. + (float(j-1)-0.5)*ddeg
1381 ll_u(j) = -10.*(deg+90.)/90.
1382 elseif (deg <= 60.)
then 1383 ll_u(j) = -10. + deg
1385 ll_u(j) = 50. - (50./30.)* (deg - 60.)
1388 ll_phi(1) = 6000. *
grav 1390 ll_phi(j)=ll_phi(j-1) - dp*sine(j) * &
1397 if ( (ll_j(jj) <=
agrid(i,j,2)) .and. (
agrid(i,j,2) <= ll_j(jj+1)) )
then 1398 delp(i,j,1)=0.5*(ll_phi(jj)+ll_phi(jj+1))
1407 ua(i,j,1) = -10.*(
agrid(i,j,2)*
todeg + 90.)/90.
1411 ua(i,j,1) = 50. - (50./30.)* (
agrid(i,j,2)*
todeg - 60.)
1418 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)
1423 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
1425 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
1426 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,
ng, gridstruct%nested, domain)
1437 if ( is_master() )
write(*,*)
'Initialzing case-9: soliton cyclones...' 1459 p1(:) = grid(i ,j ,1:2)
1460 p2(:) = grid(i,j+1 ,1:2)
1463 utmp =
ubar*exp(-(r/r0)**2)
1471 p1(:) = grid(i, j,1:2)
1472 p2(:) = grid(i+1,j,1:2)
1475 utmp =
ubar*exp(-(r/r0)**2)
1488 delp(:,:,z) = delp(:,:,1)
1495 call init_winds(
ubar, u,v,ua,va,uc,vc, initwindscase, npx, npy,
ng, ndims, nregions, gridstruct%nested, gridstruct, domain, tile)
1504 ps(i,j) = delp(i,j,1)
1518 if (.not.hydrostatic) w(:,:,:)= 0.0
1523 call hydro_eq(npz, is, ie, js, je, ps, phis, 1.e5, &
1524 delp, ak, bk, pt, delz, area,
ng, .false., hydrostatic, hybrid_z, domain)
1531 p1(2) =
pi/6. + (7.5/180.0)*
pi 1534 p2(1) =
agrid(i,j,1)
1535 p2(2) =
agrid(i,j,2)
1536 r =
min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1538 phis(i,j) =
gh0*(1.0-(r/r0))
1541 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1542 delp, ak, bk, pt, delz, area,
ng, mountain, hydrostatic, hybrid_z, domain)
1546 call surfdrv(npx, npy, gridstruct%grid_64, gridstruct%agrid_64, &
1547 gridstruct%area_64, dx, dy, dxa, dya, dxc, dyc, &
1548 gridstruct%sin_sg, phis, &
1549 flagstruct%stretch_fac, gridstruct%nested, &
1550 npx_global, domain, flagstruct%grid_number, bd)
1553 if ( hybrid_z )
then 1559 if ( is_master() )
write(*,*)
'Using const DZ' 1561 dz1(1) = ztop /
real(npz) 1562 dz1(npz) = 0.5*dz1(1)
1581 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1582 delp, ak, bk, pt, delz, area,
ng, mountain, hydrostatic, hybrid_z, domain)
1587 if (is_master()) print*,
'TEST TRACER enabled for this test case' 1590 ncnst, npz, q,
agrid(is:ie,js:je,1),
agrid(is:ie,js:je,2), 9., 9.)
1597 p1(1) = 205.*
pi/180.
1601 p2(1) =
agrid(i,j,1)
1602 p2(2) =
agrid(i,j,2)
1604 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.) .and. k > 16)
then 1622 if (cl > 0 .and. cl2 > 0)
then 1624 q, delp,ncnst,
agrid(isd:ied,jsd:jed,1),
agrid(isd:ied,jsd:jed,2))
1637 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
1646 peln(i,1,j) = log(ptop)
1647 pk(i,j,1) = ptop**
kappa 1652 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
1653 pk(i,j,k) = exp(
kappa*log(pe(i,k,j)) )
1654 peln(i,k,j) = log(pe(i,k,j))
1663 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
1671 eta(k) = 0.5*( (ak(k)+ak(k+1))/1.e5 + bk(k)+bk(k+1) )
1672 eta_v(k) = (eta(k) - eta_0)*
pi*0.5
1675 if ( .not. adiabatic )
then 1688 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j)) - 100000.
1689 q(i,j,k,
sphum) = 0.021*exp(-(
agrid(i,j,2)/pcen(2))**4.)*exp(-(ptmp/34000.)**2.)
1718 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j+1,2))**2.0
1721 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1722 vv1 = utmp*(ee2(2,i,j+1)*cos(grid(i,j+1,1)) - ee2(1,i,j+1)*sin(grid(i,j+1,1)))
1724 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1727 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1728 vv3 = utmp*(ee2(2,i,j)*cos(grid(i,j,1)) - ee2(1,i,j)*sin(grid(i,j,1)))
1730 p1(:) = grid(i ,j ,1:2)
1731 p2(:) = grid(i,j+1 ,1:2)
1733 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1736 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1737 vv2 = utmp*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1739 v(i,j,z) = 0.25*(vv1 + 2.*vv2 + vv3)
1744 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1747 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1748 uu1 = utmp*(ee1(2,i,j)*cos(grid(i,j,1)) - ee1(1,i,j)*sin(grid(i,j,1)))
1750 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i+1,j,2))**2.0
1753 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1754 uu3 = utmp*(ee1(2,i+1,j)*cos(grid(i+1,j,1)) - ee1(1,i+1,j)*sin(grid(i+1,j,1)))
1756 p1(:) = grid(i ,j ,1:2)
1757 p2(:) = grid(i+1,j ,1:2)
1759 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1762 if (-(r/r0)**2.0 > -40.0) utmp = utmp +
u1*exp(-(r/r0)**2.0)
1763 uu2 = utmp*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1765 u(i,j,z) = 0.25*(uu1 + 2.*uu2 + uu3)
1780 eta(z) = 0.5*( (ak(z)+ak(z+1))/1.e5 + bk(z)+bk(z+1) )
1782 t_mean = t_0 * eta(z)**(
rdgas*lapse_rate/
grav)
1783 if (eta_t > eta(z)) t_mean = t_mean + delta_t*(eta_t - eta(z))**5.0
1785 230
format(i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
1788 press = press + delp(is,js,zz)
1790 if (is_master())
write(*,230) z, eta(z), press/100., t_mean
1794 pt1 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1795 ( -2.0*(sin(
agrid(i,j,2))**6.0) *(cos(
agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1796 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1808 pt2 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1809 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1810 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1811 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1813 pt3 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1814 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1815 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1816 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1818 pt4 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1819 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1820 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1821 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1823 pt5 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1824 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1825 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1826 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1828 pt6 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1829 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1830 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1831 ( (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 )
1832 pt7 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1833 ( -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 ) * &
1834 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1835 ( (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 )
1836 pt8 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1837 ( -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 ) * &
1838 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1839 ( (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 )
1840 pt9 = t_mean + 0.75*(eta(z)*
pi*
ubar/
rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1841 ( -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 ) * &
1842 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1843 ( (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 )
1844 pt(i,j,z) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1851 if ( (r/r0)**2 < 40. )
then 1852 pt(i,j,z) = pt(i,j,z) + pt0*exp(-(r/r0)**2)
1859 if (is_master()) print*,
' ' 1866 pt1 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1867 ( -2.0*(sin(
agrid(i,j,2))**6.0) *(cos(
agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1868 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1880 pt2 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1881 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1882 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1883 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1885 pt3 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1886 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1887 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1888 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1890 pt4 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1891 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1892 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1893 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1895 pt5 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1896 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1897 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1898 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) -
pi/4.0 )*
radius*
omega )
1900 pt6 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1901 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1902 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1903 ( (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 )
1904 pt7 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1905 ( -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 ) * &
1906 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1907 ( (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 )
1908 pt8 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1909 ( -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 ) * &
1910 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1911 ( (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 )
1912 pt9 =
ubar* (cos( (eta_s-eta_0)*
pi/2.0 ))**(3.0/2.0) * ( &
1913 ( -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 ) * &
1914 ubar*cos( (eta_s-eta_0)*
pi/2.0 )**(3.0/2.0) + &
1915 ( (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 )
1916 phis(i,j) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1923 if ( .not.hydrostatic )
then 1929 delz(i,j,k) =
rdgas/
grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
1935 if (.not. adiabatic)
then 1941 pt(i,j,k) = pt(i,j,k)/(1. + zvir*q(i,j,k,
sphum))
1951 write(stdout(), *)
'PI:',
pi 1952 write(stdout(), *)
'PHIS:',
mpp_chksum(phis(is:ie,js:je))
1957 is,ie,js,je,isd,ied,jsd,jed,npz,ncnst,ak,bk,ptop, &
1958 pk,peln,pe,pkz,gz,phis,ps,grid,
agrid,hydrostatic, &
1959 nwat, adiabatic,
test_case == -13, domain)
1961 write(stdout(), *)
'PHIS:',
mpp_chksum(phis(is:ie,js:je))
1987 ze1(k) = ze1(k+1) + ztop/
real(npz)
1991 ze1(1) = ztop + 1.5*ztop/
real(npz)
2004 delz(i,j,k) = ze1(k+1) - ze1(k)
2005 pk(i,j,k) = pk(i,j,k+1) +
grav*delz(i,j,k)/(
cp_air*t00)*pk0
2006 pe(i,k,j) = pk(i,j,k)**(1./
kappa)
2012 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
2017 peln(i,k,j) = log(pe(i,k,j))
2026 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
2027 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2040 r0 = 0.5*(ze1(k)+ze1(k+1)) - 3.2e3
2042 r0 = (0.5*(ze1(k)+ze1(k+1)) - 3.0e3) / 2.e3
2047 p2(1) =
agrid(i,j,1)
2048 p2(2) =
agrid(i,j,2)
2051 dist = sqrt( r**2 + r0**2 ) / 3.2e3
2054 dist = sqrt( r**2 + r0**2 )
2056 if ( dist<=1. )
then 2057 q(i,j,k,1) = pk0 * pturb/pkz(i,j,k)*(cos(
pi*dist)+1.)/2.
2058 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(
pi*dist)+1.)/2.
2063 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
2080 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2083 pe1(z) = ak(z) + bk(z)*p00
2088 ze1(z) = ze1(z+1) + ztop/
real(npz)
2092 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2096 ps(i,j) = pe1(npz+1)
2104 peln(i,z,j) = log(pe1(z))
2105 pk(i,j,z) = exp(
kappa*peln(i,z,j))
2118 vort(i,j) = 0.5*(1.+cos(
pi*r/r0))
2129 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*
pi/ztop )
2132 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(
kappa*(peln(i,z+1,j)-peln(i,z,j)))
2133 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2135 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2136 q(i,j,z,1) = q(i,j,z,1) + vort(i,j)*zmid
2149 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2152 pe1(z) = ak(z) + bk(z)*p00
2157 ze1(z) = ze1(z+1) + ztop/
real(npz)
2161 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2165 ps(i,j) = pe1(npz+1)
2173 peln(i,z,j) = log(pe1(z))
2174 pk(i,j,z) = exp(
kappa*peln(i,z,j))
2187 vort(i,j) = 0.5*(1.+cos(
pi*r/r0))
2197 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*
pi/ztop )
2200 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(
kappa*(peln(i,z+1,j)-peln(i,z,j)))
2201 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2203 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2220 phis(i,j) =
grav*2.e3*exp( -(r/1500.e3)**2 )
2230 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
2236 p1(:) = grid(i ,j ,1:2)
2237 p2(:) = grid(i,j+1 ,1:2)
2241 utmp =
ubar * cos(p3(2))
2250 p1(:) = grid(i, j,1:2)
2251 p2(:) = grid(i+1,j,1:2)
2255 utmp =
ubar * cos(p3(2))
2281 p1(:) = grid(i ,j ,1:2)
2282 p2(:) = grid(i,j+1 ,1:2)
2286 utmp =
ubar * cos(p3(2))
2293 p1(:) = grid(i, j,1:2)
2294 p2(:) = grid(i+1,j,1:2)
2298 utmp =
ubar * cos(p3(2))
2319 p1(1) = (0.5-0.125) *
pi 2326 p2(:) =
agrid(i,j,1:2)
2329 p4(:) = p2(:) - p1(:)
2330 if ( abs(p4(1)) > 1.e-12 )
then 2331 zeta = asin( p4(2) / sqrt(p4(1)**2 + p4(2)**2) )
2335 if ( p4(1) <= 0. ) zeta =
pi - zeta
2337 v1 = r/uu1 * cos( zeta )
2338 v2 = r/uu2 * sin( zeta )
2339 phis(i,j) = ftop / ( 1. + v1**2 + v2**2 )
2346 if ( hybrid_z )
then 2350 elseif( npz.eq.31 .or. npz.eq.41 .or. npz.eq.51 )
then 2354 if ( is_master() )
write(*,*)
'Using const DZ' 2356 dz1(1) = ztop /
real(npz) 2361 dz1(1) =
max( 1.0e3, 3.*dz1(2) )
2367 ze1(k) = ze1(k+1) + dz1(k)
2374 call mpp_error(fatal,
'This test case is only currently setup for hybrid_z')
2380 delz(i,j,k) = ze0(i,j,k+1) - ze0(i,j,k)
2394 pe1(k) = p00*( (1.-s0/t00) + s0/t00*exp(-n2*ze1(k)/
grav) )**(1./
kappa)
2398 if ( is_master() )
write(*,*)
'Lee vortex testcase: model top (mb)=', ptop/100.
2404 bk(k) = (pe1(k) - pe1(1)) / (pe1(npz+1)-pe1(1))
2405 ak(k) = pe1(1)*(1.-bk(k))
2415 pe(i,k,j) = pk(i,j,k) ** (1./
kappa)
2416 peln(i,k,j) = log(pe(i,k,j))
2424 peln(i,1,j) = log(pe(i,1,j))
2425 pk(i,j,1) = pe(i,1,j) **
kappa 2426 ps(i,j) = pe(i,npz+1,j)
2433 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
2434 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2435 pt(i,j,k) = pkz(i,j,k)*
grav*delz(i,j,k) / (
cp_air*(pk(i,j,k)-pk(i,j,k+1)) )
2446 if (.not.hydrostatic) w(:,:,:)= 0.0
2459 dz = 12000./
real(npz)
2461 allocate(
zz0(npz+1))
2462 allocate(
pz0(npz+1))
2470 if (is_master()) print*,
'TRACER ADVECTION TEST CASE' 2471 if (is_master()) print*,
'INITIAL LEVELS' 2477 if (is_master())
write(*,*) k,
pz0(k),
zz0(k)
2490 ptmp = 0.5*(
pz0(k) +
pz0(k+1))
2497 delp(i,j,k) =
pz0(k+1)-
pz0(k)
2502 ptop = 100000.*exp(-12000.*
grav/t00/
rdgas)
2516 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
2517 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
2525 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
2526 if (dist==0) vc(i,j,k) = 0.
2532 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
2533 if (dist==0) uc(i,j,k) = 0.
2540 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
2541 if (dist==0) v(i,j,k) = 0.
2547 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
2548 if (dist==0) u(i,j,k) = 0.
2554 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
2555 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
2557 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
2558 if (dist==0) ua(i,j,k) = 0.
2559 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
2560 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
2562 va(i,j,k) = (psi2 - psi1) / (dist)
2563 if (dist==0) va(i,j,k) = 0.
2570 uc(:,:,k) = uc(:,:,1)
2571 vc(:,:,k) = vc(:,:,1)
2572 ua(:,:,k) = ua(:,:,1)
2573 va(:,:,k) = va(:,:,1)
2577 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
2585 call mpp_error(fatal,
'Value of tracer_test not implemented ')
2603 if (.not.hydrostatic) w(:,:,:)= 0.0
2607 dz = 12000./
real(npz)
2613 px = ((t00-9000.*gamma)/t00)**(1./exponent)
2617 height = 12000. - dz*
real(k-1)
2618 if (height >= 9000. )
then 2619 ak(k) = p00*((t00-height*gamma)/t00)**(1./exponent)
2622 ak(k) = (((t00-height*gamma)/t00)**(1./exponent)-1.)/(px - 1.)*px*p00
2623 bk(k) = (((t00-height*gamma)/t00)**(1./exponent)-px)/(1.-px)
2625 if (is_master())
write(*,*) k, ak(k), bk(k), height, ak(k)+bk(k)*p00
2631 p1(1) = 3.*
pi/2. ; p1(2) = 0.
2638 p2(:) =
agrid(i,j,1:2)
2641 phis(i,j) =
grav*0.5*2000.*(1. + cos(
pi*r/r0))*cos(
pi*r/zetam)**2.
2642 pe(i,npz+1,j) = p00*(1.-gamma/t00*phis(i,j)/
grav)**(1./exponent)
2647 ps(i,j) = pe(i,npz+1,j)
2654 pe(i,k,j) = ak(k) + bk(k)*ps(i,j)
2655 gz(i,j,k) = t00/gamma*(1. - (pe(i,k,j)/p00)**exponent)
2667 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
2672 pt(i,j,k) = -
grav*t00*p00/(
rdgas*gamma +
grav)/delp(i,j,k) * &
2673 ( (pe(i,k,j)/p00)**(exponent+1.) - (pe(i,k+1,j)/p00)**(exponent+1.) )
2695 pk(i,j,1) = ptop**
kappa 2697 peln(i,1,j) = log(ptop)
2704 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
2705 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
2706 peln(i,k+1,j) = log(pe(i,k+1,j))
2707 pk(i,j,k+1) = exp(
kappa*peln(i,k+1,j) )
2715 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
2722 pp0(1) = 262.0/180.*
pi 2723 pp0(2) = 35.0/180.*
pi 2730 delz(i,j,k) =
rdgas/
grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
2737 ze1(k) = ze1(k+1) - delz(is,js,k)
2741 if (is_master())
then 2743 write(6,*)
'Toy supercell winds, piecewise approximation' 2745 write(6,*)
'Toy supercell winds, tanh approximation' 2750 zm = 0.5*(ze1(k)+ze1(k+1))
2755 if ( zm .le. 2.e3 )
then 2756 utmp = 8.*(1.-cos(
pi*zm/4.e3))
2757 vtmp = 8.*sin(
pi*zm/4.e3)
2758 elseif (zm .le. 6.e3 )
then 2759 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
2769 utmp = 15.0*(1.+tanh(zm/2000. - 1.5))
2770 vtmp = 8.5*tanh(zm/1000.)
2785 if( is_master() )
then 2786 write(6,*) k, utmp, vtmp
2791 p1(:) = grid(i ,j ,1:2)
2792 p2(:) = grid(i,j+1 ,1:2)
2803 p1(:) = grid(i, j,1:2)
2804 p2(:) = grid(i+1,j,1:2)
2815 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
2816 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
2817 .true., hydrostatic, nwat, domain)
2824 zm = 0.5*(ze1(k)+ze1(k+1))
2825 ptmp = ( (zm-zc)/zc ) **2
2826 if ( ptmp < 1. )
then 2830 if ( dist < 1. )
then 2831 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
2840 call mpp_error(fatal,
' test_case 32 not yet implemented')
2866 p0(1) = 60./180. *
pi 2881 p2(2) =
agrid(i,j,2)
2882 pt1 = exp( dum*(us0*sin(p2(2)))**2 )
2884 pt2 = exp( dum*(us0*sin(p2(2)))**2 )
2886 pt3 = exp( dum*(us0*sin(p2(2)))**2 )
2888 pt4 = exp( dum*(us0*sin(p2(2)))**2 )
2890 pt5 = exp( dum*(us0*sin(p2(2)))**2 )
2892 pt6 = exp( dum*(us0*sin(p2(2)))**2 )
2893 p2(2) = grid(i+1,j,2)
2894 pt7 = exp( dum*(us0*sin(p2(2)))**2 )
2895 p2(2) = grid(i+1,j+1,2)
2896 pt8 = exp( dum*(us0*sin(p2(2)))**2 )
2897 p2(2) = grid(i,j+1,2)
2898 pt9 = exp( dum*(us0*sin(p2(2)))**2 )
2899 ptmp = t00*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
2901 ptmp = t00*exp( dum*(us0*sin(
agrid(i,j,2)))**2 )
2917 p2(1:2) =
agrid(i,j,1:2)
2919 pt1 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2922 pt2 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2925 pt3 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2928 pt4 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2931 pt5 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2932 p2(1:2) = grid(i,j,1:2)
2934 pt6 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2935 p2(1:2) = grid(i+1,j,1:2)
2937 pt7 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2938 p2(1:2) = grid(i+1,j+1,1:2)
2940 pt8 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2941 p2(1:2) = grid(i,j+1,1:2)
2943 pt9 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2944 phis(i,j) =
grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
2946 p2(1:2) =
agrid(i,j,1:2)
2948 phis(i,j) =
grav*h0*cos(p2(2))*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2965 pt1 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2968 pt2 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2971 pt3 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2974 pt4 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2977 pt5 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2979 pt6 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2981 pt7 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2983 pt8 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2985 pt9 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2986 phis(i,j) =
grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
2989 pt1 = exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2990 phis(i,j) =
grav*h0*exp(-(r/5.e3)**2)*cos(
pi*r/4.e3)**2
2999 ps(i,j) = p00*exp( -0.5*(us0*sin(
agrid(i,j,2)))**2/(
rdgas*t00)-phis(i,j)/(
rdgas*pt(i,j,1)) )
3001 peln(i,1,j) = log(ptop)
3002 pk(i,j,1) = ptop**
kappa 3009 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
3010 peln(i,k,j) = log(pe(i,k,j))
3011 pk(i,j,k) = exp(
kappa*peln(i,k,j) )
3019 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3020 delz(i,j,k) =
rdgas/
grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
3028 ze1(npz+1) = phis(i,j)/
grav 3030 ze1(k) = ze1(k+1) - delz(i,j,k)
3033 w(i,j,k) = 0.5*(ze1(k)+ze1(k+1))
3042 p1(:) = grid(i ,j, 1:2)
3043 p2(:) = grid(i,j+1, 1:2)
3048 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i-1,j,k)+w(i,j,k)) )
3054 p1(:) = grid(i, j, 1:2)
3055 p2(:) = grid(i+1,j, 1:2)
3059 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i,j-1,k)+w(i,j,k)) )
3068 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3069 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
3070 .true., hydrostatic, nwat, domain)
3095 ze1(k) = ze1(k+1) + ztop/
real(npz)
3104 call var_dz(npz, ztop, ze1)
3107 zs1(k) = 0.5*(ze1(k)+ze1(k+1))
3115 ts1(k) =
cp_air*ts1(k)*(1.+zvir*qs1(k))
3121 call balanced_k(npz, is, ie, js, je,
ng, pe1(npz+1), ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
3122 delz, zvir, ptop, ak, bk,
agrid)
3125 ps(i,j) = pe(i,npz+1,j)
3132 peln(i,k,j) = log(pe(i,k,j))
3133 pk(i,j,k) = exp(
kappa*peln(i,k,j) )
3141 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3142 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
3151 delz(i,j,k) =
rdgas/
grav*pt(i,j,k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
3158 delz(i,j,k) = ze1(k+1) - ze1(k)
3159 pt(i,j,k) = delz(i,j,k)*
grav/(
rdgas*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j)))
3168 p1(:) = grid(i ,j ,1:2)
3169 p2(:) = grid(i,j+1 ,1:2)
3173 v(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e2,ex)
3178 p1(:) = grid(i, j,1:2)
3179 p2(:) = grid(i+1,j,1:2)
3183 u(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e1,ex)
3200 zm = 0.5*(ze1(k)+ze1(k+1))
3201 ptmp = ( (zm-zc)/zc ) **2
3202 if ( ptmp < 1. )
then 3207 if ( dist < 1. )
then 3208 pt(i,j,k) = pt(i,j,k) + (pkz(i,j,k)/pk0)*pturb*cos(0.5*
pi*dist)**2
3236 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3244 peln(i,1,j) = log(pe(i,1,j))
3245 pk(i,j,1) = exp(
kappa*peln(i,1,j))
3249 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3250 peln(i,k,j) = log(pe(i,k,j))
3251 pk(i,j,k) = exp(
kappa*peln(i,k,j))
3263 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
3265 if ( dist .le. r0 )
then 3277 if (.not.hydrostatic)
then 3281 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))
3310 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
3318 peln(i,1,j) = log(pe(i,1,j))
3319 pk(i,j,1) = exp(
kappa*peln(i,1,j))
3323 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3324 peln(i,k,j) = log(pe(i,k,j))
3325 pk(i,j,k) = exp(
kappa*peln(i,k,j))
3346 p1(:) = grid(i ,j ,1:2)
3347 p2(:) = grid(i,j+1 ,1:2)
3350 utmp =
ubar*exp(-(r/r0)**2)
3358 p1(:) = grid(i, j,1:2)
3359 p2(:) = grid(i+1,j,1:2)
3362 utmp =
ubar*exp(-(r/r0)**2)
3371 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
3373 pt(i,j,k) = pt0/p00**
kappa 3375 pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
3389 q(i,j,k,:) =
agrid(i,j,1)*0.180/
pi 3395 ncnst, npz, q,
agrid(is:ie,js:je,1),
agrid(is:ie,js:je,2), 9., 9.)
3398 if ( .not. hydrostatic )
then 3402 delz(i,j,k) =
rdgas*pt(i,j,k)/
grav*log(pe(i,k,j)/pe(i,k+1,j))
3419 p0(1) = 180. *
pi / 180.
3420 p0(2) = 10. *
pi / 180.
3435 p2(:) =
agrid(i,j,1:2)
3437 ps(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3442 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3448 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3460 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3468 p2(:) = 0.5*(grid(i,j,1:2)+grid(i,j+1,1:2))
3470 ps_v(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3475 p2(:) = 0.5*(grid(i,j,1:2)+grid(i+1,j,1:2))
3477 ps_u(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3488 pe_v(i,k,j) = ak(k) + ps_v(i,j)*bk(k)
3498 pe_u(i,k,j) = ak(k) + ps_u(i,j)*bk(k)
3510 p0 = (/
pi,
pi/18. /)
3517 t00 = ts0*(1.+zvir*q00)
3522 cor = 2.*
omega*sin(p0(2))
3527 p1(:) = grid(i ,j ,1:2)
3528 p2(:) = grid(i,j+1 ,1:2)
3533 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3534 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3535 d =
max(1.e-15,sqrt(d1**2+d2**2))
3540 ptmp = 0.5*(pe_v(i,k,j)+pe_v(i,k+1,j))
3541 height = (t00/gamma)*(1.-(ptmp/ps_v(i,j))**exponent)
3542 if (height > ztrop)
then 3545 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3546 - exppr*(r/rp)**exppr*
rdgas*(t00-gamma*height) &
3547 /(exppz*height*
rdgas*(t00-gamma*height)/(
grav*zp**exppz) &
3548 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3560 p1(:) = grid(i, j,1:2)
3561 p2(:) = grid(i+1,j,1:2)
3566 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3567 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3568 d =
max(1.e-15,sqrt(d1**2+d2**2))
3573 ptmp = 0.5*(pe_u(i,k,j)+pe_u(i,k+1,j))
3574 height = (t00/gamma)*(1.-(ptmp/ps_u(i,j))**exponent)
3575 if (height > ztrop)
then 3578 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3579 - exppr*(r/rp)**exppr*
rdgas*(t00-gamma*height) &
3580 /(exppz*height*
rdgas*(t00-gamma*height)/(
grav*zp**exppz) &
3581 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3593 ttrop = t00 - gamma*ztrop
3602 ptmp = 0.5*(pe(i,k,j)+pe(i,k+1,j))
3603 height = (t00/gamma)*(1.-(ptmp/ps(i,j))**exponent)
3604 if (height > ztrop)
then 3608 q(i,j,k,1) = q00*exp(-height/zq1)*exp(-(height/zq2)**exppz)
3609 p2(:) =
agrid(i,j,1:2)
3611 pt(i,j,k) = (t00-gamma*height)/(1.d0+zvir*q(i,j,k,1))/(1.d0+exppz*
rdgas*(t00-gamma*height)*height &
3612 /(
grav*zp**exppz*(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz))))
3621 ps(i,j) = pe(i,npz+1,j)
3625 if (.not.hydrostatic)
then 3629 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))
3636 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy,
ng)
3638 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3656 call dcmip16_tc (delp, pt, u, v, q, w, delz, &
3657 is, ie, js, je, isd, ied, jsd, jed, npz, ncnst, &
3658 ak, bk, ptop, pk, peln, pe, pkz, gz, phis, &
3659 ps, grid,
agrid, hydrostatic, nwat, adiabatic)
3663 call mpp_error(fatal,
" test_case not defined" )
3669 ftop =
g_sum(domain, phis(is:ie,js:je), is, ie, js, je,
ng, area, 1)
3670 if(is_master())
write(*,*)
'mean terrain height (m)=', ftop/
grav 3674 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3675 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., mountain, &
3676 moist_phys, hydrostatic, nwat, domain, .not.hydrostatic)
3679 #ifdef COLUMN_TRACER 3680 if( ncnst>1 ) q(:,:,:,2:ncnst) = 0.0
3688 p1(:) = grid(i ,j ,1:2)
3689 p2(:) = grid(i,j+1 ,1:2)
3695 if (-(r/r0)**2.0 > -40.0) q(i,j,z,1) = exp(-(r/r0)**2.0)
3731 nullify(cubed_sphere)
3736 nullify(have_south_pole)
3737 nullify(have_north_pole)
3746 subroutine get_vorticity(isc, iec, jsc, jec ,isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
3747 integer isd, ied, jsd, jed, npz
3748 integer isc, iec, jsc, jec
3749 real,
intent(in) :: u(isd:ied, jsd:jed+1, npz), v(isd:ied+1, jsd:jed, npz)
3750 real,
intent(out) :: vort(isc:iec, jsc:jec, npz)
3751 real,
intent(IN) :: dx(isd:ied,jsd:jed+1)
3752 real,
intent(IN) :: dy(isd:ied+1,jsd:jed)
3753 real,
intent(IN) :: rarea(isd:ied,jsd:jed)
3755 real :: utmp(isc:iec, jsc:jec+1), vtmp(isc:iec+1, jsc:jec)
3761 utmp(i,j) = u(i,j,k)*dx(i,j)
3766 vtmp(i,j) = v(i,j,k)*dy(i,j)
3772 vort(i,j,k) = rarea(i,j)*(utmp(i,j)-utmp(i,j+1)-vtmp(i,j)+vtmp(i+1,j))
3779 subroutine checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, &
3780 nq, km, q, lon, lat, nx, ny, rn)
3790 integer,
intent(in):: nq
3791 integer,
intent(in):: km
3792 integer,
intent(in):: i0, i1
3793 integer,
intent(in):: j0, j1
3794 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3795 real,
intent(in):: nx
3796 real,
intent(in):: ny
3797 real,
intent(in),
optional:: rn
3798 real(kind=R_GRID),
intent(in),
dimension(i0:i1,j0:j1):: lon, lat
3799 real,
intent(out):: q(ifirst:ilast,jfirst:jlast,km,nq)
3801 real:: qt(i0:i1,j0:j1)
3809 qtmp = sin(nx*lon(i,j))*sin(ny*lat(i,j))
3810 if ( qtmp < 0. )
then 3818 if (
present(rn) )
then 3826 call random_number(ftmp)
3827 q(i,j,k,iq) = qt(i,j) + rn*ftmp
3839 q(i,j,k,iq) = qt(i,j)
3849 km, q, delp, ncnst, lon, lat)
3854 integer,
intent(in):: km
3855 integer,
intent(in):: i0, i1
3856 integer,
intent(in):: j0, j1
3857 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3858 integer,
intent(in):: ncnst
3859 real(kind=R_GRID),
intent(in),
dimension(ifirst:ilast,jfirst:jlast):: lon, lat
3860 real,
intent(inout):: q(ifirst:ilast,jfirst:jlast,km,ncnst)
3861 real,
intent(in):: delp(ifirst:ilast,jfirst:jlast,km)
3863 real:: D, k1, r, ll, sinthc, costhc, mm
3869 real,
parameter :: qcly = 4.e-6
3870 real,
parameter :: lc = 5.*
pi/3.
3871 real,
parameter :: thc =
pi/9.
3872 real,
parameter :: k2 = 1.
3882 k1 =
max(0., sin(lat(i,j))*sinthc + cos(lat(i,j))*costhc*cos(lon(i,j) - lc))
3884 d = sqrt(r*r + 2.*r*qcly)
3886 q(i,j,1,cl2) = 0.5*(qcly - q(i,j,1,cl))
3893 q(i,j,k,cl) = q(i,j,1,cl)
3894 q(i,j,k,cl2) = q(i,j,1,cl2)
3901 if (is_master())
then 3906 qcly0 =
qcly0 + (q(i,j,k,cl) + 2.*q(i,j,k,cl2))*delp(i,j,k)
3907 mm = mm + delp(i,j,k)
3912 if (is_master()) print*,
' qcly0 = ',
qcly0 3921 real,
intent(in):: ubar
3922 real,
intent(in):: r0
3923 real,
intent(in):: p1(2)
3924 real,
intent(inout):: u(isd:ied, jsd:jed+1)
3925 real,
intent(inout):: v(isd:ied+1,jsd:jed)
3926 real(kind=R_GRID),
intent(IN) :: grid(isd:ied+1,jsd:jed+1,2)
3928 real(kind=R_GRID):: p2(2), p3(2), p4(2)
3929 real(kind=R_GRID):: e1(3), e2(3), ex(3), ey(3)
3930 real:: vr, r, d2, cos_p, x1, y1
3939 p2(1) = p2(1) - p1(1)
3940 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
3948 x1 = cos(p2(2))*sin(p2(1))
3949 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
3950 d2 =
max(1.e-25, sqrt(x1**2 + y1**2))
3953 p3(1) = grid(i,j, 1) - p1(1)
3954 p3(2) = grid(i,j, 2)
3955 p4(1) = grid(i+1,j,1) - p1(1)
3956 p4(2) = grid(i+1,j,2)
3968 p2(1) = p2(1) - p1(1)
3969 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
3976 x1 = cos(p2(2))*sin(p2(1))
3977 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
3978 d2 =
max(1.e-25, sqrt(x1**2 + y1**2))
3981 p3(1) = grid(i,j, 1) - p1(1)
3982 p3(2) = grid(i,j, 2)
3983 p4(1) = grid(i,j+1,1) - p1(1)
3984 p4(2) = grid(i,j+1,2)
3994 real function gh_jet(npy, lat_in)
3995 integer,
intent(in):: npy
3996 real,
intent(in):: lat_in
3997 real lat, lon, dp, uu
4004 dp =
pi /
real(jm-1)
4014 lat = -
pi/2. + (
real(j-1)-0.5)*dp
4016 ft = 2.*
omega*sin(lat)
4041 real function u_jet(lat)
4043 real umax, en, ph0, ph1
4048 en = exp( -4./(ph1-ph0)**2 )
4050 if ( lat>ph0 .and. lat<ph1 )
then 4051 u_jet = (umax/en)*exp( 1./( (lat-ph0)*(lat-ph1) ) )
4058 real,
intent(OUT) :: B(isd:ied,jsd:jed)
4059 real,
intent(IN) :: agrid(isd:ied,jsd:jed,2)
4067 if (sin(agrid(i,j,2)) > 0.)
then 4068 myc = sin(agrid(i,j,1))
4069 yy = (cos(agrid(i,j,2))/sin(agrid(i,j,2)))**2
4070 myb =
gh0*yy*exp(1.-yy)
4088 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4089 real ,
intent(IN) :: time_since_start
4095 tday = time_since_start/86400.0
4096 if (tday >= 20.)
then 4097 aoft(2) = 0.5*(1.-cos(0.25*
pi*(tday-20)))
4098 if (tday == 24)
aoft(2) = 1.0
4099 elseif (tday <= 4.)
then 4100 aoft(2) = 0.5*(1.-cos(0.25*
pi*tday))
4101 elseif (tday <= 16.)
then 4104 aoft(2) = 0.5*(1.+cos(0.25*
pi*(tday-16.)))
4109 phis(i,j) = amean*
case9_b(i,j)
4122 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4139 subroutine case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain)
4141 real,
intent(INOUT) :: delp(isd:ied,jsd:jed,npz)
4142 real,
intent(INOUT) :: uc(isd:ied+1,jsd:jed,npz)
4143 real,
intent(INOUT) :: vc(isd:ied,jsd:jed+1,npz)
4144 real,
intent(INOUT) :: u(isd:ied,jsd:jed+1,npz)
4145 real,
intent(INOUT) :: v(isd:ied+1,jsd:jed,npz)
4146 real,
intent(INOUT) :: ua(isd:ied,jsd:jed,npz)
4147 real,
intent(INOUT) :: va(isd:ied,jsd:jed,npz)
4148 real,
intent(INOUT) :: pe(is-1:ie+1, npz+1,js-1:je+1)
4149 real,
intent(IN) :: time, dt
4150 real,
intent(INOUT) :: ptop
4151 integer,
intent(IN) :: npx, npy, npz
4153 type(
domain2d),
intent(INOUT) :: domain
4160 real :: s, l, dt2, v0, phase
4161 real :: ull, vll, lonp
4162 real :: p0(2), elon(3), elat(3)
4164 real :: psi(isd:ied,jsd:jed)
4165 real :: psi_b(isd:ied+1,jsd:jed+1)
4166 real :: dist, psi1, psi2
4171 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3), pt(2), p1(2), p2(2), p3(2), rperiod, timefac, t00
4175 real(kind=R_GRID),
pointer,
dimension(:,:,:) ::
agrid, grid
4176 real,
pointer,
dimension(:,:) :: dx, dxa, dy, dya, dxc, dyc
4178 agrid => gridstruct%agrid_64
4179 grid => gridstruct%grid_64
4182 dxa => gridstruct%dxa
4183 dxc => gridstruct%dxc
4185 dya => gridstruct%dya
4186 dyc => gridstruct%dyc
4188 period =
real( 12*24*3600 ) 4193 phase =
pi*time/period
4206 omega0 = 23000.*
pi/period
4209 ptop = 100000.*exp(-12000.*
grav/t00/
rdgas)
4214 s =
min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*
pi)))
4215 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(
agrid(i,j,1)-period*(time+dt2))*cos(
agrid(i,j,2))* &
4216 cos(period*(time+dt2))*sin(s*0.5*
pi)
4224 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4241 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
4242 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
4251 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
4252 if (dist==0) vc(i,j,k) = 0.
4258 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
4259 if (dist==0) uc(i,j,k) = 0.
4266 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
4267 if (dist==0) v(i,j,k) = 0.
4273 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
4274 if (dist==0) u(i,j,k) = 0.
4280 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
4281 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
4283 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
4284 if (dist==0) ua(i,j,k) = 0.
4285 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
4286 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
4288 va(i,j,k) = (psi2 - psi1) / (dist)
4289 if (dist==0) va(i,j,k) = 0.
4295 omega0 = 23000.*
pi/period
4300 s =
min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*
pi)))
4301 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(
agrid(i,j,1)-period*(time+dt2))*cos(
agrid(i,j,2))* &
4302 cos(period*(time+dt2))*sin(s*0.5*
pi)
4310 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4320 p1(:) = grid(i ,j ,1:2)
4321 p2(:) = grid(i,j+1 ,1:2)
4325 l = p3(1) - 2.*
pi*time/period
4326 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(
pi*time/period) + 2.*
pi*
radius/period*cos(p3(2))
4327 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(
pi*time/period)
4333 p1(:) = grid(i, j,1:2)
4334 p2(:) = grid(i+1,j,1:2)
4338 l = p3(1) - 2.*
pi*time/period
4339 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(
pi*time/period) + 2.*
pi*
radius/period*cos(p3(2))
4340 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(
pi*time/period)
4363 call dtoa( u(:,:,1), v(:,:,1),ua(:,:,1),va(:,:,1),dx,dy,dxa,dya,dxc,dyc,npx,npy,
ng)
4365 call atoc(ua(:,:,1),va(:,:,1),uc(:,:,1),vc(:,:,1),dx,dy,dxa,dya,npx,npy,
ng, gridstruct%nested, domain)
4370 ua(i,j,k) = ua(i,j,1)
4375 va(i,j,k) = va(i,j,1)
4383 vc(i,j,k) = vc(i,j,1)
4388 uc(i,j,k) = uc(i,j,1)
4401 pe(i,k,j) = pe(i,k,j) + dt*omega0*
grav*pe(i,k,j)/
rdgas/300./k_cell* &
4402 (-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)))* &
4403 sin(
pi*
zz0(k)/12000.)*cos(phase)
4411 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4423 vtmp = -
radius * omega0 *
pi / k_cell / 12000. * &
4424 cos(
agrid(i,j,2)) * sin(k_cell *
agrid(i,j,2)) * &
4425 sin(
pi*
zz0(k)/12000.)*cos(phase)
4434 uc(:,:,k) = uc(:,:,1)
4435 vc(:,:,k) = vc(:,:,1)
4436 ua(:,:,k) = ua(:,:,1)
4437 va(:,:,k) = va(:,:,1)
4441 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
4458 subroutine get_stats(dt, dtout, nt, maxnt, ndays, u,v,pt,delp,q,phis, ps, &
4459 uc,vc, ua,va, npx, npy, npz, ncnst, ndims, nregions, &
4460 gridstruct, stats_lun, consv_lun, monitorFreq, tile, &
4462 integer,
intent(IN) :: nt, maxnt
4463 real ,
intent(IN) :: dt, dtout, ndays
4464 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
4465 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
4466 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
4467 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
4468 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
4469 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4470 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
4471 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
4472 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
4473 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
4474 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
4475 integer,
intent(IN) :: npx, npy, npz, ncnst, tile
4476 integer,
intent(IN) :: ndims
4477 integer,
intent(IN) :: nregions
4478 integer,
intent(IN) :: stats_lun
4479 integer,
intent(IN) :: consv_lun
4480 integer,
intent(IN) :: monitorfreq
4482 type(
domain2d),
intent(INOUT) :: domain
4483 logical,
intent(IN) :: nested
4488 real :: pmin, pmin1, uamin1, vamin1
4489 real :: pmax, pmax1, uamax1, vamax1
4490 real(kind=4) :: arr_r4(5)
4491 real :: tmass0, tvort0, tener0, tke0
4492 real :: tmass, tvort, tener, tke
4493 real :: temp(is:ie,js:je)
4494 integer :: i0, j0, k0, n0
4495 integer :: i, j, k, n, iq
4497 real :: psmo, vtx, p, w_p, p0
4498 real :: x1,y1,z1,x2,y2,z2,ang
4500 real :: p1(2), p2(2), p3(2), r, r0, dist, heading
4502 real :: uc0(isd:ied+1,jsd:jed ,npz)
4503 real :: vc0(isd:ied ,jsd:jed+1,npz)
4508 real,
save,
allocatable,
dimension(:,:,:) :: u0, v0
4509 real ::
up(isd:ied ,jsd:jed+1,npz)
4510 real :: vp(isd:ied+1,jsd:jed ,npz)
4512 real,
dimension(:,:,:),
pointer :: grid,
agrid 4513 real,
dimension(:,:),
pointer :: area,
f0, dx, dy, dxa, dya, dxc, dyc
4515 grid => gridstruct%grid
4516 agrid=> gridstruct%agrid
4518 area => gridstruct%area
4523 dxa => gridstruct%dxa
4524 dya => gridstruct%dya
4525 dxc => gridstruct%dxc
4526 dyc => gridstruct%dyc
4529 if (nt == 0 .and. is_master()) print*,
'INITIALIZING GET_STATS' 4532 myday = ndays*((float(nt)/float(maxnt)))
4534 #if defined(SW_DYNAMICS) 4543 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
4545 if (p /= 0.0) w_p = vtx/p
4547 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
4557 dist = 2.0*
pi*
radius* ((float(nt)/float(maxnt)))
4563 p2(1) =
agrid(i,j,1)
4564 p2(2) =
agrid(i,j,2)
4567 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(
pi*r/r0))
4569 phi0(i,j,1) = phis(i,j)
4576 call pmxn(delp(:,:,1), npx, npy, nregions, tile, gridstruct, pmin1, pmax1, i0, j0, n0)
4580 call get_scalar_stats( delp(:,:,1),
phi0(:,:,1), npx, npy, ndims, nregions, &
4581 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4588 arr_r4(5) = linf_norm
4601 200
format(i6.6,a,i6.6,a,e21.14)
4602 201
format(
' ',a,e21.14,
' ',e21.14)
4603 202
format(
' ',a,i4.4,
'x',i4.4,
'x',i4.4)
4605 if ( (is_master()) .and. mod(nt,monitorfreq)==0 )
then 4606 write(*,200) nt,
' step of ', maxnt,
' DAY ', myday
4607 write(*,201)
'Height MAX : ', pmax1
4608 write(*,201)
'Height MIN : ', pmin1
4609 write(*,202)
'HGT MAX location : ', i0, j0, n0
4611 write(*,201)
'Height L1_norm : ', l1_norm
4612 write(*,201)
'Height L2_norm : ', l2_norm
4613 write(*,201)
'Height Linf_norm : ', linf_norm
4618 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy,
ng)
4619 call pmxn(ua(:,:,1), npx, npy, nregions, tile, gridstruct, pmin1, pmax1, i0, j0, n0)
4621 call get_vector_stats( ua(:,:,1),
ua0(:,:,1), va(:,:,1),
va0(:,:,1), npx, npy, ndims, nregions, &
4622 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4628 arr_r4(5) = linf_norm
4630 if ( (is_master()) .and. mod(nt,monitorfreq)==0)
then 4631 write(*,201)
'UV MAX : ', pmax1
4632 write(*,201)
'UV MIN : ', pmin1
4633 write(*,202)
'UV MAX location : ', i0, j0, n0
4635 write(*,201)
'UV L1_norm : ', l1_norm
4636 write(*,201)
'UV L2_norm : ', l2_norm
4637 write(*,201)
'UV Linf_norm : ', linf_norm
4642 200
format(i6.6,a,i6.6,a,e10.4)
4643 201
format(
' ',a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4)
4644 202
format(
' ',a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4,
' ',e10.4)
4645 203
format(
' ',a,i3.3,a,e10.4,
' ',e10.4,
' ',i4.4,
'x',i4.4,
'x',i4.4,
'x',i4.4)
4647 if(is_master())
write(*,200) nt,
' step of ', maxnt,
' DAY ', myday
4650 psmo =
globalsum(ps(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4651 if(is_master())
write(*,*)
' Total surface pressure =', 0.01*psmo
4652 call pmxn(ps, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4653 if (is_master())
then 4654 write(*,201)
'PS MAX|MIN : ', 0.01*pmax, 0.01*pmin, i0, j0, n0
4665 call pmxn(pt(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4666 pmin1 =
min(pmin, pmin1)
4667 pmax1 =
max(pmax, pmax1)
4668 if (pmax1 == pmax) k0 = k
4670 if (is_master())
then 4671 write(*,201)
'PT MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4674 #if defined(DEBUG_TEST_CASES) 4675 if(is_master())
write(*,*)
' ' 4683 call pmxn(pt(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4684 pmin1 =
min(pmin, pmin1)
4685 pmax1 =
max(pmax, pmax1)
4686 if (is_master())
then 4687 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) )
4690 if(is_master())
write(*,*)
' ' 4701 call pmxn(delp(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4702 pmin1 =
min(pmin, pmin1)
4703 pmax1 =
max(pmax, pmax1)
4704 if (pmax1 == pmax) k0 = k
4706 if (is_master())
then 4707 write(*,201)
'Delp MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4718 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)
4719 call pmxn(ua(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4720 uamin1 =
min(pmin, uamin1)
4721 uamax1 =
max(pmax, uamax1)
4722 if (uamax1 == pmax) k0 = k
4724 if (is_master())
then 4725 write(*,201)
'U MAX|MIN : ', uamax1, uamin1, i0, j0, k0, n0
4735 call pmxn(va(:,:,k), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4736 vamin1 =
min(pmin, vamin1)
4737 vamax1 =
max(pmax, vamax1)
4738 if (vamax1 == pmax) k0 = k
4740 if (is_master())
then 4741 write(*,201)
'V MAX|MIN : ', vamax1, vamin1, i0, j0, k0, n0
4752 call pmxn(q(isd,jsd,k,1), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4753 pmin1 =
min(pmin, pmin1)
4754 pmax1 =
max(pmax, pmax1)
4755 if (pmax1 == pmax) k0 = k
4757 if (is_master())
then 4758 write(*,201)
'Q MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4770 call pmxn(q(isd,jsd,k,iq), npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
4771 pmin1 =
min(pmin, pmin1)
4772 pmax1 =
max(pmax, pmax1)
4773 if (pmax1 == pmax) k0 = k
4775 if (is_master())
then 4776 write(*,203)
'TR',iq-1,
' MAX|MIN : ', pmax1, pmin1, i0, j0, k0, n0
4784 call get_vector_stats( ua(:,:,22),
ua0(:,:,22), va(:,:,22),
va0(:,:,22), npx, npy, ndims, nregions, &
4785 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
4786 if (is_master())
then 4787 write(*,201)
'UV(850) L1_norm : ', l1_norm
4788 write(*,201)
'UV(850) L2_norm : ', l2_norm
4789 write(*,201)
'UV(850) Linf_norm : ', linf_norm
4797 #if defined(SW_DYNAMICS) 4805 temp(:,:) = delp(is:ie,js:je,k)
4806 tmass0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4807 tmass = tmass + tmass0
4811 allocate(u0(isd:ied,jsd:jed+1,npz))
4812 allocate(v0(isd:ied+1,jsd:jed,npz))
4821 call dtoa(
up(isd,jsd,k), vp(isd,jsd,k), ua, va, dx,dy, dxa, dya, dxc, dyc, npx, npy,
ng)
4822 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.)
4826 temp(i,j) = ( uc0(i,j,k)*uc0(i,j,k) + uc0(i+1,j,k)*uc0(i+1,j,k) + &
4827 vc0(i,j,k)*vc0(i,j,k) + vc0(i,j+1,k)*vc0(i,j+1,k) )
4830 tke0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4836 temp(i,j) = 0.5 * (delp(i,j,k)/
grav) * temp(i,j)
4837 temp(i,j) = temp(i,j) + &
4838 grav*((delp(i,j,k)/
grav + phis(i,j))*(delp(i,j,k)/
grav + phis(i,j))) - &
4842 tener0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4843 tener = tener + tener0
4849 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)) - &
4850 (u(i,j+1,k)*dx(i,j+1) - u(i,j,k)*dx(i,j)) )
4851 temp(i,j) = (
grav*(temp(i,j)*temp(i,j))/delp(i,j,k) )
4854 tvort0 =
globalsum(temp, npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4855 tvort = tvort + tvort0
4871 #if defined(SW_DYNAMICS) 4874 myrec = myday*86400.0/dtout + 1
4876 if (is_master())
write(consv_lun,rec=myrec) arr_r4(1:4)
4877 #if defined(SW_DYNAMICS) 4878 if ( (is_master()) .and. mod(nt,monitorfreq)==0)
then 4880 if ( (is_master()) )
then 4882 write(*,201)
'MASS TOTAL : ', tmass
4885 write(*,201)
'Kinetic Energy KE : ', tke
4886 write(*,201)
'ENERGY TOTAL : ', tener
4888 write(*,201)
'ENSTR TOTAL : ', tvort
4913 real ,
intent(IN) :: p1(2), p2(2)
4914 real ,
intent(IN) :: dist
4915 real ,
intent(IN) :: heading
4916 real ,
intent(OUT) :: p3(2)
4922 p3(2) = asin( (cos(heading)*cos(p1(2))*sin(pha)) + (sin(p1(2))*cos(pha)) )
4923 dp = atan2( sin(heading)*sin(pha)*cos(p1(2)) , cos(pha) - sin(p1(2))*sin(p3(2)) )
4924 p3(1) = mod( (p1(1)-
pi)-dp+
pi , 2.*
pi )
4941 vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
4942 integer,
intent(IN) :: npx, npy
4943 integer,
intent(IN) :: ndims
4944 integer,
intent(IN) :: nregions, tile
4945 real ,
intent(IN) :: var(isd:ied,jsd:jed)
4946 real ,
intent(IN) :: varT(isd:ied,jsd:jed)
4947 real ,
intent(OUT) :: vmin
4948 real ,
intent(OUT) :: vmax
4949 real ,
intent(OUT) :: L1_norm
4950 real ,
intent(OUT) :: L2_norm
4951 real ,
intent(OUT) :: Linf_norm
4953 type(fv_grid_type),
target :: gridstruct
4962 real :: varSUM, varSUM2, varMAX
4964 real :: vminT, vmaxT, vmeanT, vvarT
4965 integer :: i0, j0, n0
4967 real,
dimension(:,:,:),
pointer :: grid, agrid
4968 real,
dimension(:,:),
pointer :: area
4970 grid => gridstruct%grid
4971 agrid=> gridstruct%agrid
4973 area => gridstruct%area
4992 vmean =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4993 vmeant =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
4994 vmean = vmean / (4.0*
pi)
4995 vmeant = vmeant / (4.0*
pi)
4997 call pmxn(var, npx, npy, nregions, tile, gridstruct, vmin , vmax , i0, j0, n0)
4998 call pmxn(vart, npx, npy, nregions, tile, gridstruct, vmint, vmaxt, i0, j0, n0)
4999 call pmxn(var-vart, npx, npy, nregions, tile, gridstruct, pdiffmn, pdiffmx, i0, j0, n0)
5001 vmax = (vmax - vmaxt) / (vmaxt-vmint)
5002 vmin = (vmin - vmint) / (vmaxt-vmint)
5004 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5005 varsum2 =
globalsum(vart(is:ie,js:je)**2., npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5006 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)
5007 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)
5008 l1_norm = l1_norm/varsum
5009 l2_norm = sqrt(l2_norm)/sqrt(varsum2)
5011 call pmxn(abs(vart), npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5013 call pmxn(abs(var-vart), npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5014 linf_norm = vmax/varmax
5029 npx, npy, ndims, nregions, &
5030 vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
5031 integer,
intent(IN) :: npx, npy
5032 integer,
intent(IN) :: ndims
5033 integer,
intent(IN) :: nregions, tile
5034 real ,
intent(IN) :: varU(isd:ied,jsd:jed)
5035 real ,
intent(IN) :: varUT(isd:ied,jsd:jed)
5036 real ,
intent(IN) :: varV(isd:ied,jsd:jed)
5037 real ,
intent(IN) :: varVT(isd:ied,jsd:jed)
5038 real ,
intent(OUT) :: vmin
5039 real ,
intent(OUT) :: vmax
5040 real ,
intent(OUT) :: L1_norm
5041 real ,
intent(OUT) :: L2_norm
5042 real ,
intent(OUT) :: Linf_norm
5044 real :: var(isd:ied,jsd:jed)
5045 real :: varT(isd:ied,jsd:jed)
5053 real :: varSUM, varSUM2, varMAX
5055 real :: vminT, vmaxT, vmeanT, vvarT
5057 integer :: i0, j0, n0
5059 type(fv_grid_type),
target :: gridstruct
5061 real,
dimension(:,:,:),
pointer :: grid, agrid
5062 real,
dimension(:,:),
pointer :: area
5064 grid => gridstruct%grid
5065 agrid=> gridstruct%agrid
5067 area => gridstruct%area
5088 var(i,j) = sqrt( (varu(i,j)-varut(i,j))**2. + &
5089 (varv(i,j)-varvt(i,j))**2. )
5090 vart(i,j) = sqrt( varut(i,j)*varut(i,j) + &
5091 varvt(i,j)*varvt(i,j) )
5094 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5095 l1_norm =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5096 l1_norm = l1_norm/varsum
5098 call pmxn(vart, npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5100 call pmxn(var, npx, npy, nregions, tile, gridstruct, vmin, vmax, i0, j0, n0)
5101 linf_norm = vmax/varmax
5105 var(i,j) = ( (varu(i,j)-varut(i,j))**2. + &
5106 (varv(i,j)-varvt(i,j))**2. )
5107 vart(i,j) = ( varut(i,j)*varut(i,j) + &
5108 varvt(i,j)*varvt(i,j) )
5111 varsum =
globalsum(vart(is:ie,js:je), npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5112 l2_norm =
globalsum(var(is:ie,js:je) , npx, npy, is,ie, js,je, isd, ied, jsd, jed, gridstruct, tile)
5113 l2_norm = sqrt(l2_norm)/sqrt(varsum)
5127 real,
intent(IN) :: ndt
5128 integer,
intent(IN) :: n_split
5129 integer,
intent(IN) :: npx, npy, npz, tile
5130 logical,
OPTIONAL,
intent(IN) :: noprint
5131 real ,
intent(IN) :: uc(isd:ied+1,jsd:jed ,npz)
5132 real ,
intent(IN) :: vc(isd:ied ,jsd:jed+1,npz)
5134 real :: ideal_c=0.06
5135 real :: tolerance= 1.e-3
5136 real :: dt_inc, dt_orig
5137 real :: meancy, mincy, maxcy, meancx, mincx, maxcx
5146 real,
dimension(:,:),
pointer :: dxc, dyc
5148 dxc => gridstruct%dxc
5149 dyc => gridstruct%dyc
5151 dt = ndt/
real(n_split)
5153 300
format(i4.4,
' ',i4.4,
' ',i4.4,
' ',i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
5159 do while(.not. ideal)
5171 mincx =
min(mincx, abs( (dt/dxc(i,j))*uc(i,j,k) ))
5172 maxcx =
max(maxcx, abs( (dt/dxc(i,j))*uc(i,j,k) ))
5173 meancx = meancx + abs( (dt/dxc(i,j))*uc(i,j,k) )
5175 if (abs( (dt/dxc(i,j))*uc(i,j,k) ) > 1.0)
then 5177 write(*,300) i,j,k,tile, abs( (dt/dxc(i,j))*uc(i,j,k) ), dt, dxc(i,j), uc(i,j,k), counter
5185 mincy =
min(mincy, abs( (dt/dyc(i,j))*vc(i,j,k) ))
5186 maxcy =
max(maxcy, abs( (dt/dyc(i,j))*vc(i,j,k) ))
5187 meancy = meancy + abs( (dt/dyc(i,j))*vc(i,j,k) )
5189 if (abs( (dt/dyc(i,j))*vc(i,j,k) ) > 1.0)
then 5191 write(*,300) i,j,k,tile, abs( (dt/dyc(i,j))*vc(i,j,k) ), dt, dyc(i,j), vc(i,j,k), counter
5199 call mp_reduce_max(maxcx)
5200 call mp_reduce_max(maxcy)
5203 call mp_reduce_max(mincx)
5204 call mp_reduce_max(mincy)
5207 call mp_reduce_sum(meancx)
5208 call mp_reduce_sum(meancy)
5209 meancx = meancx/(6.0*dble(npx)*dble(npy-1))
5210 meancy = meancy/(6.0*dble(npx-1)*dble(npy))
5222 if ( (.not.
present(noprint)) .and. (is_master()) )
then 5224 print*,
'--------------------------------------------' 5225 print*,
'Y-dir Courant number MIN : ', mincy
5226 print*,
'Y-dir Courant number MAX : ', maxcy
5228 print*,
'X-dir Courant number MIN : ', mincx
5229 print*,
'X-dir Courant number MAX : ', maxcx
5231 print*,
'X-dir Courant number MEAN : ', meancx
5232 print*,
'Y-dir Courant number MEAN : ', meancy
5234 print*,
'NDT: ', ndt
5235 print*,
'n_split: ', n_split
5238 print*,
'--------------------------------------------' 5252 subroutine pmxn(p, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
5253 integer,
intent(IN) :: npx
5254 integer,
intent(IN) :: npy
5255 integer,
intent(IN) :: nregions, tile
5256 real ,
intent(IN) :: p(isd:ied,jsd:jed)
5257 type(fv_grid_type),
intent(IN),
target :: gridstruct
5258 real ,
intent(OUT) :: pmin
5259 real ,
intent(OUT) :: pmax
5260 integer,
intent(OUT) :: i0
5261 integer,
intent(OUT) :: j0
5262 integer,
intent(OUT) :: n0
5268 real,
pointer,
dimension(:,:,:) :: agrid, grid
5269 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
5270 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
5271 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
5272 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5274 logical,
pointer :: cubed_sphere, latlon
5276 logical,
pointer :: have_south_pole, have_north_pole
5278 integer,
pointer :: ntiles_g
5279 real,
pointer :: acapN, acapS, globalarea
5281 grid => gridstruct%grid
5282 agrid=> gridstruct%agrid
5284 area => gridstruct%area
5285 rarea => gridstruct%rarea
5290 ee1 => gridstruct%ee1
5291 ee2 => gridstruct%ee2
5294 en1 => gridstruct%en1
5295 en2 => gridstruct%en2
5299 dxa => gridstruct%dxa
5300 dya => gridstruct%dya
5301 rdxa => gridstruct%rdxa
5302 rdya => gridstruct%rdya
5303 dxc => gridstruct%dxc
5304 dyc => gridstruct%dyc
5306 cubed_sphere => gridstruct%cubed_sphere
5307 latlon => gridstruct%latlon
5309 have_south_pole => gridstruct%have_south_pole
5310 have_north_pole => gridstruct%have_north_pole
5312 ntiles_g => gridstruct%ntiles_g
5313 acapn => gridstruct%acapN
5314 acaps => gridstruct%acapS
5315 globalarea => gridstruct%globalarea
5326 if (temp > pmax)
then 5330 elseif (temp < pmin)
then 5337 call mp_reduce_max(temp)
5338 if (temp /= pmax)
then 5344 call mp_reduce_max(i0)
5345 call mp_reduce_max(j0)
5346 call mp_reduce_max(n0)
5349 call mp_reduce_max(pmin)
5365 subroutine output_ncdf(dt, nt, maxnt, nout, u,v,pt,delp,q,phis,ps, uc,vc, ua,va, &
5366 omga, npx, npy, npz, ng, ncnst, ndims, nregions, ncid, &
5367 npx_p1_id, npy_p1_id, npx_id, npy_id, npz_id, ntiles_id, ncnst_id, nt_id, &
5368 phis_id, delp_id, ps_id, pt_id, pv_id, om_id, u_id, v_id, q_id, tracers_ids, &
5369 lats_id, lons_id, gridstruct, flagstruct)
5370 real,
intent(IN) :: dt
5371 integer,
intent(IN) :: nt, maxnt
5372 integer,
intent(INOUT) :: nout
5374 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
5375 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
5376 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
5377 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
5378 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
5380 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
5381 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
5383 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
5384 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
5385 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
5386 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
5387 real ,
intent(INOUT) :: omga(isd:ied ,jsd:jed ,npz)
5389 integer,
intent(IN) :: npx, npy, npz
5390 integer,
intent(IN) :: ng, ncnst
5391 integer,
intent(IN) :: ndims
5392 integer,
intent(IN) :: nregions
5393 integer,
intent(IN) :: ncid
5394 integer,
intent(IN) :: npx_p1_id, npy_p1_id, npx_id, npy_id, npz_id, ncnst_id
5395 integer,
intent(IN) :: ntiles_id, nt_id
5396 integer,
intent(IN) :: phis_id, delp_id, ps_id, pt_id, pv_id, u_id, v_id, q_id
5397 integer,
intent(IN) :: om_id
5398 integer,
intent(IN) :: tracers_ids(ncnst-1)
5399 integer,
intent(IN) :: lats_id, lons_id
5401 type(fv_grid_type),
target :: gridstruct
5402 type(fv_flags_type),
intent(IN) :: flagstruct
5404 real,
allocatable :: tmp(:,:,:)
5405 real,
allocatable :: tmpA(:,:,:)
5406 #if defined(SW_DYNAMICS) 5407 real,
allocatable :: ut(:,:,:)
5408 real,
allocatable :: vt(:,:,:)
5410 real,
allocatable :: ut(:,:,:,:)
5411 real,
allocatable :: vt(:,:,:,:)
5412 real,
allocatable :: tmpA_3d(:,:,:,:)
5414 real,
allocatable :: vort(:,:)
5421 real :: utmp, vtmp, r, r0, dist, heading
5422 integer :: i,j,k,n,iq,nreg
5425 real :: x1,y1,z1,x2,y2,z2,ang
5427 real,
pointer,
dimension(:,:,:) :: agrid, grid
5428 real,
pointer,
dimension(:,:) :: area, rarea
5429 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5431 grid => gridstruct%grid
5432 agrid => gridstruct%agrid
5434 area => gridstruct%area
5435 rarea => gridstruct%rarea
5439 dxa => gridstruct%dxa
5440 dya => gridstruct%dya
5441 rdxa => gridstruct%rdxa
5442 rdya => gridstruct%rdya
5443 dxc => gridstruct%dxc
5444 dyc => gridstruct%dyc
5446 allocate( tmp(npx ,npy ,nregions) )
5447 allocate( tmpa(npx-1,npy-1,nregions) )
5448 #if defined(SW_DYNAMICS) 5449 allocate( ut(npx-1,npy-1,nregions) )
5450 allocate( vt(npx-1,npy-1,nregions) )
5452 allocate( ut(npx-1,npy-1,npz,nregions) )
5453 allocate( vt(npx-1,npy-1,npz,nregions) )
5454 allocate( tmpa_3d(npx-1,npy-1,npz,nregions) )
5456 allocate( vort(isd:ied,jsd:jed) )
5461 tmp(is:ie+1,js:je+1,tile) = grid(is:ie+1,js:je+1,2)
5462 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)
5463 tmp(is:ie+1,js:je+1,tile) = grid(is:ie+1,js:je+1,1)
5464 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)
5467 #if defined(SW_DYNAMICS) 5469 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)/
grav 5478 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
5479 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0) /
grav 5493 dist = 2.0*
pi*
radius* ((float(nt)/float(maxnt)))
5498 p2(1) = agrid(i,j,1)
5499 p2(2) = agrid(i,j,2)
5502 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(
pi*r/r0))
5504 phi0(i,j,1) = phis(i,j)
5516 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
5518 if (p /= 0.0) w_p = vtx/p
5519 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
5524 tmpa(is:ie,js:je,tile) =
phi0(is:ie,js:je,1)
5525 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5526 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)
5528 call wrtvar_ncdf(ncid, ps_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5534 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)) - &
5535 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
5536 vort(i,j) =
grav*vort(i,j)/delp(i,j,1)
5539 tmpa(is:ie,js:je,tile) = vort(is:ie,js:je)
5540 call wrtvar_ncdf(ncid, pv_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5543 call cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, 1, 1, gridstruct%grid_type, gridstruct%nested, flagstruct%c2l_ord, bd)
5546 ut(i,j,tile) = ua(i,j,1)
5547 vt(i,j,tile) = va(i,j,1)
5551 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)
5552 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)
5554 if ((
test_case >= 2) .and. (nt==0) )
then 5555 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/
grav 5556 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa, 3)
5561 tmpa_3d(is:ie,js:je,1:npz,tile) = q(is:ie,js:je,1:npz,1)
5562 call wrtvar_ncdf(ncid, q_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5566 tmpa_3d(is:ie,js:je,1:npz,tile) = q(is:ie,js:je,1:npz,iq)
5567 call wrtvar_ncdf(ncid, tracers_ids(iq-1), nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5571 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/
grav 5572 call wrtvar_ncdf(ncid, phis_id, nout, is,ie, js,je, npx, npy, 1, nregions, tmpa, 3)
5575 tmpa(is:ie,js:je,tile) = ps(is:ie,js:je)
5576 call wrtvar_ncdf(ncid, ps_id, nout, is,ie, js,je, npx, npy, 1, nregions, tmpa, 3)
5578 tmpa_3d(is:ie,js:je,k,tile) = delp(is:ie,js:je,k)/
grav 5580 call wrtvar_ncdf(ncid, delp_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5584 tmpa_3d(is:ie,js:je,k,tile) = pt(is:ie,js:je,k)
5586 call wrtvar_ncdf(ncid, pt_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5589 call cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, npz, gridstruct%grid_type, gridstruct%nested, flagstruct%c2l_ord)
5593 ut(i,j,k,tile) = ua(i,j,k)
5594 vt(i,j,k,tile) = va(i,j,k)
5598 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)
5599 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)
5606 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)) - &
5607 (u(i,j+1,k)*dx(i,j+1) - u(i,j,k)*dx(i,j)) )
5611 call wrtvar_ncdf(ncid, pv_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5617 tmpa_3d(i,j,k,tile) = omga(i,j,k)
5621 call wrtvar_ncdf(ncid, om_id, nout, is,ie, js,je, npx, npy, npz, nregions, tmpa_3d, 4)
5627 #if defined(SW_DYNAMICS) 5633 deallocate( tmpa_3d )
5651 end subroutine output_ncdf
5662 subroutine output(dt, nt, maxnt, nout, u,v,pt,delp,q,phis,ps, uc,vc, ua,va, &
5663 npx, npy, npz, ng, ncnst, ndims, nregions, phis_lun, phi_lun, &
5664 pt_lun, pv_lun, uv_lun, gridstruct)
5666 real,
intent(IN) :: dt
5667 integer,
intent(IN) :: nt, maxnt
5668 integer,
intent(INOUT) :: nout
5670 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
5671 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
5672 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
5673 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
5674 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
5676 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
5677 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
5679 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
5680 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
5681 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
5682 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
5684 integer,
intent(IN) :: npx, npy, npz
5685 integer,
intent(IN) :: ng, ncnst
5686 integer,
intent(IN) :: ndims
5687 integer,
intent(IN) :: nregions
5688 integer,
intent(IN) :: phis_lun, phi_lun, pt_lun, pv_lun, uv_lun
5690 type(fv_grid_type),
target :: gridstruct
5692 real :: tmp(1-ng:npx +ng,1-ng:npy +ng,1:nregions)
5693 real :: tmpA(1-ng:npx-1+ng,1-ng:npy-1+ng,1:nregions)
5699 real :: ut(1:npx,1:npy,1:nregions)
5700 real :: vt(1:npx,1:npy,1:nregions)
5701 real :: utmp, vtmp, r, r0, dist, heading
5702 integer :: i,j,k,n,nreg
5703 real :: vort(isd:ied,jsd:jed)
5706 real :: x1,y1,z1,x2,y2,z2,ang
5708 real,
pointer,
dimension(:,:,:) :: agrid, grid
5709 real,
pointer,
dimension(:,:) :: area, rarea
5710 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
5712 grid => gridstruct%grid
5713 agrid => gridstruct%agrid
5715 area => gridstruct%area
5719 dxa => gridstruct%dxa
5720 dya => gridstruct%dya
5721 rdxa => gridstruct%rdxa
5722 rdya => gridstruct%rdya
5723 dxc => gridstruct%dxc
5724 dyc => gridstruct%dyc
5726 cubed_sphere => gridstruct%cubed_sphere
5730 #if defined(SW_DYNAMICS) 5732 call atob_s(delp(:,:,1)/
grav, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5733 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)/
grav 5742 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
5743 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0) /
grav 5757 dist = 2.0*
pi*
radius* ((float(nt)/float(maxnt)))
5762 p2(1) = agrid(i,j,1)
5763 p2(2) = agrid(i,j,2)
5766 phi0(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(
pi*r/r0))
5768 phi0(i,j,1) = phis(i,j)
5780 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
5782 if (p /= 0.0) w_p = vtx/p
5783 phi0(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*(nt*dt/86400.0)) )
5788 call atob_s(
phi0(:,:,1), tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5789 tmpa(is:ie,js:je,tile) =
phi0(is:ie,js:je,1)
5790 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5791 call atob_s(delp(:,:,1), tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5792 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,1)
5795 call wrt2d(phi_lun, nout, is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5801 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)) - &
5802 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
5803 vort(i,j) =
grav*vort(i,j)/delp(i,j,1)
5806 call atob_s(vort, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5807 call wrt2d(pv_lun, nout, is,ie+1, js,je+1, npx+1, npy+1, nregions, tmp(1:npx,1:npy,1:nregions))
5810 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy, ng)
5812 if (cubed_sphere)
then 5821 if (cubed_sphere)
call rotate_winds(utmp,vtmp, p1,p2,p3,p4, agrid(i,j,1:2), 2, 2)
5828 call wrt2d(uv_lun, 2*(nout-1) + 1, is,ie, js,je, npx, npy, nregions, ut(1:npx-1,1:npy-1,1:nregions))
5829 call wrt2d(uv_lun, 2*(nout-1) + 2, is,ie, js,je, npx, npy, nregions, vt(1:npx-1,1:npy-1,1:nregions))
5831 if ((
test_case >= 2) .and. (nt==0) )
then 5832 call atob_s(phis/
grav, tmp(isd:ied+1,jsd:jed+1,tile), npx,npy, dxa, dya, gridstruct%nested)
5834 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/
grav 5835 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5841 tmpa(is:ie,js:je,tile) = phis(is:ie,js:je)/
grav 5842 call wrt2d(phis_lun, nout , is,ie, js,je, npx, npy, nregions, tmpa(1:npx-1,1:npy-1,1:nregions))
5852 tmpa(is:ie,js:je,tile) = ps(is:ie,js:je)
5853 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))
5855 tmpa(is:ie,js:je,tile) = delp(is:ie,js:je,k)/
grav 5856 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))
5861 tmpa(is:ie,js:je,tile) = pt(is:ie,js:je,k)
5862 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))
5867 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)
5869 if (cubed_sphere)
then 5878 if (cubed_sphere)
call rotate_winds(utmp,vtmp, p1,p2,p3,p4, agrid(i,j,1:2), 2, 2)
5884 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))
5885 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))
5903 nullify(cubed_sphere)
5905 end subroutine output
5914 subroutine wrtvar_ncdf(ncid, varid, nrec, i1,i2, j1,j2, npx, npy, npz, ntiles, p, ndims)
5915 #include <netcdf.inc> 5916 integer,
intent(IN) :: ncid, varid
5917 integer,
intent(IN) :: nrec
5918 integer,
intent(IN) :: i1,i2,j1,j2
5919 integer,
intent(IN) :: npx
5920 integer,
intent(IN) :: npy
5921 integer,
intent(IN) :: npz
5922 integer,
intent(IN) :: ntiles
5923 real ,
intent(IN) :: p(npx-1,npy-1,npz,ntiles)
5924 integer,
intent(IN) :: ndims
5927 real(kind=4),
allocatable :: p_R4(:,:,:,:)
5929 integer :: istart(ndims+1), icount(ndims+1)
5931 allocate( p_r4(npx-1,npy-1,npz,ntiles) )
5934 p_r4(i1:i2,j1:j2,1:npz,tile) = p(i1:i2,j1:j2,1:npz,tile)
5935 call mp_gather(p_r4, i1,i2, j1,j2, npx-1, npy-1, npz, ntiles)
5938 istart(ndims+1) = nrec
5942 if (ndims == 3) icount(3) = ntiles
5943 if (ndims == 4) icount(4) = ntiles
5946 if (is_master())
then 5947 error = nf_put_vara_real(ncid, varid, istart, icount, p_r4)
5952 end subroutine wrtvar_ncdf
5961 subroutine wrt2d(iout, nrec, i1,i2, j1,j2, npx, npy, nregions, p)
5962 integer,
intent(IN) :: iout
5963 integer,
intent(IN) :: nrec
5964 integer,
intent(IN) :: i1,i2,j1,j2
5965 integer,
intent(IN) :: npx
5966 integer,
intent(IN) :: npy
5967 integer,
intent(IN) :: nregions
5968 real ,
intent(IN) :: p(npx-1,npy-1,nregions)
5970 real(kind=4) :: p_R4(npx-1,npy-1,nregions)
5976 p_r4(i,j,n) = p(i,j,n)
5981 call mp_gather(p_r4, i1,i2, j1,j2, npx-1, npy-1, nregions)
5983 if (is_master())
then 5984 write(iout,rec=nrec) p_r4(1:npx-1,1:npy-1,1:nregions)
5987 end subroutine wrt2d
5996 subroutine init_double_periodic(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
5997 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, &
5998 mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd)
6002 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6003 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6004 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
6005 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6006 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6007 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
6009 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
6011 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
6012 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
6013 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
6014 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
6015 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
6016 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6017 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6018 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6019 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6020 real ,
intent(inout) :: delz(bd%isd:,bd%jsd:,1:)
6021 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
6023 real ,
intent(inout) :: ak(npz+1)
6024 real ,
intent(inout) :: bk(npz+1)
6026 integer,
intent(IN) :: npx, npy, npz
6027 integer,
intent(IN) ::
ng, ncnst, nwat
6028 integer,
intent(IN) :: ndims
6029 integer,
intent(IN) :: nregions
6031 real,
intent(IN) :: dry_mass
6032 logical,
intent(IN) :: mountain
6033 logical,
intent(IN) :: moist_phys
6034 logical,
intent(IN) :: hydrostatic, hybrid_z
6035 integer,
intent(INOUT) :: ks
6036 integer,
intent(INOUT),
target :: tile_in
6037 real,
intent(INOUT) :: ptop
6039 type(
domain2d),
intent(IN),
target :: domain_in
6044 real,
dimension(bd%is:bd%ie):: pm, qs
6045 real,
dimension(1:npz):: pk1, ts1, qs1
6047 real :: dist, r0, f0_const, prf, rgrav
6048 real :: ptmp, ze, zc, zm, utmp, vtmp
6049 real :: t00, p00, xmax, xc, xx, yy, pk0, pturb, ztop
6053 integer :: i, j, k, m, icenter, jcenter
6055 real,
pointer,
dimension(:,:,:) ::
agrid, grid
6056 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
6057 real,
pointer,
dimension(:,:) :: rarea, fc,
f0 6058 real,
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
6059 real,
pointer,
dimension(:,:,:,:) :: ew, es
6060 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
6062 logical,
pointer :: cubed_sphere, latlon
6065 integer,
pointer :: tile
6067 logical,
pointer :: have_south_pole, have_north_pole
6069 integer,
pointer :: ntiles_g
6070 real,
pointer :: acapn, acaps, globalarea
6072 real(kind=R_GRID),
pointer :: dx_const, dy_const
6074 integer :: is, ie, js, je
6075 integer :: isd, ied, jsd, jed
6086 agrid => gridstruct%agrid
6087 grid => gridstruct%grid
6089 area => gridstruct%area_64
6093 dxa => gridstruct%dxa
6094 dya => gridstruct%dya
6095 rdxa => gridstruct%rdxa
6096 rdya => gridstruct%rdya
6097 dxc => gridstruct%dxc
6098 dyc => gridstruct%dyc
6104 dx_const => flagstruct%dx_const
6105 dy_const => flagstruct%dy_const
6110 have_south_pole => gridstruct%have_south_pole
6111 have_north_pole => gridstruct%have_north_pole
6113 ntiles_g => gridstruct%ntiles_g
6114 acapn => gridstruct%acapN
6115 acaps => gridstruct%acapS
6116 globalarea => gridstruct%globalarea
6118 f0_const = 2.*
omega*sin(flagstruct%deglat/180.*
pi)
6139 if (j>0 .and. j<5)
then 6141 if (i>0 .and. i<5)
then 6154 r0 = 5.*sqrt(dx_const**2 + dy_const**2)
6159 dist=(i-icenter)*dx_const*(i-icenter)*dx_const &
6160 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6161 dist=
min(r0,sqrt(dist))
6162 phis(i,j)=1500.*(1. - (dist/r0))
6183 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
6184 delp, ak, bk, pt, delz, area,
ng, .false., hydrostatic, hybrid_z, domain)
6188 r0 = 100.*sqrt(dx_const**2 + dy_const**2)
6194 dist = (i-icenter)*dx_const*(i-icenter)*dx_const &
6195 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6196 dist =
min(r0, sqrt(dist))
6198 prf = ak(k) + ps(i,j)*bk(k)
6199 if ( prf > 100.e2 )
then 6200 pt(i,j,k) = pt(i,j,k) + 0.01*(1. - (dist/r0)) * prf/ps(i,j)
6206 if ( hydrostatic )
then 6207 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6208 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6209 moist_phys, .true., nwat , domain)
6212 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6213 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6214 moist_phys, hydrostatic, nwat, domain, .true. )
6221 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6223 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6225 q(i,j,k,1) =
max(2.e-6, 0.8*pm(i)/ps(i,j)*qs(i) )
6241 if ( .not. hydrostatic ) w(:,:,:) = 0.
6253 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6262 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6268 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6269 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6270 moist_phys, .false., nwat, domain)
6273 r0 = 5.*
max(dx_const, dy_const)
6282 zm = ze - 0.5*delz(i,j,k)
6283 ze = ze - delz(i,j,k)
6284 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + &
6287 if ( dist <= r0 )
then 6288 pt(i,j,k) = pt(i,j,k) + 5.*(1.-dist/r0)
6311 ze1(k) = ze1(k+1) + ztop/
real(npz)
6325 delz(i,j,k) = ze1(k+1) - ze1(k)
6326 pk(i,j,k) = pk(i,j,k+1) +
grav*delz(i,j,k)/(
cp_air*t00)*pk0
6327 pe(i,k,j) = pk(i,j,k)**(1./
kappa)
6333 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
6338 peln(i,k,j) = log(pe(i,k,j))
6347 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
6348 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6359 zm = (0.5*(ze1(k)+ze1(k+1))-3.e3) / 2.e3
6363 xx = (dx_const * (0.5+
real(i-1)) - xc) / 4.e3
6364 yy = (dy_const * (0.5+
real(j-1)) - xc) / 4.e3
6365 dist = sqrt( xx**2 + yy**2 + zm**2 )
6366 if ( dist<=1. )
then 6367 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(
pi*dist)+1.)/2.
6370 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
6385 pk(i,j,1) = ptop**
kappa 6387 peln(i,1,j) = log(ptop)
6394 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6395 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6396 peln(i,k+1,j) = log(pe(i,k+1,j))
6397 pk(i,j,k+1) = exp(
kappa*peln(i,k+1,j) )
6405 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
6418 delz(i,j,k) =
rdgas/
grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6425 ze1(k) = ze1(k+1) - delz(is,js,k)
6429 zm = 0.5*(ze1(k)+ze1(k+1))
6430 utmp = us0*tanh(zm/3.e3)
6438 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6439 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6440 .true., hydrostatic, nwat, domain)
6446 icenter = (npx-1)/3 + 1
6447 jcenter = (npy-1)/2 + 1
6449 zm = 0.5*(ze1(k)+ze1(k+1))
6450 ptmp = ( (zm-zc)/zc ) **2
6451 if ( ptmp < 1. )
then 6454 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6455 if ( dist < 1. )
then 6456 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6474 pk(i,j,1) = ptop**
kappa 6476 peln(i,1,j) = log(ptop)
6483 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6484 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6485 peln(i,k+1,j) = log(pe(i,k+1,j))
6486 pk(i,j,k+1) = exp(
kappa*peln(i,k+1,j) )
6494 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
6506 delz(i,j,k) =
rdgas/
grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6513 ze1(k) = ze1(k+1) - delz(is,js,k)
6519 zm = 0.5*(ze1(k)+ze1(k+1))
6520 if ( zm .le. 2.e3 )
then 6521 utmp = 8.*(1.-cos(
pi*zm/4.e3))
6522 vtmp = 8.*sin(
pi*zm/4.e3)
6523 elseif (zm .le. 6.e3 )
then 6524 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
6533 u(i,j,k) = utmp - 8.
6539 v(i,j,k) = vtmp - 4.
6545 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6546 pe, peln, pk, pkz,
kappa, q,
ng, ncnst, area, dry_mass, .false., .false., &
6547 .true., hydrostatic, nwat, domain)
6553 icenter = (npx-1)/2 + 1
6554 jcenter = (npy-1)/2 + 1
6556 zm = 0.5*(ze1(k)+ze1(k+1))
6557 ptmp = ( (zm-zc)/zc ) **2
6558 if ( ptmp < 1. )
then 6561 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6562 if ( dist < 1. )
then 6563 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6585 if (.not.hybrid_z)
call mpp_error(fatal,
'hybrid_z must be .TRUE.')
6590 call mpp_error(fatal,
'npz must be == 101 ')
6603 peln(i,npz+1,j) = log(p00)
6610 peln(i,k,j) = peln(i,k+1,j) +
grav*delz(i,j,k)/(
rdgas*t00)
6611 pe(i,k,j) = exp(peln(i,k,j))
6612 pk(i,j,k) = pe(i,k,j)**
kappa 6621 if ( is_master() )
write(*,*)
'LES testcase: computed model top (mb)=', ptop/100.
6626 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
6627 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6635 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6637 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6639 if ( pm(i) > 100.e2 )
then 6640 q(i,j,k,1) = 0.9*qs(i)
6657 zm = 0.5*(ze0(i,j,k)+ze0(i,j,k+1))
6658 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + (zm-zc)**2
6660 if ( dist <= r0 )
then 6661 pt(i,j,k) = pt(i,j,k) + 2.0*(1.-dist/r0)
6699 nullify(have_south_pole)
6700 nullify(have_north_pole)
6712 integer,
intent(in):: km
6713 real,
intent(in):: p00
6714 real,
intent(inout),
dimension(km+1):: pe
6715 real,
intent(in),
dimension(km+1):: ze
6718 real,
intent(out),
dimension(km):: pt, qz
6720 integer,
parameter:: nx = 5
6721 real,
parameter:: qst = 1.0e-6
6722 real,
parameter:: qv0 = 1.4e-2
6723 real,
parameter:: ztr = 12.e3
6724 real,
parameter:: ttr = 213.
6725 real,
parameter:: ptr = 343.
6726 real,
parameter:: pt0 = 300.
6727 real,
dimension(km):: zs, rh, temp, dp, dp0
6728 real,
dimension(km+1):: peln, pk
6729 real:: qs, zvir, fac_z, pk0, temp1, pm
6734 if ( (is_master()) )
then 6735 write(*,*)
'Computing sounding for HIWPP super-cell test using p00=', p00
6742 zs(k) = 0.5*(ze(k)+ze(k+1))
6744 if ( zs(k) .gt. ztr )
then 6746 pt(k) = ptr*exp(
grav*(zs(k)-ztr)/(
cp_air*ttr))
6749 fac_z = (zs(k)/ztr)**1.25
6750 pt(k) = pt0 + (ptr-pt0)* fac_z
6751 rh(k) = 1. - 0.75 * fac_z
6753 qz(k) = qv0 - (qv0-qst)*fac_z
6755 if ( is_master() )
write(*,*) zs(k), pt(k), qz(k)
6760 #ifdef USE_MOIST_P00 6767 peln(km+1) = log(p00)
6772 pk(k) = pk(k+1) -
grav*(ze(k)-ze(k+1))/(
cp_air*pt(k)*(1.+zvir*qz(k)))
6773 peln(k) = log(pk(k)) /
kappa 6774 pe(k) = exp(peln(k))
6777 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
6778 temp(k) = pt(k)*pm**
kappa 6780 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6781 qz(k) =
min( qv0, rh(k)*qs )
6782 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
6789 peln(km+1) = log(p00)
6793 pk(k) = pk(k+1) -
grav*(ze(k)-ze(k+1))/(
cp_air*pt(k))
6794 peln(k) = log(pk(k)) /
kappa 6795 pe(k) = exp(peln(k))
6798 dp0(k) = pe(k+1) - pe(k)
6799 pm = dp0(k)/(peln(k+1)-peln(k))
6800 temp(k) = pt(k)*pm**
kappa 6802 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6803 qz(k) =
min( qv0, rh(k)*qs )
6809 dp(k) = dp0(k)*(1. + qz(k))
6810 pe(k+1) = pe(k) + dp(k)
6813 pk(km+1) = pe(km+1)**
kappa 6814 peln(km+1) = log(pe(km+1))
6818 pk(k) = pk(k+1) -
grav*(ze(k)-ze(k+1))/(
cp_air*pt(k)*(1.+zvir*qz(k)))
6819 peln(k) = log(pk(k)) /
kappa 6820 pe(k) = exp(peln(k))
6823 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
6824 temp(k) = pt(k)*pm**
kappa 6826 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
6827 qz(k) =
min( qv0, rh(k)*qs )
6828 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
6833 if ( is_master() )
then 6834 write(*,*)
'Super_K: computed ptop (mb)=', 0.01*pe(1),
' PS=', 0.01*pe(km+1)
6835 call prt_m1(
'1D Sounding T0', temp, 1, km, 1, 1, 0, 1, 1.)
6840 subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
6841 delz, zvir, ptop, ak, bk, agrid)
6842 integer,
intent(in):: is, ie, js, je, ng, km
6843 real,
intent(in),
dimension(km ):: ts1, qs1, uz1, dudz
6844 real,
intent(in),
dimension(km+1):: ze1
6845 real,
intent(in):: zvir, ps0
6846 real,
intent(inout):: ptop
6847 real(kind=R_GRID),
intent(in):: agrid(is-ng:ie+ng,js-ng:je+ng,2)
6848 real,
intent(inout),
dimension(km+1):: ak, bk
6849 real,
intent(inout),
dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delz
6850 real,
intent(out),
dimension(is:ie,js:je,km+1):: pk
6852 real,
intent(inout),
dimension(is-1:ie+1,km+1,js-1:je+1):: pe
6854 integer,
parameter:: nt=5
6855 integer,
parameter:: nlat=1001
6856 real,
dimension(nlat,km):: pt2, pky, dzc
6857 real,
dimension(nlat,km+1):: pk2, pe2, peln2, pte
6858 real,
dimension(km+1):: pe1
6859 real:: lat(nlat), latc(nlat-1)
6860 real:: fac_y, dlat, dz0, pk0, tmp1, tmp2, tmp3, pint
6861 integer::i,j,k,n, jj, k1
6865 dz0 = ze1(km) - ze1(km+1)
6868 dlat = 0.5*
pi/
real(nlat-1)
6870 lat(j) = dlat*
real(j-1)
6872 dzc(j,k) = ze1(k) - ze1(k+1)
6876 latc(j) = 0.5*(lat(j)+lat(j+1))
6885 if ( is_master() )
then 6887 call prt_m1(
'Super_K PT0', pt2, 1, nlat, 1, km, 0, 1, tmp1)
6894 call ppme(pt2, pte, dzc, nlat, km)
6897 tmp1 = 0.5*(pte(j-1,k ) + pte(j,k ))
6898 tmp3 = 0.5*(pte(j-1,k+1) + pte(j,k+1))
6899 pt2(j,k) = pt2(j-1,k) + dlat/(2.*
grav)*sin(2.*latc(j-1))*uz1(k)* &
6900 ( uz1(k)*(tmp1-tmp3)/dzc(j,k) - (pt2(j-1,k)+pt2(j,k))*dudz(k) )
6903 if ( is_master() )
then 6904 call prt_m1(
'Super_K PT', pt2, 1, nlat, 1, km, 0, 1, pk0/
cp_air)
6910 pk2(1,km+1) = ps0**
kappa 6912 pk2(j,km+1) = pk2(j-1,km+1) - dlat*uz1(km)*uz1(km)*sin(2.*latc(j-1)) &
6913 / (pt2(j-1,km) + pt2(j,km))
6918 pk2(j,k) = pk2(j,k+1) -
grav*dzc(j,k)/pt2(j,k)
6924 peln2(j,k) = log(pk2(j,k)) /
kappa 6925 pe2(j,k) = exp(peln2(j,k))
6931 pky(j,k) = (pk2(j,k+1)-pk2(j,k))/(
kappa*(peln2(j,k+1)-peln2(j,k)))
6932 pt2(j,k) = pt2(j,k)*pky(j,k)/(
cp_air*(1.+zvir*qs1(k)))
6940 if ( is_master() )
then 6941 write(*,*)
'SuperK ptop at EQ=', 0.01*pe1(1),
'new ptop=', 0.01*ptop
6942 call prt_m1(
'Super_K pe', pe2, 1, nlat, 1, km+1, 0, 1, 0.01)
6943 call prt_m1(
'Super_K Temp', pt2, 1, nlat, 1, km, 0, 1, 1.)
6950 if (abs(agrid(i,j,2))>=lat(jj) .and. abs(agrid(i,j,2))<=lat(jj+1) )
then 6952 fac_y = (abs(agrid(i,j,2))-lat(jj)) / dlat
6954 pt(i, j,k) = pt2(jj, k) + fac_y*(pt2(jj+1, k)-pt2(jj,k))
6957 pe(i,k,j) = pe2(jj,k) + fac_y*(pe2(jj+1,k)-pe2(jj,k))
6980 bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint)
6981 ak(k) = pe1(k) - bk(k) * pe1(km+1)
6982 if ( is_master() )
write(*,*) k, ak(k), bk(k)
6995 subroutine superk_u(km, zz, um, dudz)
6996 integer,
intent(in):: km
6997 real,
intent(in):: zz(km)
6998 real,
intent(out):: um(km), dudz(km)
7000 real,
parameter:: zs = 5.e3
7001 real,
parameter:: us = 30.
7008 if ( zz(k) .gt. zs+1.e3 )
then 7011 elseif ( abs(zz(k)-zs) .le. 1.e3 )
then 7012 um(k) = us*(-4./5. + 3.*zz(k)/zs - 5./4.*(zz(k)/zs)**2)
7013 dudz(k) = us/zs*(3. - 5./2.*zz(k)/zs)
7022 um(k) = us*tanh( zz(k)/zs ) - uc
7023 dudz(k) = (us/zs)/cosh(zz(k)/zs)**2
7031 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7032 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7033 hydrostatic, nwat, adiabatic, do_pert, domain)
7035 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7036 real,
intent(IN) :: ptop
7037 real,
intent(IN),
dimension(npz+1) :: ak, bk
7038 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7039 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w, delz
7040 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7041 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7042 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7043 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7044 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7045 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7046 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7047 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7048 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7049 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7050 logical,
intent(IN) :: hydrostatic,adiabatic,do_pert
7051 type(domain2d),
intent(INOUT) :: domain
7053 real,
parameter :: p0 = 1.e5
7054 real,
parameter :: u0 = 35.
7055 real,
parameter :: b = 2.
7056 real,
parameter :: KK = 3.
7057 real,
parameter :: Te = 310.
7058 real,
parameter :: Tp = 240.
7059 real,
parameter :: T0 = 0.5*(te + tp)
7060 real,
parameter :: up = 1.
7061 real,
parameter :: zp = 1.5e4
7062 real(kind=R_GRID),
parameter :: lamp =
pi/9.
7063 real(kind=R_GRID),
parameter :: phip = 2.*lamp
7064 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7065 real,
parameter :: Rp =
radius/10.
7066 real,
parameter :: lapse = 5.e-3
7067 real,
parameter :: dT = 4.8e5
7068 real,
parameter :: phiW = 2.*
pi/9.
7069 real,
parameter :: pW = 34000.
7070 real,
parameter :: q0 = .018
7071 real,
parameter :: qt = 1.e-12
7072 real,
parameter :: ptrop = 1.e4
7074 real,
parameter :: zconv = 1.e-6
7079 integer :: i,j,k,iter, sphum, cl, cl2, n
7080 real :: p,z,z0,ziter,piter,titer,uu,vv,pl,pt_u,pt_v
7081 real(kind=R_GRID),
dimension(2) :: pa
7082 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7083 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2
7084 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7085 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2
7086 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7105 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7116 peln(i,1,j) = log(ptop)
7117 pk(i,j,1) = ptop**
kappa 7121 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7124 pk(i,j,k) = exp(
kappa*log(pe(i,k,j)))
7125 peln(i,k,j) = log(pe(i,k,j))
7133 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
7155 z = ziter + (piter - p)*rdgrav*titer/piter
7161 if (abs(z - ziter) < zconv)
exit 7172 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7183 peln_u(i,j) = log(p0)
7198 p = ak(k) + ps_u(i,j)*bk(k)
7207 z = ziter + (piter - p)*rdgrav*titer/piter
7208 if (abs(z - ziter) < zconv)
exit 7211 pt_u = rrdgrav * ( z - gz_u(i,j) ) / (peln_u(i,j) - pl)
7218 u(i,j,k) = u1(i,j)*uu
7231 peln_v(i,j) = log(p0)
7246 p = ak(k) + ps_v(i,j)*bk(k)
7255 z = ziter + (piter - p)*rdgrav*titer/piter
7256 if (abs(z - ziter) < zconv)
exit 7259 pt_v = rrdgrav * ( z - gz_v(i,j) ) / (peln_v(i,j) - pl)
7265 v(i,j,k) = v1(i,j)*uu
7283 if (.not. adiabatic)
then 7288 p = delp(i,j,k)/(peln(i,k+1,j) - peln(i,k,j))
7289 q(i,j,k,sphum) =
dcmip16_bc_sphum(p,ps(i,j),agrid(i,j,2),agrid(i,j,1))
7291 pt(i,j,k) = pt(i,j,k) / ( 1. + zvir*q(i,j,k,sphum))
7299 if (cl > 0 .and. cl2 > 0)
then 7301 q, delp,nq,agrid(isd,jsd,1),agrid(isd,jsd,2))
7306 if (.not. hydrostatic)
then 7311 delz(i,j,k) = gz(i,j,k) - gz(i,j,k+1)
7322 real,
intent(IN) :: z
7323 real(kind=R_GRID),
intent(IN) :: lat
7324 real :: it, t1, t2, tr, zsc
7326 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7328 tr = ( 1. - 2.*zsc**2.) * exp(-zsc**2. )
7330 t1 = (1./t0)*exp(lapse*z/t0) + (t0 - tp)/(t0*tp) * tr
7331 t2 = 0.5* ( kk + 2.) * (te - tp)/(te*tp) * tr
7339 real,
intent(IN) :: z
7340 real(kind=R_GRID),
intent(IN) :: lat
7341 real :: it, ti1, ti2, tir
7343 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7346 ti1 = 1./lapse* (exp(lapse*z/t0) - 1.) + tir*(t0-tp)/(t0*tp)
7347 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7355 real,
intent(IN) :: z, t
7356 real(kind=R_GRID),
intent(IN) :: lat
7357 real :: tir, ti2, uu, ur
7360 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7362 uu =
grav*kk/
radius * ti2 * ( cos(lat)**(int(kk)-1) - cos(lat)**(int(kk)+1) ) * t
7371 real,
intent(IN) :: z
7372 real(kind=R_GRID),
intent(IN) :: lat, lon
7374 real(kind=R_GRID) :: dst, pphere(2)
7377 zz =
max(1. - 3.*zrat*zrat + 2.*zrat*zrat*zrat, 0.)
7379 pphere = (/ lon, lat /)
7388 real,
intent(IN) :: p, ps
7389 real(kind=R_GRID),
intent(IN) :: lat, lon
7404 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7405 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7406 hydrostatic, nwat, adiabatic)
7408 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7409 real,
intent(IN) :: ptop
7410 real,
intent(IN),
dimension(npz+1) :: ak, bk
7411 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7412 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w, delz
7413 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7414 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7415 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7416 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7417 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7418 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7419 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7420 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7421 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7422 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7423 logical,
intent(IN) :: hydrostatic,adiabatic
7425 real,
parameter :: zt = 15000
7426 real,
parameter :: q0 = 0.021
7427 real,
parameter :: qt = 1.e-11
7428 real,
parameter :: T0 = 302.15
7429 real,
parameter :: Tv0 = 302.15*(1.+0.608*q0)
7430 real,
parameter :: Ts = 302.15
7431 real,
parameter :: zq1 = 3000.
7432 real,
parameter :: zq2 = 8000.
7433 real,
parameter :: lapse = 7.e-3
7434 real,
parameter :: Tvt = tv0 - lapse*zt
7435 real,
parameter :: pb = 101500.
7436 real,
parameter :: ptt = pb*(tvt/tv0)**(
grav/
rdgas/lapse)
7437 real(kind=R_GRID),
parameter :: lamp =
pi 7438 real(kind=R_GRID),
parameter :: phip =
pi/18.
7439 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7440 real,
parameter :: dp = 1115.
7441 real,
parameter :: rp = 282000.
7442 real,
parameter :: zp = 7000.
7443 real,
parameter :: fc = 2.*
omega*sin(phip)
7445 real,
parameter :: zconv = 1.e-6
7449 integer :: i,j,k,iter, sphum, cl, cl2, n
7450 real :: p,z,z0,ziter,piter,titer,uu,vv,pl, r
7451 real(kind=R_GRID),
dimension(2) :: pa
7452 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7453 real,
dimension(is:ie,js:je) :: rc
7454 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2, rc_u
7455 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7456 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2, rc_v
7457 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7475 ps(i,j) = pb - dp*exp( -sqrt((rc(i,j)/rp)**3) )
7483 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7494 peln(i,1,j) = log(ptop)
7495 pk(i,j,1) = ptop**
kappa 7499 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7502 pk(i,j,k) = exp(
kappa*log(pe(i,k,j)))
7503 peln(i,k,j) = log(pe(i,k,j))
7511 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(
kappa*(peln(i,k+1,j)-peln(i,k,j)))
7533 z = ziter + (piter - p)*rdgrav*titer/piter
7539 if (abs(z - ziter) < zconv)
exit 7550 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7568 p_u(i,j) = pb - dp*exp( -sqrt((rc_u(i,j)/rp)**3) )
7569 peln_u(i,j) = log(p_u(i,j))
7570 ps_u(i,j) = p_u(i,j)
7577 p = ak(k) + ps_u(i,j)*bk(k)
7586 z = ziter + (piter - p)*rdgrav*titer/piter
7587 if (abs(z - ziter) < zconv)
exit 7591 u(i,j,k) = u1(i,j)*uu + u2(i,j)*vv
7611 p_v(i,j) = pb - dp*exp( - sqrt((rc_v(i,j)/rp)**3) )
7612 peln_v(i,j) = log(p_v(i,j))
7613 ps_v(i,j) = p_v(i,j)
7620 p = ak(k) + ps_v(i,j)*bk(k)
7629 z = ziter + (piter - p)*rdgrav*titer/piter
7630 if (abs(z - ziter) < zconv)
exit 7634 v(i,j,k) = v1(i,j)*uu + v2(i,j)*vv
7652 if (.not. adiabatic)
then 7657 z = 0.5*(gz(i,j,k) + gz(i,j,k+1))
7665 if (.not. hydrostatic)
then 7670 delz(i,j,k) = gz(i,j,k) - gz(i,j,k+1)
7681 real,
intent(IN) :: z, r
7682 real :: tv, term1, term2
7690 term1 =
grav*zp*zp* ( 1. - pb/dp * exp( sqrt(r/rp)**3 + (z/zp)**2 ) )
7691 term2 = 2*
rdgas*tv*z
7699 real,
intent(IN) :: z, r
7703 exp(
grav/(
rdgas*lapse) * log( (tv0-lapse*z)/tv0) )
7712 real,
intent(IN) :: z, r
7713 real(kind=R_GRID),
intent(IN) :: lon, lat
7714 real,
intent(OUT) :: uu, vv
7715 real :: rfac, Tvrd, vt, fr5, d1, d2, d
7716 real(kind=R_GRID) :: dst, pphere(2)
7724 rfac = sqrt(r/rp)**3
7727 tvrd = (tv0 - lapse*z)*
rdgas 7729 vt = -fr5 + sqrt( fr5**2 - (1.5 * rfac * tvrd) / &
7730 ( 1. + 2*tvrd*z/(
grav*zp**2) - pb/dp*exp( rfac + (z/zp)**2) ) )
7732 d1 = sin(phip)*cos(lat) - cos(phip)*sin(lat)*cos(lon - lamp)
7733 d2 = cos(phip)*sin(lon - lamp)
7734 d =
max(1.e-25,sqrt(d1*d1 + d2*d2))
7743 real,
intent(IN) :: z
7754 subroutine init_latlon(u,v,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
7755 gridstruct, npx, npy, npz, ng, ncnst, ndims, nregions, dry_mass, &
7756 mountain, moist_phys, hybrid_z, delz, ze0, domain_in, tile_in)
7758 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
7759 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
7760 real ,
intent(INOUT) :: pt(isd:ied ,jsd:jed ,npz)
7761 real ,
intent(INOUT) :: delp(isd:ied ,jsd:jed ,npz)
7762 real ,
intent(INOUT) :: q(isd:ied ,jsd:jed ,npz, ncnst)
7764 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
7766 real ,
intent(INOUT) :: ps(isd:ied ,jsd:jed )
7767 real ,
intent(INOUT) :: pe(is-1:ie+1,npz+1,js-1:je+1)
7768 real ,
intent(INOUT) :: pk(is:ie ,js:je ,npz+1)
7769 real ,
intent(INOUT) :: peln(is :ie ,npz+1 ,js:je)
7770 real ,
intent(INOUT) :: pkz(is:ie ,js:je ,npz )
7771 real ,
intent(INOUT) :: uc(isd:ied+1,jsd:jed ,npz)
7772 real ,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1,npz)
7773 real ,
intent(INOUT) :: ua(isd:ied ,jsd:jed ,npz)
7774 real ,
intent(INOUT) :: va(isd:ied ,jsd:jed ,npz)
7775 real ,
intent(inout) :: delz(isd:,jsd:,1:)
7776 real ,
intent(inout) :: ze0(is:,js:,1:)
7778 real ,
intent(IN) :: ak(npz+1)
7779 real ,
intent(IN) :: bk(npz+1)
7781 integer,
intent(IN) :: npx, npy, npz
7782 integer,
intent(IN) ::
ng, ncnst
7783 integer,
intent(IN) :: ndims
7784 integer,
intent(IN) :: nregions
7785 integer,
target,
intent(IN):: tile_in
7787 real,
intent(IN) :: dry_mass
7788 logical,
intent(IN) :: mountain
7789 logical,
intent(IN) :: moist_phys
7790 logical,
intent(IN) :: hybrid_z
7793 type(
domain2d),
intent(IN),
target :: domain_in
7795 real,
pointer,
dimension(:,:,:) ::
agrid, grid
7796 real,
pointer,
dimension(:,:) :: area, rarea, fc,
f0 7797 real,
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
7798 real,
pointer,
dimension(:,:,:,:) :: ew, es
7799 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
7801 logical,
pointer :: cubed_sphere, latlon
7804 integer,
pointer :: tile
7806 logical,
pointer :: have_south_pole, have_north_pole
7808 integer,
pointer :: ntiles_g
7809 real,
pointer :: acapn, acaps, globalarea
7811 real(kind=R_GRID) :: p1(2), p2(2)
7815 agrid => gridstruct%agrid
7816 grid => gridstruct%grid
7818 area => gridstruct%area
7822 dxa => gridstruct%dxa
7823 dya => gridstruct%dya
7824 rdxa => gridstruct%rdxa
7825 rdya => gridstruct%rdya
7826 dxc => gridstruct%dxc
7827 dyc => gridstruct%dyc
7832 ntiles_g => gridstruct%ntiles_g
7833 acapn => gridstruct%acapN
7834 acaps => gridstruct%acapS
7835 globalarea => gridstruct%globalarea
7840 have_south_pole => gridstruct%have_south_pole
7841 have_north_pole => gridstruct%have_north_pole
7845 fc(i,j) = 2.*
omega*( -cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) &
7846 +sin(grid(i,j,2))*cos(
alpha) )
7867 p2(1) =
agrid(i,j,1)
7868 p2(2) =
agrid(i,j,2)
7871 delp(i,j,1) = phis(i,j) + 0.5*(1.0+cos(
pi*r/r0))
7873 delp(i,j,1) = phis(i,j)
7924 nullify(have_south_pole)
7925 nullify(have_north_pole)
7938 real,
intent(INOUT) :: UBar
7939 real,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
7940 real,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
7941 real,
intent(INOUT) :: uc(isd:ied+1,jsd:jed )
7942 real,
intent(INOUT) :: vc(isd:ied ,jsd:jed+1)
7943 real,
intent(INOUT) :: ua(isd:ied ,jsd:jed )
7944 real,
intent(INOUT) :: va(isd:ied ,jsd:jed )
7945 integer,
intent(IN) :: defOnGrid
7946 type(fv_grid_type),
intent(IN),
target :: gridstruct
7948 real :: p1(2),p2(2),p3(2),p4(2), pt(2)
7949 real :: e1(3), e2(3), ex(3), ey(3)
7955 real :: psi_b(isd:ied+1,jsd:jed+1), psi(isd:ied,jsd:jed), psi1, psi2
7957 real,
dimension(:,:,:),
pointer :: grid, agrid
7958 real,
dimension(:,:),
pointer :: area, dx, dy, dxc, dyc
7960 grid => gridstruct%grid
7961 agrid=> gridstruct%agrid
7963 area => gridstruct%area
7966 dxc => gridstruct%dxc
7967 dyc => gridstruct%dyc
7973 psi(i,j) = (-1.0 * ubar *
radius *( sin(agrid(i,j,2)) *cos(
alpha) - &
7974 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
7979 psi_b(i,j) = (-1.0 * ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
7980 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
7984 if ( defongrid == 1 )
then 7988 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
7989 if (dist==0) vc(i,j) = 0.
7995 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
7996 if (dist==0) uc(i,j) = 0.
8004 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
8005 if (dist==0) v(i,j) = 0.
8011 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
8012 if (dist==0) u(i,j) = 0.
8019 subroutine d2a2c(im,jm,km, ifirst,ilast, jfirst,jlast, ng, nested, &
8020 u,v, ua,va, uc,vc, gridstruct, domain)
8023 integer,
intent(IN) :: im,jm,km
8024 integer,
intent(IN) :: ifirst,ilast
8025 integer,
intent(IN) :: jfirst,jlast
8026 integer,
intent(IN) :: ng
8027 logical,
intent(IN) :: nested
8028 type(fv_grid_type),
intent(IN),
target :: gridstruct
8029 type(domain2d),
intent(INOUT) :: domain
8044 real ,
intent(inout) :: u(isd:ied,jsd:jed+1)
8045 real ,
intent(inout) :: v(isd:ied+1,jsd:jed)
8046 real ,
intent(inout) :: ua(isd:ied,jsd:jed)
8047 real ,
intent(inout) :: va(isd:ied,jsd:jed)
8048 real ,
intent(inout) :: uc(isd:ied+1,jsd:jed)
8049 real ,
intent(inout) :: vc(isd:ied,jsd:jed+1)
8054 real :: sinlon(im,jm)
8055 real :: coslon(im,jm)
8056 real :: sinl5(im,jm)
8057 real :: cosl5(im,jm)
8059 real :: tmp1(jsd:jed+1)
8060 real :: tmp2(jsd:jed)
8061 real :: tmp3(jsd:jed)
8063 real mag,mag1,mag2, ang,ang1,ang2
8065 integer i, j, k, im2
8079 real,
pointer,
dimension(:,:,:) :: agrid, grid
8080 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
8081 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
8082 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
8083 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
8085 logical,
pointer :: cubed_sphere, latlon
8087 logical,
pointer :: have_south_pole, have_north_pole
8089 integer,
pointer :: ntiles_g
8090 real,
pointer :: acapN, acapS, globalarea
8092 grid => gridstruct%grid
8093 agrid=> gridstruct%agrid
8095 area => gridstruct%area
8096 rarea => gridstruct%rarea
8101 ee1 => gridstruct%ee1
8102 ee2 => gridstruct%ee2
8105 en1 => gridstruct%en1
8106 en2 => gridstruct%en2
8110 dxa => gridstruct%dxa
8111 dya => gridstruct%dya
8112 rdxa => gridstruct%rdxa
8113 rdya => gridstruct%rdya
8114 dxc => gridstruct%dxc
8115 dyc => gridstruct%dyc
8117 cubed_sphere => gridstruct%cubed_sphere
8118 latlon => gridstruct%latlon
8120 have_south_pole => gridstruct%have_south_pole
8121 have_north_pole => gridstruct%have_north_pole
8123 ntiles_g => gridstruct%ntiles_g
8124 acapn => gridstruct%acapN
8125 acaps => gridstruct%acapS
8126 globalarea => gridstruct%globalarea
8128 if (cubed_sphere)
then 8130 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,im,jm,ng)
8131 if (.not. nested)
call fill_corners(ua, va, im, jm, vector=.true., agrid=.true.)
8132 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,im,jm,ng, nested, domain, nocomm=.true.)
8133 if (.not. nested)
call fill_corners(uc, vc, im, jm, vector=.true., cgrid=.true.)
8145 js2gcp1 = jfirst-ng-1
8151 jn2gsp1 = jlast+ng-1
8153 if (have_south_pole)
then 8161 if (have_north_pole)
then 8171 if ( ng == 1 .AND. ng > 1 )
THEN 8174 js2gc1 = jfirst-ng+1
8175 if (have_south_pole) js2gc1 = 2
8180 if ((have_south_pole) .or. (have_north_pole))
then 8182 call vpol5(u(1:im,:), v(1:im,:), im, jm, &
8183 coslon, sinlon, cosl5, sinl5, ng, ng, jfirst, jlast )
8184 call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, v(:,:))
8187 call dtoa(u, v, ua, va, dx,dy,dxa,dya,dxc,dyc,im, jm, ng)
8188 if (.not. nested)
call fill_corners(ua, va, im, jm, vector=.true., agrid=.true.)
8190 if ( have_south_pole )
then 8195 us = us + (ua(i+im2,2)-ua(i,2))*sinlon(i,2) &
8196 + (va(i,2)-va(i+im2,2))*coslon(i,2)
8197 vs = vs + (ua(i+im2,2)-ua(i,2))*coslon(i,2) &
8198 + (va(i+im2,2)-va(i,2))*sinlon(i,2)
8204 ua(i,1) = -us*sinlon(i,1) - vs*coslon(i,1)
8205 va(i,1) = us*coslon(i,1) - vs*sinlon(i,1)
8206 ua(i+im2,1) = -ua(i,1)
8207 va(i+im2,1) = -va(i,1)
8210 ua(im+1,1) = ua(1 ,1)
8211 va(im+1,1) = va(1 ,1)
8214 if ( have_north_pole )
then 8220 un = un + (ua(i+im2,j)-ua(i,j))*sinlon(i,j) &
8221 + (va(i+im2,j)-va(i,j))*coslon(i,j)
8222 vn = vn + (ua(i,j)-ua(i+im2,j))*coslon(i,j) &
8223 + (va(i+im2,j)-va(i,j))*sinlon(i,j)
8229 ua(i,jm) = -un*sinlon(i,jm) + vn*coslon(i,jm)
8230 va(i,jm) = -un*coslon(i,jm) - vn*sinlon(i,jm)
8231 ua(i+im2,jm) = -ua(i,jm)
8232 va(i+im2,jm) = -va(i,jm)
8234 ua(0 ,jm) = ua(im,jm)
8235 ua(im+1,jm) = ua(1 ,jm)
8236 va(im+1,jm) = va(1 ,jm)
8239 if (latlon)
call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, ua(:,:))
8240 if (latlon)
call mp_ghost_ew(im,jm,1,1, ifirst,ilast, jfirst,jlast, 1,1, ng,ng, ng,ng, va(:,:))
8243 call atoc(ua, va, uc, vc, dx,dy,dxa,dya,im, jm, ng, nested, domain, nocomm=.true.)
8247 if (.not. nested)
call fill_corners(uc, vc, im, jm, vector=.true., cgrid=.true.)
8251 end subroutine d2a2c 8254 subroutine atob_s(qin, qout, npx, npy, dxa, dya, nested, cubed_sphere, altInterp)
8258 integer,
intent(IN) :: npx, npy
8259 real ,
intent(IN) :: qin(isd:ied ,jsd:jed )
8260 real ,
intent(OUT) :: qout(isd:ied+1,jsd:jed+1)
8261 integer,
OPTIONAL,
intent(IN) :: altInterp
8262 logical,
intent(IN) :: nested, cubed_sphere
8263 real,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8267 real :: tmp1j(jsd:jed+1)
8268 real :: tmp2j(jsd:jed+1)
8269 real :: tmp3j(jsd:jed+1)
8270 real :: tmp1i(isd:ied+1)
8271 real :: tmp2i(isd:ied+1)
8272 real :: tmp3i(isd:ied+1)
8273 real :: tmpq(isd:ied ,jsd:jed )
8274 real :: tmpq1(isd:ied+1,jsd:jed+1)
8275 real :: tmpq2(isd:ied+1,jsd:jed+1)
8277 if (
present(altinterp))
then 8279 tmpq(:,:) = qin(:,:)
8281 if (.not. nested)
call fill_corners(tmpq , npx, npy, fill=xdir,
agrid=.true.)
8287 if (.not. nested)
call fill_corners(tmpq , npx, npy, fill=ydir,
agrid=.true.)
8290 tmp1j(jsd:jed) = 0.0
8291 tmp2j(jsd:jed) = tmpq(i,jsd:jed)
8292 tmp3j(jsd:jed) = dya(i,jsd:jed)
8294 tmpq2(i,jsd:jed) = tmp1j(jsd:jed)
8299 tmp1j(:) = tmpq1(i,:)
8300 tmp2j(:) = tmpq1(i,:)
8303 tmpq1(i,:) = tmp1j(:)
8308 tmp1i(:) = tmpq2(:,j)
8309 tmp2i(:) = tmpq2(:,j)
8312 tmpq2(:,j) = tmp1i(:)
8318 qout(i,j) = 0.5 * (tmpq1(i,j) + tmpq2(i,j))
8323 if (cubed_sphere .and. .not. nested)
then 8326 if ( (is==i) .and. (js==j) )
then 8327 qout(i,j) = (1./3.) * (qin(i,j) + qin(i-1,j) + qin(i,j-1))
8332 if ( (ie+1==i) .and. (js==j) )
then 8333 qout(i,j) = (1./3.) * (qin(i-1,j) + qin(i-1,j-1) + qin(i,j))
8338 if ( (is==i) .and. (je+1==j) )
then 8339 qout(i,j) = (1./3.) * (qin(i,j-1) + qin(i-1,j-1) + qin(i,j))
8344 if ( (ie+1==i) .and. (je+1==j) )
then 8345 qout(i,j) = (1./3.) * (qin(i-1,j-1) + qin(i,j-1) + qin(i-1,j))
8353 qout(i,j) = 0.25 * (qin(i-1,j) + qin(i-1,j-1) + &
8354 qin(i ,j) + qin(i ,j-1))
8358 if (.not. nested)
then 8361 if ( (is==i) .and. (js==j) )
then 8362 qout(i,j) = (1./3.) * (qin(i,j) + qin(i-1,j) + qin(i,j-1))
8367 if ( (ie+1==i) .and. (js==j) )
then 8368 qout(i,j) = (1./3.) * (qin(i-1,j) + qin(i-1,j-1) + qin(i,j))
8373 if ( (is==i) .and. (je+1==j) )
then 8374 qout(i,j) = (1./3.) * (qin(i,j-1) + qin(i-1,j-1) + qin(i,j))
8379 if ( (ie+1==i) .and. (je+1==j) )
then 8380 qout(i,j) = (1./3.) * (qin(i-1,j-1) + qin(i,j-1) + qin(i-1,j))
8396 subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, nested, domain)
8399 integer,
intent(IN) :: npx, npy, ng
8400 real ,
intent(IN) :: uin(isd:ied ,jsd:jed )
8401 real ,
intent(IN) :: vin(isd:ied ,jsd:jed )
8402 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed+1)
8403 real ,
intent(OUT) :: vout(isd:ied+1,jsd:jed )
8404 logical,
intent(IN) :: nested
8405 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8406 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dxc
8407 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dyc
8408 type(domain2d),
intent(INOUT) :: domain
8412 real :: tmp1i(isd:ied+1)
8413 real :: tmp2i(isd:ied)
8414 real :: tmp3i(isd:ied)
8415 real :: tmp1j(jsd:jed+1)
8416 real :: tmp2j(jsd:jed)
8417 real :: tmp3j(jsd:jed)
8421 tmp2i(:) = vin(:,j)*dxa(:,j)
8424 vout(:,j) = tmp1i(:)/dxc(:,j)
8428 tmp2j(:) = uin(i,:)*dya(i,:)
8431 uout(i,:) = tmp1j(:)/dyc(i,:)
8434 if (.not. nested)
call fill_corners(uout, vout, npx, npy, vector=.true., dgrid=.true.)
8445 subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng)
8447 integer,
intent(IN) :: npx, npy, ng
8448 real ,
intent(IN) :: uin(isd:ied ,jsd:jed+1)
8449 real ,
intent(IN) :: vin(isd:ied+1,jsd:jed )
8450 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed )
8451 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed )
8452 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dx, dyc
8453 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dy, dxc
8454 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8458 real :: tmp1i(isd:ied+1)
8459 real :: tmp2i(isd:ied+1)
8460 real :: tmp3i(isd:ied+1)
8461 real :: tmp1j(jsd:jed+1)
8462 real :: tmp2j(jsd:jed+1)
8463 real :: tmp3j(jsd:jed+1)
8470 uout(i,j) = 0.5*(uin(i,j)*dx(i,j)+uin(i,j+1)*dx(i,j+1))/dxa(i,j)
8471 vout(i,j) = 0.5*(vin(i,j)*dy(i,j)+vin(i+1,j)*dy(i+1,j))/dya(i,j)
8477 tmp2j(:) = uin(i,:)*dyc(i,:)
8480 uout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dya(i,jsd:jed)
8484 tmp2i(:) = vin(:,j)*dxc(:,j)
8487 vout(isd:ied,j) = tmp1i(isd+1:ied+1)/dxa(isd:ied,j)
8501 subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, nested, domain, noComm)
8504 integer,
intent(IN) :: npx, npy, ng
8505 real ,
intent(IN) :: uin(isd:ied ,jsd:jed )
8506 real ,
intent(IN) :: vin(isd:ied ,jsd:jed )
8507 real ,
intent(OUT) :: uout(isd:ied+1,jsd:jed )
8508 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed+1)
8509 logical,
intent(IN) :: nested
8510 logical,
OPTIONAL,
intent(IN) :: noComm
8511 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dx
8512 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dy
8513 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8514 type(domain2d),
intent(INOUT) :: domain
8519 real :: tmp1i(isd:ied+1)
8520 real :: tmp2i(isd:ied)
8521 real :: tmp3i(isd:ied)
8522 real :: tmp1j(jsd:jed+1)
8523 real :: tmp2j(jsd:jed)
8524 real :: tmp3j(jsd:jed)
8526 #if !defined(ALT_INTERP) 8531 uout(i,j) = ( uin(i,j)*dxa(i,j) + uin(i-1,j)*dxa(i-1,j) ) &
8532 / ( dxa(i,j) + dxa(i-1,j) )
8537 vout(i,j) = ( vin(i,j)*dya(i,j) + vin(i,j-1)*dya(i,j-1) ) &
8538 / ( dya(i,j) + dya(i,j-1) )
8550 vout(i,:) = tmp1j(:)
8557 tmp2i(:) = uin(:,j)*dya(:,j)
8560 uout(:,j) = tmp1i(:)/dy(:,j)
8564 tmp2j(:) = vin(i,:)*dxa(i,:)
8567 vout(i,:) = tmp1j(:)/dx(i,:)
8570 if (cubed_sphere .and. .not. nested)
then 8571 csfac = cos(30.0*
pi/180.0)
8573 if ( (is==1) .and. (js==1) )
then 8576 uout(i,j)=uout(i,j)*csfac
8577 uout(i,j-1)=uout(i,j-1)*csfac
8578 vout(i,j)=vout(i,j)*csfac
8579 vout(i-1,j)=vout(i-1,j)*csfac
8581 if ( (is==1) .and. (je==npy-1) )
then 8584 uout(i,j)=uout(i,j)*csfac
8585 uout(i,j+1)=uout(i,j+1)*csfac
8586 vout(i,j+1)=vout(i,j+1)*csfac
8587 vout(i-1,j+1)=vout(i-1,j+1)*csfac
8589 if ( (ie==npx-1) .and. (je==npy-1) )
then 8592 uout(i+1,j)=uout(i+1,j)*csfac
8593 uout(i+1,j+1)=uout(i+1,j+1)*csfac
8594 vout(i,j+1)=vout(i,j+1)*csfac
8595 vout(i+1,j+1)=vout(i+1,j+1)*csfac
8597 if ( (ie==npx-1) .and. (js==1) )
then 8600 uout(i+1,j)=uout(i+1,j)*csfac
8601 uout(i+1,j-1)=uout(i+1,j-1)*csfac
8602 vout(i,j)=vout(i,j)*csfac
8603 vout(i+1,j)=vout(i+1,j)*csfac
8609 if (
present(nocomm))
then 8610 if (.not. nocomm)
call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
8612 call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
8614 if (.not. nested)
call fill_corners(uout, vout, npx, npy, vector=.true., cgrid=.true.)
8626 subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng)
8629 integer,
intent(IN) :: npx, npy, ng
8630 real ,
intent(IN) :: uin(isd:ied+1,jsd:jed )
8631 real ,
intent(IN) :: vin(isd:ied ,jsd:jed+1)
8632 real ,
intent(OUT) :: uout(isd:ied ,jsd:jed )
8633 real ,
intent(OUT) :: vout(isd:ied ,jsd:jed )
8634 real ,
intent(IN),
dimension(isd:ied+1,jsd:jed) :: dxc, dy
8635 real ,
intent(IN),
dimension(isd:ied,jsd:jed+1) :: dyc, dx
8636 real ,
intent(IN),
dimension(isd:ied,jsd:jed) :: dxa, dya
8640 real :: tmp1i(isd:ied+1)
8641 real :: tmp2i(isd:ied+1)
8642 real :: tmp3i(isd:ied+1)
8643 real :: tmp1j(jsd:jed+1)
8644 real :: tmp2j(jsd:jed+1)
8645 real :: tmp3j(jsd:jed+1)
8659 tmp2j(:) = vin(i,:)*dx(i,:)
8662 vout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dxa(i,jsd:jed)
8666 tmp2i(:) = uin(:,j)*dy(:,j)
8669 uout(isd:ied,j) = tmp1i(isd+1:ied+1)/dya(isd:ied,j)
8682 subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
8685 integer,
intent(IN) :: ndims
8686 real ,
intent(INOUT) :: myU
8687 real ,
intent(INOUT) :: myV
8688 real(kind=R_GRID) ,
intent(IN) :: p1(ndims)
8689 real(kind=R_GRID) ,
intent(IN) :: p2(ndims)
8690 real(kind=R_GRID) ,
intent(IN) :: p3(ndims)
8691 real(kind=R_GRID) ,
intent(IN) :: p4(ndims)
8692 real(kind=R_GRID) ,
intent(IN) :: t1(ndims)
8693 integer,
intent(IN) :: dir
8695 real(kind=R_GRID) :: ee1(3), ee2(3), ee3(3), elon(3), elat(3)
8697 real :: g11, g12, g21, g22
8703 elon(1) = -sin(t1(1) -
pi)
8704 elon(2) = cos(t1(1) -
pi)
8706 elat(1) = -sin(t1(2))*cos(t1(1) -
pi)
8707 elat(2) = -sin(t1(2))*sin(t1(1) -
pi)
8708 elat(3) = cos(t1(2))
8716 newu = myu*g11 + myv*g12
8717 newv = myu*g21 + myv*g22
8719 newu = ( myu*g22 - myv*g12)/(g11*g22 - g21*g12)
8720 newv = (-myu*g21 + myv*g11)/(g11*g22 - g21*g12)
8729 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1)
8730 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed )
8731 integer,
intent(IN) :: npx, npy
8732 type(domain2d),
intent(INOUT) :: domain
8747 real ,
intent(INOUT) :: u(isd:ied ,jsd:jed+1,npz)
8748 real ,
intent(INOUT) :: v(isd:ied+1,jsd:jed ,npz)
8749 integer,
intent(IN) :: npx, npy, npz
8750 type(domain2d),
intent(INOUT) :: domain
8765 real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
result (gsum)
8767 integer,
intent(IN) :: npx, npy
8768 integer,
intent(IN) :: ifirst, ilast
8769 integer,
intent(IN) :: jfirst, jlast
8770 integer,
intent(IN) :: isd, ied
8771 integer,
intent(IN) :: jsd, jed, tile
8772 real ,
intent(IN) :: p(ifirst:ilast,jfirst:jlast)
8778 real,
allocatable :: p_r8(:,:,:)
8780 real,
pointer,
dimension(:,:,:) ::
agrid, grid
8781 real,
pointer,
dimension(:,:) :: area, rarea, fc,
f0 8782 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
8784 logical,
pointer :: cubed_sphere, latlon
8786 logical,
pointer :: have_south_pole, have_north_pole
8788 integer,
pointer :: ntiles_g
8789 real,
pointer :: acapn, acaps, globalarea
8791 grid => gridstruct%grid
8792 agrid=> gridstruct%agrid
8794 area => gridstruct%area
8795 rarea => gridstruct%rarea
8802 dxa => gridstruct%dxa
8803 dya => gridstruct%dya
8804 rdxa => gridstruct%rdxa
8805 rdya => gridstruct%rdya
8806 dxc => gridstruct%dxc
8807 dyc => gridstruct%dyc
8809 cubed_sphere => gridstruct%cubed_sphere
8810 latlon => gridstruct%latlon
8812 have_south_pole => gridstruct%have_south_pole
8813 have_north_pole => gridstruct%have_north_pole
8815 ntiles_g => gridstruct%ntiles_g
8816 acapn => gridstruct%acapN
8817 acaps => gridstruct%acapS
8818 globalarea => gridstruct%globalarea
8820 allocate(p_r8(npx-1,npy-1,ntiles_g))
8827 gsum = gsum + p(1,1)*acaps
8828 gsum = gsum + p(1,npy-1)*acapn
8831 gsum = gsum + p(i,j)*cos(
agrid(i,j,2))
8839 p_r8(i,j,n) = p(i,j)*area(i,j)
8843 call mp_gather(p_r8, ifirst,ilast, jfirst,jlast, npx-1, npy-1, ntiles_g)
8844 if (is_master())
then 8848 gsum = gsum + p_r8(i,j,n)
8852 gsum = gsum/globalarea
8864 real(kind=R_GRID),
intent(in):: p1(2), p2(2), p3(2)
8865 real(kind=R_GRID),
intent(out):: uvect(3)
8868 real(kind=R_GRID) :: xyz1(3), xyz2(3), xyz3(3)
8875 uvect(n) = xyz3(n)-xyz1(n)
8888 integer,
intent(in):: np
8889 real(kind=R_GRID),
intent(inout):: e(3,np)
8895 pdot = sqrt(e(1,n)**2+e(2,n)**2+e(3,n)**2)
8897 e(k,n) = e(k,n) / pdot
8907 subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
8908 kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
8911 integer,
intent(in):: im, jm, km, nq
8912 integer,
intent(in):: ifirst, ilast
8913 integer,
intent(in):: jfirst, jlast
8914 integer,
intent(in):: kfirst, klast
8915 integer,
intent(in):: ng_e
8916 integer,
intent(in):: ng_w
8917 integer,
intent(in):: ng_s
8918 integer,
intent(in):: ng_n
8919 real,
intent(inout):: q_ghst(ifirst-ng_w:ilast+ng_e,jfirst-ng_s:jlast+ng_n,kfirst:klast,nq)
8920 real,
optional,
intent(in):: q(ifirst:ilast,jfirst:jlast,kfirst:klast,nq)
8934 if (
present(q))
then 8935 q_ghst(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq) = &
8936 q(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq)
8942 do j=jfirst-ng_s,jlast+ng_n
8944 q_ghst(ifirst-i,j,k,n) = q_ghst(ilast-i+1,j,k,n)
8947 q_ghst(ilast+i,j,k,n) = q_ghst(ifirst+i-1,j,k,n)
8970 integer,
intent(in):: ifirst,ilast
8971 real,
intent(out) :: qout(ifirst:)
8972 real,
intent(in) :: qin(ifirst:)
8973 real,
intent(in) :: dx(ifirst:)
8974 integer,
intent(in):: order
8977 real :: dm(ifirst:ilast),qmax,qmin
8978 real :: r3, da1, da2, a6da, a6, al, ar
8979 real :: qLa, qLb1, qLb2
8988 qout(i) = 0.5 * (qin(i-1) + qin(i))
8990 elseif (order==2)
then 8993 qout(i) = (dx(i-1)*qin(i-1) + dx(i)*qin(i))/(dx(i-1)+dx(i))
8995 elseif (order==3)
then 8998 do i=ifirst+1,ilast-1
8999 dm(i) = 0.25*(qin(i+1) - qin(i-1))
9004 do i=ifirst+1,ilast-1
9005 qmax =
max(qin(i-1),qin(i),qin(i+1)) - qin(i)
9006 qmin = qin(i) -
min(qin(i-1),qin(i),qin(i+1))
9007 dm(i) = sign(
min(abs(dm(i)),qmin,qmax),dm(i))
9010 do i=ifirst+1,ilast-1
9011 qout(i) = 0.5*(qin(i-1)+qin(i)) + r3*(dm(i-1) - dm(i))
9018 qout(ifirst+1) = 0.5 * (qin(ifirst) + qin(ifirst+1))
9019 qout(ilast) = 0.5 * (qin(ilast-1) + qin(ilast))
9021 elseif (order==4)
then 9024 do i=ifirst+1,ilast-1
9025 dm(i) = ( (2.*dx(i-1) + dx(i) ) / &
9026 ( dx(i+1) + dx(i) ) ) * ( qin(i+1) - qin(i) ) + &
9027 ( (dx(i) + 2.*dx(i+1)) / &
9028 (dx(i-1) + dx(i) ) ) * ( qin(i) - qin(i-1) )
9029 dm(i) = ( dx(i) / ( dx(i-1) + dx(i) + dx(i+1) ) ) * dm(i)
9030 if ( (qin(i+1)-qin(i))*(qin(i)-qin(i-1)) > 0.)
then 9031 dm(i) = sign(
min( abs(dm(i)), 2.*abs(qin(i)-qin(i-1)), 2.*abs(qin(i+1)-qin(i)) ) , dm(i) )
9037 do i=ifirst+2,ilast-1
9038 qla = ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) - &
9039 ( (dx(i+1) + dx(i)) / (2.*dx(i) + dx(i-1)) )
9040 qla = ( (2.*dx(i) * dx(i-1)) / (dx(i-1) + dx(i)) ) * qla * &
9042 qlb1 = dx(i-1) * ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) * &
9044 qlb2 = dx(i) * ( (dx(i) + dx(i+1)) / (dx(i-1) + 2.*dx(i)) ) * &
9047 qout(i) = 1. / ( dx(i-2) + dx(i-1) + dx(i) + dx(i+1) )
9048 qout(i) = qout(i) * ( qla - qlb1 + qlb2 )
9049 qout(i) = qin(i-1) + ( dx(i-1) / ( dx(i-1) + dx(i) ) ) * (qin(i) - qin(i-1)) + qout(i)
9052 elseif (order==5)
then 9055 do i=ifirst+1,ilast-1
9056 x = float(i-(ifirst+1))*float(ilast-ifirst+1-1)/float(ilast-ifirst-1)
9057 qout(i) = qin(ifirst+nint(x)) + (x - nint(x)) * (qin(ifirst+nint(x+1)) - qin(ifirst+nint(x)))
9078 subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, &
9079 ng_d, ng_s, jfirst, jlast)
9086 integer,
intent(in):: ng_s, ng_d
9087 real,
intent(in):: coslon(im,jm), sinlon(im,jm)
9088 real,
intent(in):: cosl5(im,jm),sinl5(im,jm)
9089 real,
intent(in):: u(im,jfirst-ng_d:jlast+ng_s)
9092 real,
intent(inout):: v(im,jfirst-ng_d:jlast+ng_d)
9109 real uanp(im), uasp(im), vanp(im), vasp(im)
9110 real un, vn, us, vs, r2im
9113 r2im = 0.5d0/dble(im)
9118 if ( jfirst-ng_d <= 1 )
then 9120 uasp(i) = u(i, 2) + u(i,3)
9124 vasp(i) = v(i, 2) + v(i+1,2)
9126 vasp(im) = v(im,2) + v(1,2)
9132 us = us + (uasp(i+imh)-uasp(i))*sinlon(i,1) &
9133 + (vasp(i)-vasp(i+imh))*coslon(i,1)
9134 vs = vs + (uasp(i+imh)-uasp(i))*coslon(i,1) &
9135 + (vasp(i+imh)-vasp(i))*sinlon(i,1)
9143 v(i, 1) = us*cosl5(i,1) - vs*sinl5(i,1)
9144 v(i+imh,1) = -v(i,1)
9149 if ( jlast+ng_d >= jm )
then 9152 uanp(i) = u(i,jm-1) + u(i,jm)
9156 vanp(i) = v(i,jm-1) + v(i+1,jm-1)
9158 vanp(im) = v(im,jm-1) + v(1,jm-1)
9165 un = un + (uanp(i+imh)-uanp(i))*sinlon(i,jm) &
9166 + (vanp(i+imh)-vanp(i))*coslon(i,jm)
9167 vn = vn + (uanp(i)-uanp(i+imh))*coslon(i,jm) &
9168 + (vanp(i+imh)-vanp(i))*sinlon(i,jm)
9176 v(i, jm) = -un*cosl5(i,jm) - vn*sinl5(i,jm)
9177 v(i+imh,jm) = -v(i,jm)
9182 end subroutine vpol5 9184 subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
9186 character(len=*),
intent(in):: qname
9187 integer,
intent(in):: is, ie, js, je
9188 integer,
intent(in):: n_g, km
9189 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
9190 real,
intent(in):: fac
9201 if( q(i,j,k) < qmin )
then 9203 elseif( q(i,j,k) > qmax )
then 9210 write(*,*) qname,
' max = ', qmax*fac,
' min = ', qmin*fac
9214 subroutine var_dz(km, ztop, ze)
9215 integer,
intent(in):: km
9216 real,
intent(in):: ztop
9217 real,
intent(out),
dimension(km+1):: ze
9219 real,
dimension(km):: dz, s_fac
9230 s_fac(k) = 1.05 * s_fac(k+1)
9232 s_fac(4) = 1.1*s_fac(5)
9233 s_fac(3) = 1.2*s_fac(4)
9234 s_fac(2) = 1.3*s_fac(3)
9235 s_fac(1) = 1.5*s_fac(2)
9239 sum1 = sum1 + s_fac(k)
9245 dz(k) = s_fac(k) * dz0
9250 ze(k) = ze(k+1) + dz(k)
9255 dz(k) = dz(k) * (ztop/ze(1))
9259 ze(k) = ze(k+1) + dz(k)
9263 call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
9265 if ( is_master() )
then 9266 write(*,*)
'var_dz: model top (km)=', ztop*0.001
9268 dz(k) = ze(k) - ze(k+1)
9269 write(*,*) k, 0.5*(ze(k)+ze(k+1)),
'dz=', dz(k)
9275 subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
9276 integer,
intent(in):: is, ie, js, je, km
9277 integer,
intent(in):: ntimes, i, j
9278 real,
intent(inout):: ze(is:ie,js:je,km+1)
9280 real,
parameter:: df = 0.25
9283 integer k, n, k1, k2
9287 dz(k) = ze(i,j,k+1) - ze(i,j,k)
9296 flux(k) = df*(dz(k) - dz(k-1))
9300 dz(k) = dz(k) - flux(k) + flux(k+1)
9305 ze(i,j,k) = ze(i,j,k+1) - dz(k)
subroutine, public eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, hydrostatic, moist)
real, parameter, public radius
Radius of the Earth [m].
integer, parameter, public model_atmos
subroutine mp_update_dwinds_3d(u, v, npx, npy, npz, domain)
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)
integer, public tracer_test
real, parameter, public omega
Rotation rate of the Earth [1/s].
subroutine, public get_unit_vect2(e1, e2, uc)
integer, parameter initwindscase6
real, parameter, public ptop_min
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)
integer, parameter, public dgrid_ne
integer, parameter initwindscase2
real(kind=8), parameter, public pi_8
Ratio of circle circumference to diameter [N/A].
subroutine, public case9_forcing1(phis, time_since_start)
subroutine var_dz(km, ztop, ze)
real(kind=kind_real), parameter f0
Coriolis parameter at southern boundary.
real function dcmip16_tc_pressure(z, r)
subroutine superk_u(km, zz, um, dudz)
integer, parameter initwindscase5
real, dimension(:,:,:), allocatable ua0
real, dimension(:), allocatable gh_table
integer, parameter, public up
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)
subroutine mp_update_dwinds_2d(u, v, npx, npy, domain)
subroutine, public case9_forcing2(phis)
real function dcmip16_bc_uwind_pert(z, lat, lon)
real function dcmip16_bc_temperature(z, lat)
real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
subroutine init_winds(UBar, u, v, ua, va, uc, vc, defOnGrid, npx, npy, ng, ndims, nregions, nested, gridstruct, domain, tile)
void mid_pt_sphere(const double *p1, const double *p2, double *pm)
real, parameter, public rdgas
Gas constant for dry air [J/kg/deg].
subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng)
subroutine get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
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 gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
subroutine superk_sounding(km, pe, p00, ze, pt, qz)
subroutine dcmip16_tc_uwind_pert(z, r, lon, lat, uu, vv)
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
subroutine get_pt_on_great_circle(p1, p2, dist, heading, p3)
subroutine init_latlon_winds(UBar, u, v, ua, va, uc, vc, defOnGrid, gridstruct)
subroutine, public project_sphere_v(np, f, e)
subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
real, dimension(:,:,:), allocatable va0
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 pmxn(p, npx, npy, nregions, tile, gridstruct, pmin, pmax, i0, j0, n0)
subroutine get_case9_b(B, agrid)
integer, parameter, public agrid
subroutine, public hybrid_z_dz(km, dz, ztop, s_rate)
subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, ng_d, ng_s, jfirst, jlast)
subroutine interp_left_edge_1d(qout, qin, dx, ifirst, ilast, order)
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)
integer, parameter initwindscase9
real, dimension(:), allocatable, public zz0
integer, parameter, public f_p
real function gh_jet(npy, lat_in)
subroutine, public cart_to_latlon(np, q, xs, ys)
real, parameter, public rvgas
Gas constant for water vapor [J/kg/deg].
subroutine, public check_courant_numbers(uc, vc, ndt, n_split, gridstruct, npx, npy, npz, tile, noPrint)
real(kind=kind_real), parameter u1
subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, nested, domain, noComm)
subroutine get_vector_stats(varU, varUT, varV, varVT, npx, npy, ndims, nregions, vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
integer, parameter initwindscase0
integer, public wind_field
integer, parameter, public ng
subroutine, public hydro_eq(km, is, ie, js, je, ps, hs, drym, delp, ak, bk, pt, delz, area, ng, mountain, hydrostatic, hybrid_z, domain)
real, public soliton_size
real function dcmip16_bc_sphum(p, ps, lat, lon)
subroutine d2a2c(im, jm, km, ifirst, ilast, jfirst, jlast, ng, nested, u, v, ua, va, uc, vc, gridstruct, domain)
real(fp), parameter, public one
real, parameter, public cp_air
Specific heat capacity of dry air at constant pressure [J/kg/deg].
integer, parameter interporder
integer, parameter, public r_grid
subroutine rankine_vortex(ubar, r0, p1, u, v, grid)
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)
real function, public great_circle_dist(q1, q2, radius)
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng)
integer, public nsolitons
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(z, T, lat)
real, public soliton_umax
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)
real, dimension(:,:), allocatable case9_b
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)
real, dimension(:,:,:), allocatable phi0
real, parameter, public grav
Acceleration due to gravity [m/s^2].
subroutine get_scalar_stats(var, varT, npx, npy, ndims, nregions, vmin, vmax, L1_norm, L2_norm, Linf_norm, gridstruct, tile)
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
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)
subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, nested, domain)
subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, delz, zvir, ptop, ak, bk, agrid)
real function dcmip16_bc_pressure(z, lat)
void normalize_vect(double *e)
logical, public bubble_do
real, dimension(:), allocatable lats_table
subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
integer, parameter initwindscase1
subroutine terminator_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, km, q, delp, ncnst, lon, lat)
subroutine, public set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
subroutine, public compute_dz_l101(km, ztop, dz)
real function dcmip16_tc_temperature(z, r)
subroutine, public compute_dz_l32(km, ztop, dz)
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
integer, public test_case
real function, public inner_prod(v1, v2)
subroutine, public ppme(p, qe, delp, im, km)
real, dimension(:), allocatable, public pz0
void latlon2xyz(int size, const double *lon, const double *lat, double *x, double *y, double *z)
real function dcmip16_tc_sphum(z)
integer, parameter, public scalar_pair
real, parameter, public kappa
RDGAS / CP_AIR [dimensionless].
integer, parameter, public cgrid_ne
subroutine, public case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain)
subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
real(fp), parameter, public pi