100 use constants_mod
, only: cnst_radius=>
radius, pi=>pi_8, omega, grav, kappa, rdgas, cp_air, rvgas
103 domain_decomp, fill_corners, xdir, ydir, &
104 mp_stop, mp_reduce_sum, mp_reduce_max, mp_gather, mp_bcst
114 use mpp_mod
, only: mpp_error, fatal, mpp_root_pe, mpp_broadcast, mpp_sum
115 use mpp_mod
, only: stdlog, input_nml_file
116 use fms_mod
, only: check_nml_error, close_file, open_namelist_file
117 use mpp_domains_mod
, only: mpp_update_domains, domain2d
118 use mpp_parameter_mod
, only: agrid_param=>agrid,cgrid_ne_param=>cgrid_ne, &
123 use mpp_mod
, only: mpp_pe, mpp_chksum, stdout
126 use tracer_manager_mod
, only: get_tracer_index
127 use field_manager_mod
, only: model_atmos
178 real(kind=R_GRID),
parameter ::
radius = cnst_radius
179 real(kind=R_GRID),
parameter ::
one = 1.d0
203 real,
allocatable,
dimension(:) ::
pz0,
zz0 218 real ,
allocatable ::
phi0(:,:,:)
219 real ,
allocatable ::
ua0(:,:,:)
220 real ,
allocatable ::
va0(:,:,:)
236 public :: output, output_ncdf
254 subroutine init_winds(UBar, u,v,ua,va,uc,vc, defOnGrid, npx, npy, ng, ndims, nregions, bounded_domain, gridstruct, domain, tile, bd)
258 real ,
intent(INOUT) :: UBar
259 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
260 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed )
261 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed )
262 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
263 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed )
264 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed )
265 integer,
intent(IN) :: defOnGrid
266 integer,
intent(IN) :: npx, npy
267 integer,
intent(IN) :: ng
268 integer,
intent(IN) :: ndims
269 integer,
intent(IN) :: nregions
270 logical,
intent(IN) :: bounded_domain
272 type(domain2d),
intent(INOUT) :: domain
273 integer,
intent(IN) :: tile
275 real(kind=R_GRID) :: p1(2), p2(2), p3(2), p4(2), pt(2)
276 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3)
282 real :: psi_b(bd%isd:bd%ied+1,bd%jsd:bd%jed+1), psi(bd%isd:bd%ied,bd%jsd:bd%jed), psi1, psi2
283 integer :: is2, ie2, js2, je2
285 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
286 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
287 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
288 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
289 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
291 logical,
pointer :: cubed_sphere, latlon
293 logical,
pointer :: have_south_pole, have_north_pole
295 integer,
pointer :: ntiles_g
296 real,
pointer :: acapN, acapS, globalarea
298 integer :: is, ie, js, je
299 integer :: isd, ied, jsd, jed
301 grid => gridstruct%grid_64
302 agrid=> gridstruct%agrid_64
304 area => gridstruct%area
305 rarea => gridstruct%rarea
310 ee1 => gridstruct%ee1
311 ee2 => gridstruct%ee2
314 en1 => gridstruct%en1
315 en2 => gridstruct%en2
319 dxa => gridstruct%dxa
320 dya => gridstruct%dya
321 rdxa => gridstruct%rdxa
322 rdya => gridstruct%rdya
323 dxc => gridstruct%dxc
324 dyc => gridstruct%dyc
326 cubed_sphere => gridstruct%cubed_sphere
327 latlon => gridstruct%latlon
329 have_south_pole => gridstruct%have_south_pole
330 have_north_pole => gridstruct%have_north_pole
332 ntiles_g => gridstruct%ntiles_g
333 acapn => gridstruct%acapN
334 acaps => gridstruct%acapS
335 globalarea => gridstruct%globalarea
346 if (bounded_domain)
then 362 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)
368 psi(i,j) = (-1.0 * ubar *
radius *( sin(agrid(i,j,2)) *cos(
alpha) - &
369 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
372 call mpp_update_domains( psi, domain )
375 psi_b(i,j) = (-1.0 * ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
376 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
380 if ( (cubed_sphere) .and. (defongrid==0) )
then 384 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
385 if (dist==0) vc(i,j) = 0.
391 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
392 if (dist==0) uc(i,j) = 0.
395 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
396 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
400 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
401 if (dist==0) v(i,j) = 0.
407 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
408 if (dist==0) u(i,j) = 0.
414 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
415 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
417 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
418 if (dist==0) ua(i,j) = 0.
419 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
420 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
422 va(i,j) = (psi2 - psi1) / (dist)
423 if (dist==0) va(i,j) = 0.
427 elseif ( (cubed_sphere) .and. (defongrid==1) )
then 431 vc(i,j) = (psi_b(i+1,j)-psi_b(i,j))/dist
432 if (dist==0) vc(i,j) = 0.
438 uc(i,j) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
439 if (dist==0) uc(i,j) = 0.
442 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
443 call fill_corners(uc, vc, npx, npy, vector=.true., cgrid=.true.)
444 call ctoa(uc,vc,ua,va,dx, dy, dxc,dyc,dxa,dya,npx,npy,ng, bd)
445 call atod(ua,va,u ,v ,dxa, dya,dxc,dyc,npx,npy,ng, bounded_domain, domain, bd)
448 elseif ( (cubed_sphere) .and. (defongrid==2) )
then 452 v(i,j) = (psi(i,j)-psi(i-1,j))/dist
453 if (dist==0) v(i,j) = 0.
459 u(i,j) = -1.0*(psi(i,j)-psi(i,j-1))/dist
460 if (dist==0) u(i,j) = 0.
464 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng, bd)
465 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, bounded_domain, domain, bd)
466 elseif ( (cubed_sphere) .and. (defongrid==3) )
then 469 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
470 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
472 ua(i,j) = -1.0 * (psi2 - psi1) / (dist)
473 if (dist==0) ua(i,j) = 0.
474 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
475 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
477 va(i,j) = (psi2 - psi1) / (dist)
478 if (dist==0) va(i,j) = 0.
481 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
482 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, bounded_domain, domain, bd)
483 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, bounded_domain,domain, bd)
484 elseif ( (latlon) .or. (defongrid==4) )
then 488 ua(i,j) = ubar * ( cos(agrid(i,j,2))*cos(
alpha) + &
489 sin(agrid(i,j,2))*cos(agrid(i,j,1))*sin(
alpha) )
490 va(i,j) = -ubar * sin(agrid(i,j,1))*sin(
alpha)
495 if (cubed_sphere)
call rotate_winds(ua(i,j), va(i,j), p1,p2,p3,p4, agrid(i,j,1:2), 2, 1)
497 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
498 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
500 if ( (tile==1) .and.(i==1) ) print*, ua(i,j), -1.0 * (psi2 - psi1) / (dist)
504 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
505 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, bounded_domain, domain, bd)
506 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, bounded_domain, domain, bd)
507 elseif ( (latlon) .or. (defongrid==5) )
then 512 p1(:) = grid(i ,j ,1:2)
513 p2(:) = grid(i,j+1 ,1:2)
517 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
518 sin(pt(2))*cos(pt(1))*sin(
alpha) )
519 vtmp = -ubar * sin(pt(1))*sin(
alpha)
526 p1(:) = grid(i ,j ,1:2)
527 p2(:) = grid(i+1,j ,1:2)
531 utmp = ubar * ( cos(pt(2))*cos(
alpha) + &
532 sin(pt(2))*cos(pt(1))*sin(
alpha) )
533 vtmp = -ubar * sin(pt(1))*sin(
alpha)
539 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng, bd)
540 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, bounded_domain, domain, bd)
561 subroutine init_case(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
562 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, &
563 dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, adiabatic, &
564 ks, npx_global, ptop, domain_in, tile_in, bd)
567 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
568 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
569 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
570 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
571 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
572 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
574 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
576 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
577 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
578 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
579 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
580 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
582 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
583 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
584 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
585 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
586 real ,
intent(inout) :: delz(bd%is:,bd%js:,1:)
587 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
589 real ,
intent(inout) :: ak(npz+1)
590 real ,
intent(inout) :: bk(npz+1)
592 integer,
intent(IN) :: npx, npy, npz
593 integer,
intent(IN) :: ng, ncnst, nwat
594 integer,
intent(IN) :: ndims
595 integer,
intent(IN) :: nregions
597 real,
intent(IN) :: dry_mass
598 logical,
intent(IN) :: mountain
599 logical,
intent(IN) :: moist_phys
600 logical,
intent(IN) :: hydrostatic
601 logical,
intent(IN) :: hybrid_z
602 logical,
intent(IN) :: adiabatic
603 integer,
intent(IN) :: ks
608 integer,
intent(IN) :: npx_global
609 integer,
intent(IN),
target :: tile_in
610 real,
intent(INOUT) :: ptop
612 type(domain2d),
intent(IN),
target :: domain_in
614 real :: tmp(1-ng:npx +ng,1-ng:npy +ng,1:nregions)
615 real :: tmp1(1 :npx ,1 :npy ,1:nregions)
617 real(kind=R_GRID) :: p0(2)
618 real(kind=R_GRID) :: p1(2)
619 real(kind=R_GRID) :: p2(2)
620 real(kind=R_GRID) :: p3(2)
621 real(kind=R_GRID) :: p4(2)
622 real(kind=R_GRID) :: pa(2)
623 real(kind=R_GRID) :: pb(2)
624 real(kind=R_GRID) :: pcen(2)
625 real(kind=R_GRID) :: e1(3), e2(3), e3(3), ex(3), ey(3)
626 real :: dist, r, r1, r2, r0, omg, A, B, C
627 integer :: i,j,k,nreg,z,zz
628 integer :: i0,j0,n0, nt
629 real :: utmp,vtmp,ftmp
632 integer,
parameter :: jm = 5761
639 real :: ddeg, deg, DDP, DP, ph5
644 real :: x1,y1,z1,x2,y2,z2,ang
646 integer :: initWindsCase
658 real :: grad(bd%isd:bd%ied ,bd%jsd:bd%jed,2)
659 real :: div0(bd%isd:bd%ied ,bd%jsd:bd%jed )
660 real :: vor0(bd%isd:bd%ied ,bd%jsd:bd%jed )
661 real :: divg(bd%isd:bd%ied ,bd%jsd:bd%jed )
662 real :: vort(bd%isd:bd%ied ,bd%jsd:bd%jed )
663 real :: ztop, rgrav, p00, pturb, zmid, pk0, t00
664 real :: dz1(npz), ppt(npz)
665 real :: ze1(npz+1), pe1(npz+1)
668 character(len=80) :: oflnm, hgtflnm
669 integer :: is2, ie2, js2, je2
671 real :: psi(bd%isd:bd%ied,bd%jsd:bd%jed)
672 real :: psi_b(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
676 real :: eta(npz), eta_0, eta_s, eta_t
677 real :: eta_v(npz), press, anti_rot
678 real :: T_0, T_mean, delta_T, lapse_rate, n2, zeta, s0
679 real :: pt1,pt2,pt3,pt4,pt5,pt6, pt7, pt8, pt9, u1, pt0
680 real :: uu1, uu2, uu3, vv1, vv2, vv3
683 real wbuffer(npy+2,npz)
684 real sbuffer(npx+2,npz)
686 real :: gz(bd%isd:bd%ied,bd%jsd:bd%jed,npz+1), zt, zdist
693 real,
dimension(npz):: pk1, ts1, qs1, uz1, zs1, dudz
695 real(kind=R_GRID):: pp0(2)
700 real :: omega0, k_cell, z0, H, px
701 real :: d1, d2, p1p(2), rt, s
702 real :: wind_alpha, period, h0, rm, zp3(3), dz3(3), k0, lp
706 real,
dimension(npz+1) :: pe0, gz0, ue, ve, we, pte, qe
707 real :: d, cor, exppr, exppz, gamma, Ts0, q00, exponent, ztrop, height, zp, rp
708 real :: qtrop, ttrop, zq1, zq2
709 real :: dum, dum1, dum2, dum3, dum4, dum5, dum6, ptmp, uetmp, vetmp
710 real :: pe_u(bd%is:bd%ie,npz+1,bd%js:bd%je+1)
711 real :: pe_v(bd%is:bd%ie+1,npz+1,bd%js:bd%je)
712 real :: ps_u(bd%is:bd%ie,bd%js:bd%je+1)
713 real :: ps_v(bd%is:bd%ie+1,bd%js:bd%je)
718 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
719 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
720 real,
pointer,
dimension(:,:) :: rarea, fC, f0
721 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
722 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
723 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
725 logical,
pointer :: cubed_sphere, latlon
727 type(domain2d),
pointer :: domain
728 integer,
pointer :: tile
730 logical,
pointer :: have_south_pole, have_north_pole
732 integer,
pointer :: ntiles_g
733 real,
pointer :: acapN, acapS, globalarea
735 integer :: is, ie, js, je
736 integer :: isd, ied, jsd, jed
747 grid => gridstruct%grid_64
748 agrid=> gridstruct%agrid_64
750 area => gridstruct%area_64
751 rarea => gridstruct%rarea
756 ee1 => gridstruct%ee1
757 ee2 => gridstruct%ee2
760 en1 => gridstruct%en1
761 en2 => gridstruct%en2
765 dxa => gridstruct%dxa
766 dya => gridstruct%dya
767 rdxa => gridstruct%rdxa
768 rdya => gridstruct%rdya
769 dxc => gridstruct%dxc
770 dyc => gridstruct%dyc
772 cubed_sphere => gridstruct%cubed_sphere
773 latlon => gridstruct%latlon
778 have_south_pole => gridstruct%have_south_pole
779 have_north_pole => gridstruct%have_north_pole
781 ntiles_g => gridstruct%ntiles_g
782 acapn => gridstruct%acapN
783 acaps => gridstruct%acapS
784 globalarea => gridstruct%globalarea
786 if (gridstruct%bounded_domain)
then 800 f0(:,:) = huge(dummy)
801 fc(:,:) = huge(dummy)
804 fc(i,j) = 2.*omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
805 sin(grid(i,j,2))*cos(
alpha) )
810 f0(i,j) = 2.*omega*( -1.*cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) + &
811 sin(agrid(i,j,2))*cos(
alpha) )
814 call mpp_update_domains( f0, domain )
815 if (cubed_sphere)
call fill_corners(f0, npx, npy, ydir)
817 delp(isd:is-1,jsd:js-1,1:npz)=0.
818 delp(isd:is-1,je+1:jed,1:npz)=0.
819 delp(ie+1:ied,jsd:js-1,1:npz)=0.
820 delp(ie+1:ied,je+1:jed,1:npz)=0.
822 #if defined(SW_DYNAMICS) 832 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
833 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
836 call init_winds(
ubar, u,v,ua,va,uc,vc, 1, npx, npy, ng, ndims, nregions, gridstruct%bounded_domain, gridstruct, domain, tile)
841 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
842 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
843 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)
849 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
850 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
858 call get_scalar_stats( divg, div0, npx, npy, ndims, nregions, &
859 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
860 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)
861 201
format(
' ',a,e21.14,
' ',e21.14)
862 202
format(
' ',a,i4.4,
'x',i4.4,
'x',i4.4)
863 if ( is_master() )
then 864 write(*,*)
' Error Norms of Analytical Divergence field C-Winds initialized' 865 write(*,201)
'Divergence MAX error : ', pmax
866 write(*,201)
'Divergence MIN error : ', pmin
867 write(*,201)
'Divergence L1_norm : ', l1_norm
868 write(*,201)
'Divergence L2_norm : ', l2_norm
869 write(*,201)
'Divergence Linf_norm : ', linf_norm
872 call init_winds(
ubar, u,v,ua,va,uc,vc, 3, npx, npy, ng, ndims, nregions, gridstruct%bounded_domain, gridstruct, domain, tile, bd)
876 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
877 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
878 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)
884 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
885 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
891 call get_scalar_stats( divg, div0, npx, npy, ndims, nregions, &
892 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
893 if ( is_master() )
then 894 write(*,*)
' Error Norms of Analytical Divergence field A-Winds initialized' 895 write(*,201)
'Divergence MAX error : ', pmax
896 write(*,201)
'Divergence MIN error : ', pmin
897 write(*,201)
'Divergence L1_norm : ', l1_norm
898 write(*,201)
'Divergence L2_norm : ', l2_norm
899 write(*,201)
'Divergence Linf_norm : ', linf_norm
902 call init_winds(
ubar, u,v,ua,va,uc,vc, 2, npx, npy, ng, ndims, nregions, gridstruct%bounded_domain, gridstruct, domain, tile, bd)
908 divg(i,j) = (rarea(i,j)) * ( (uc(i+1,j,1)*dy(i+1,j) - uc(i,j,1)*dy(i,j)) + &
909 (vc(i,j+1,1)*dx(i,j+1) - vc(i,j,1)*dx(i,j)) )
910 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)
916 vort(i,j) = (rarea(i,j)) * ( (v(i+1,j,1)*dy(i+1,j) - v(i,j,1)*dy(i,j)) - &
917 (u(i,j+1,1)*dx(i,j+1) - u(i,j,1)*dx(i,j)) )
921 call get_scalar_stats( divg, div0, npx, npy, ndims, nregions, &
922 pmin, pmax, l1_norm, l2_norm, linf_norm, gridstruct, tile)
923 if ( is_master() )
then 924 write(*,*)
' Error Norms of Analytical Divergence field D-Winds initialized' 925 write(*,201)
'Divergence MAX error : ', pmax
926 write(*,201)
'Divergence MIN error : ', pmin
927 write(*,201)
'Divergence L1_norm : ', l1_norm
928 write(*,201)
'Divergence L2_norm : ', l2_norm
929 write(*,201)
'Divergence Linf_norm : ', linf_norm
943 vtx = ((3.0*sqrt(2.0))/2.0) * (( 1.0/cosh(p) )**2.0) * tanh(p)
945 if (p /= 0.0) w_p = vtx/p
946 delp(i,j,1) = 1.0 - tanh( (p/
rgamma) * sin(x1 - w_p*0.0) )
947 ua(i,j,1) = w_p*(sin(
lat0)*cos(agrid(i,j,2)) + cos(
lat0)*cos(agrid(i,j,1) -
lon0)*sin(agrid(i,j,2)))
948 va(i,j,1) = w_p*cos(
lat0)*sin(agrid(i,j,1) -
lon0)
949 ua(i,j,1) = ua(i,j,1)*
radius/86400.0
950 va(i,j,1) = va(i,j,1)*
radius/86400.0
956 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)
960 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
961 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
963 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
964 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
965 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
980 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
982 delp(i,j,1) = phis(i,j)
1000 p2(1) = agrid(i,j,1)
1001 p2(2) = agrid(i,j,2)
1003 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.))
then 1021 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
1022 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
1024 ( -1.*cos(grid(i ,j ,1))*cos(grid(i ,j ,2))*sin(
alpha) + &
1025 sin(grid(i ,j ,2))*cos(
alpha) ) ** 2.0
1027 ( -1.*cos(grid(i+1,j ,1))*cos(grid(i+1,j ,2))*sin(
alpha) + &
1028 sin(grid(i+1,j ,2))*cos(
alpha) ) ** 2.0
1030 ( -1.*cos(grid(i+1,j+1,1))*cos(grid(i+1,j+1,2))*sin(
alpha) + &
1031 sin(grid(i+1,j+1,2))*cos(
alpha) ) ** 2.0
1033 ( -1.*cos(grid(i,j+1,1))*cos(grid(i,j+1,2))*sin(
alpha) + &
1034 sin(grid(i,j+1,2))*cos(
alpha) ) ** 2.0
1035 delp(i,j,1) = (0.25*(pt1+pt2+pt3+pt4) + 3.*pt5) / 4.
1038 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
1039 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2.0
1060 p2(1) = agrid(i,j,1)
1061 p2(2) = agrid(i,j,2)
1064 delp(i,j,1) = phis(i,j) +
gh0*0.5*(1.0+cos(pi*r/r0))
1066 delp(i,j,1) = phis(i,j)
1069 delp(i,j,1) = delp(i,j,1) + grav*2.e3
1080 p1(:) = grid(i ,j ,1:2)
1081 p2(:) = grid(i,j+1 ,1:2)
1085 utmp =
ubar * cos(p3(2))
1092 p1(:) = grid(i, j,1:2)
1093 p2(:) = grid(i+1,j,1:2)
1097 utmp =
ubar * cos(p3(2))
1106 fc(i,j) = 2.*anti_rot*sin(grid(i,j,2))
1111 f0(i,j) = 2.*anti_rot*sin(agrid(i,j,2))
1140 p1(1) = pi*1.5 - ddeg
1144 p2(1) = pi*1.5 + ddeg
1148 #ifndef SINGULAR_VORTEX 1180 p2(1) = agrid(i,j,1)
1181 p2(2) = agrid(i,j,2)
1182 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1184 phis(i,j) = 2000.0*grav*(1.0-(r/r0))
1190 ( -1.*cos(agrid(i ,j ,1))*cos(agrid(i ,j ,2))*sin(
alpha) + &
1191 sin(agrid(i ,j ,2))*cos(
alpha) ) ** 2 - phis(i,j)
1203 a = 0.5*omg*(2.*omega+omg)*(cos(agrid(i,j,2))**2) + &
1204 0.25*rk*rk*(cos(agrid(i,j,2))**(r+r)) * &
1205 ( (r+1)*(cos(agrid(i,j,2))**2) + (2.*r*r-r-2.) - &
1206 2.*(r*r)*cos(agrid(i,j,2))**(-2.) )
1207 b = (2.*(omega+omg)*rk / ((r+1)*(r+2))) * (cos(agrid(i,j,2))**r) * &
1208 ( (r*r+2.*r+2.) - ((r+1.)*cos(agrid(i,j,2)))**2 )
1209 c = 0.25*rk*rk*(cos(agrid(i,j,2))**(2.*r)) * ( &
1210 (r+1) * (cos(agrid(i,j,2))**2.) - (r+2.) )
1211 delp(i,j,1) =
gh0 +
radius*
radius*(a+b*cos(r*agrid(i,j,1))+c*cos(2.*r*agrid(i,j,1)))
1212 delp(i,j,1) = delp(i,j,1) - phis(i,j)
1217 p1(:) = grid(i ,j ,1:2)
1218 p2(:) = grid(i,j+1 ,1:2)
1222 utmp =
radius*omg*cos(p3(2)) + &
1223 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1224 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1230 p1(:) = grid(i, j,1:2)
1231 p2(:) = grid(i+1,j,1:2)
1235 utmp =
radius*omg*cos(p3(2)) + &
1236 radius*rk*(cos(p3(2))**(r-1))*(r*sin(p3(2))**2-cos(p3(2))**2)*cos(r*p3(1))
1237 vtmp = -
radius*rk*r*sin(p3(2))*sin(r*p3(1))*cos(p3(2))**(r-1)
1242 call dtoa( u, v,ua,va,dx,dy,dxa,dya,dxc,dyc,npx,npy,ng,bd)
1244 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
1262 pt1 =
gh_jet(npy, agrid(i,j,2))
1271 pt6 =
gh_jet(npy, grid(i, j, 2))
1272 pt7 =
gh_jet(npy, grid(i+1,j, 2))
1273 pt8 =
gh_jet(npy, grid(i+1,j+1,2))
1274 pt9 =
gh_jet(npy, grid(i ,j+1,2))
1275 ftmp = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1277 delp(i,j,1) = ftmp + 120.*grav*cos(agrid(i,j,2)) * &
1278 exp( -(3.*(agrid(i,j,1)-pi))**2 ) * exp( -(15.*(agrid(i,j,2)-pi/4.))**2 )
1284 p1(:) = agrid(i,j,1:2)
1287 if ( r < 3.*r0 )
then 1288 delp(i,j,1) = delp(i,j,1) + 1000.*grav*exp(-(r/r0)**2)
1297 p2(:) = grid(i,j+1,1:2)
1298 vv1 =
u_jet(p2(2))*(ee2(2,i,j+1)*cos(p2(1)) - ee2(1,i,j+1)*sin(p2(1)))
1299 p1(:) = grid(i,j,1:2)
1300 vv3 =
u_jet(p1(2))*(ee2(2,i,j)*cos(p1(1)) - ee2(1,i,j)*sin(p1(1)))
1303 vv2 =
u_jet(pa(2))*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1305 v(i,j,1) = 0.25*(vv1 + 2.*vv2 + vv3)
1312 p1(:) = grid(i,j,1:2)
1313 uu1 =
u_jet(p1(2))*(ee1(2,i,j)*cos(p1(1)) - ee1(1,i,j)*sin(p1(1)))
1314 p2(:) = grid(i+1,j,1:2)
1315 uu3 =
u_jet(p2(2))*(ee1(2,i+1,j)*cos(p2(1)) - ee1(1,i+1,j)*sin(p2(1)))
1318 uu2 =
u_jet(pa(2))*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1320 u(i,j,1) = 0.25*(uu1 + 2.*uu2 + uu3)
1327 call get_vorticity(is, ie, js, je, isd, ied, jsd, jed, npz, u, v, q(is:ie,js:je,:,1), dx, dy, rarea)
1330 fc(i,j) = 2.*omega*( -1.*cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) + &
1331 sin(grid(i,j,2))*cos(
alpha) )
1336 f0(i,j) = 2.*omega*( -1.*cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) + &
1337 sin(agrid(i,j,2))*cos(
alpha) )
1340 call mpp_update_domains( f0, domain )
1341 if (cubed_sphere)
call fill_corners(f0, npx, npy, ydir)
1344 q(i,j,npz,1) = ( q(i,j,npz,1) + f0(i,j) ) / delp(i,j,npz) * 1.e6
1362 p2(1) = agrid(i,j,1)
1363 p2(2) = agrid(i,j,2)
1364 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1366 phis(i,j) = 2000.0*grav*(1.0-(r/r0))
1381 if ( is_master() )
write(*,*)
'Initialzing case-8: soliton twin cycolne...' 1404 p1(:) = grid(i ,j ,1:2)
1405 p2(:) = grid(i,j+1 ,1:2)
1408 utmp =
ubar*exp(-(r/r0)**2)
1416 p1(:) = grid(i, j,1:2)
1417 p2(:) = grid(i+1,j,1:2)
1420 utmp =
ubar*exp(-(r/r0)**2)
1433 p1(:) = grid(i ,j ,1:2)
1434 p2(:) = grid(i,j+1 ,1:2)
1437 utmp =
ubar*exp(-(r/r0)**2)
1445 p1(:) = grid(i, j,1:2)
1446 p2(:) = grid(i+1,j,1:2)
1449 utmp =
ubar*exp(-(r/r0)**2)
1464 ph5 = -0.5*pi + (dble(j-1)-0.5)*ddp
1465 ll_j(j) = -0.5*pi + (dble(j-1)*ddp)
1471 cosp(j) = (sine(j+1)-sine(j)) / dp
1474 cose(j) = 0.5 * (cosp(j-1) + cosp(j))
1477 ddeg = 180./float(jm-1)
1479 deg = -90. + (float(j-1)-0.5)*ddeg
1481 ll_u(j) = -10.*(deg+90.)/90.
1482 elseif (deg <= 60.)
then 1483 ll_u(j) = -10. + deg
1485 ll_u(j) = 50. - (50./30.)* (deg - 60.)
1488 ll_phi(1) = 6000. * grav
1490 ll_phi(j)=ll_phi(j-1) - dp*sine(j) * &
1491 (
radius*2.*omega + ll_u(j)/cose(j))*ll_u(j)
1497 if ( (ll_j(jj) <= agrid(i,j,2)) .and. (agrid(i,j,2) <= ll_j(jj+1)) )
then 1498 delp(i,j,1)=0.5*(ll_phi(jj)+ll_phi(jj+1))
1506 if (agrid(i,j,2)*
todeg <= 0.0)
then 1507 ua(i,j,1) = -10.*(agrid(i,j,2)*
todeg + 90.)/90.
1508 elseif (agrid(i,j,2)*
todeg <= 60.0)
then 1509 ua(i,j,1) = -10. + agrid(i,j,2)*
todeg 1511 ua(i,j,1) = 50. - (50./30.)* (agrid(i,j,2)*
todeg - 60.)
1518 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)
1522 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
1523 call atoc(ua,va,uc,vc,dx,dy,dxa,dya,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
1524 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
1525 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
1526 call atod(ua,va, u, v,dxa, dya,dxc,dyc,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
1537 if ( is_master() )
write(*,*)
'Initialzing case-9: soliton cyclones...' 1559 p1(:) = grid(i ,j ,1:2)
1560 p2(:) = grid(i,j+1 ,1:2)
1563 utmp =
ubar*exp(-(r/r0)**2)
1571 p1(:) = grid(i, j,1:2)
1572 p2(:) = grid(i+1,j,1:2)
1575 utmp =
ubar*exp(-(r/r0)**2)
1588 delp(:,:,z) = delp(:,:,1)
1591 call mpp_update_domains( delp, domain )
1592 call mpp_update_domains( phis, domain )
1595 call init_winds(
ubar, u,v,ua,va,uc,vc, initwindscase, npx, npy, ng, ndims, nregions, gridstruct%bounded_domain, gridstruct, domain, tile, bd)
1604 ps(i,j) = delp(i,j,1)
1618 if (.not.hydrostatic) w(:,:,:)= 0.0
1623 call hydro_eq(npz, is, ie, js, je, ps, phis, 1.e5, &
1624 delp, ak, bk, pt, delz, area, ng, .false., hydrostatic, hybrid_z, domain)
1631 p1(2) = pi/6. + (7.5/180.0)*pi
1634 p2(1) = agrid(i,j,1)
1635 p2(2) = agrid(i,j,2)
1636 r = min(r0*r0, (p2(1)-p1(1))*(p2(1)-p1(1)) + (p2(2)-p1(2))*(p2(2)-p1(2)) )
1638 phis(i,j) =
gh0*(1.0-(r/r0))
1641 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1642 delp, ak, bk, pt, delz, area, ng, mountain, hydrostatic, hybrid_z, domain)
1646 call surfdrv(npx, npy, gridstruct%grid_64, gridstruct%agrid_64, &
1647 gridstruct%area_64, dx, dy, dxa, dya, dxc, dyc, &
1648 gridstruct%sin_sg, phis, &
1649 flagstruct%stretch_fac, gridstruct%nested, gridstruct%bounded_domain, &
1650 npx_global, domain, flagstruct%grid_number, bd)
1651 call mpp_update_domains( phis, domain )
1653 if ( hybrid_z )
then 1659 if ( is_master() )
write(*,*)
'Using const DZ' 1661 dz1(1) = ztop /
real(npz) 1662 dz1(npz) = 0.5*dz1(1)
1669 call set_hybrid_z(is, ie, js, je, ng, npz, ztop, dz1, rgrav, &
1681 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
1682 delp, ak, bk, pt, delz, area, ng, mountain, hydrostatic, hybrid_z, domain)
1687 if (is_master()) print*,
'TEST TRACER enabled for this test case' 1690 ncnst, npz, q, agrid(is:ie,js:je,1), agrid(is:ie,js:je,2), 9., 9.)
1697 p1(1) = 205.*pi/180.
1701 p2(1) = agrid(i,j,1)
1702 p2(2) = agrid(i,j,2)
1704 if (r < r0 .and. .not.( abs(p1(2)-p2(2)) < 1./18. .and. p2(1)-p1(1) < 5./36.) .and. k > 16)
then 1720 cl = get_tracer_index(model_atmos,
'cl')
1721 cl2 = get_tracer_index(model_atmos,
'cl2')
1722 if (cl > 0 .and. cl2 > 0)
then 1724 q, delp,ncnst,agrid(isd:ied,jsd:jed,1),agrid(isd:ied,jsd:jed,2),bd)
1725 call mpp_update_domains(q,domain)
1737 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
1746 peln(i,1,j) = log(ptop)
1747 pk(i,j,1) = ptop**kappa
1752 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
1753 pk(i,j,k) = exp( kappa*log(pe(i,k,j)) )
1754 peln(i,k,j) = log(pe(i,k,j))
1763 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
1771 eta(k) = 0.5*( (ak(k)+ak(k+1))/1.e5 + bk(k)+bk(k+1) )
1772 eta_v(k) = (eta(k) - eta_0)*pi*0.5
1775 if ( .not. adiabatic )
then 1777 sphum = get_tracer_index(model_atmos,
'sphum')
1788 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j)) - 100000.
1789 q(i,j,k,
sphum) = 0.021*exp(-(agrid(i,j,2)/pcen(2))**4.)*exp(-(ptmp/34000.)**2.)
1818 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j+1,2))**2.0
1821 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1822 vv1 = utmp*(ee2(2,i,j+1)*cos(grid(i,j+1,1)) - ee2(1,i,j+1)*sin(grid(i,j+1,1)))
1824 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1827 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1828 vv3 = utmp*(ee2(2,i,j)*cos(grid(i,j,1)) - ee2(1,i,j)*sin(grid(i,j,1)))
1830 p1(:) = grid(i ,j ,1:2)
1831 p2(:) = grid(i,j+1 ,1:2)
1833 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1836 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1837 vv2 = utmp*(ew(2,i,j,2)*cos(pa(1)) - ew(1,i,j,2)*sin(pa(1)))
1839 v(i,j,z) = 0.25*(vv1 + 2.*vv2 + vv3)
1844 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i,j,2))**2.0
1847 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1848 uu1 = utmp*(ee1(2,i,j)*cos(grid(i,j,1)) - ee1(1,i,j)*sin(grid(i,j,1)))
1850 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*grid(i+1,j,2))**2.0
1853 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1854 uu3 = utmp*(ee1(2,i+1,j)*cos(grid(i+1,j,1)) - ee1(1,i+1,j)*sin(grid(i+1,j,1)))
1856 p1(:) = grid(i ,j ,1:2)
1857 p2(:) = grid(i+1,j ,1:2)
1859 utmp =
ubar * cos(eta_v(z))**(3.0/2.0) * sin(2.0*pa(2))**2.0
1862 if (-(r/r0)**2.0 > -40.0) utmp = utmp + u1*exp(-(r/r0)**2.0)
1863 uu2 = utmp*(es(2,i,j,1)*cos(pa(1)) - es(1,i,j,1)*sin(pa(1)))
1865 u(i,j,z) = 0.25*(uu1 + 2.*uu2 + uu3)
1880 eta(z) = 0.5*( (ak(z)+ak(z+1))/1.e5 + bk(z)+bk(z+1) )
1882 t_mean = t_0 * eta(z)**(rdgas*lapse_rate/grav)
1883 if (eta_t > eta(z)) t_mean = t_mean + delta_t*(eta_t - eta(z))**5.0
1885 230
format(i4.4,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14,
' ',e21.14)
1888 press = press + delp(is,js,zz)
1890 if (is_master())
write(*,230) z, eta(z), press/100., t_mean
1894 pt1 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1895 ( -2.0*(sin(agrid(i,j,2))**6.0) *(cos(agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1896 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1897 ( (8.0/5.0)*(cos(agrid(i,j,2))**3.0)*(sin(agrid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1908 pt2 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1909 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1910 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1911 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1913 pt3 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1914 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1915 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1916 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1918 pt4 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1919 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1920 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1921 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1923 pt5 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1924 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1925 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1926 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1928 pt6 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1929 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1930 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1931 ( (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 )
1932 pt7 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1933 ( -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 ) * &
1934 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1935 ( (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 )
1936 pt8 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1937 ( -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 ) * &
1938 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1939 ( (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 )
1940 pt9 = t_mean + 0.75*(eta(z)*pi*
ubar/rdgas)*sin(eta_v(z))*sqrt(cos(eta_v(z))) * ( &
1941 ( -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 ) * &
1942 2.0*
ubar*cos(eta_v(z))**(3.0/2.0) + &
1943 ( (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 )
1944 pt(i,j,z) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
1951 if ( (r/r0)**2 < 40. )
then 1952 pt(i,j,z) = pt(i,j,z) + pt0*exp(-(r/r0)**2)
1959 if (is_master()) print*,
' ' 1966 pt1 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1967 ( -2.0*(sin(agrid(i,j,2))**6.0) *(cos(agrid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1968 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1969 ( (8.0/5.0)*(cos(agrid(i,j,2))**3.0)*(sin(agrid(i,j,2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1980 pt2 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1981 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1982 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1983 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1985 pt3 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1986 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1987 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1988 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1990 pt4 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1991 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1992 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1993 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
1995 pt5 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
1996 ( -2.0*(sin(p1(2))**6.0) *(cos(p1(2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
1997 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
1998 ( (8.0/5.0)*(cos(p1(2))**3.0)*(sin(p1(2))**2.0 + 2.0/3.0) - pi/4.0 )*
radius*omega )
2000 pt6 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
2001 ( -2.0*(sin(grid(i,j,2))**6.0) *(cos(grid(i,j,2))**2.0 + 1.0/3.0) + 10.0/63.0 ) * &
2002 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
2003 ( (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 )
2004 pt7 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
2005 ( -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 ) * &
2006 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
2007 ( (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 )
2008 pt8 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
2009 ( -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 ) * &
2010 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
2011 ( (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 )
2012 pt9 =
ubar* (cos( (eta_s-eta_0)*pi/2.0 ))**(3.0/2.0) * ( &
2013 ( -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 ) * &
2014 ubar*cos( (eta_s-eta_0)*pi/2.0 )**(3.0/2.0) + &
2015 ( (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 )
2016 phis(i,j) = 0.25*pt1 + 0.125*(pt2+pt3+pt4+pt5) + 0.0625*(pt6+pt7+pt8+pt9)
2023 if ( .not.hydrostatic )
then 2029 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
2035 if (.not. adiabatic)
then 2036 zvir = rvgas/rdgas - 1.
2041 pt(i,j,k) = pt(i,j,k)/(1. + zvir*q(i,j,k,
sphum))
2051 write(stdout(), *)
'PI:', pi
2052 write(stdout(), *)
'PHIS:', mpp_chksum(phis(is:ie,js:je))
2057 is,ie,js,je,isd,ied,jsd,jed,npz,ncnst,ak,bk,ptop, &
2058 pk,peln,pe,pkz,gz,phis,ps,grid,agrid,hydrostatic, &
2059 nwat, adiabatic,
test_case == -13, domain, bd)
2061 write(stdout(), *)
'PHIS:', mpp_chksum(phis(is:ie,js:je))
2087 ze1(k) = ze1(k+1) + ztop/
real(npz)
2091 ze1(1) = ztop + 1.5*ztop/
real(npz)
2104 delz(i,j,k) = ze1(k+1) - ze1(k)
2105 pk(i,j,k) = pk(i,j,k+1) + grav*delz(i,j,k)/(cp_air*t00)*pk0
2106 pe(i,k,j) = pk(i,j,k)**(1./kappa)
2112 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
2117 peln(i,k,j) = log(pe(i,k,j))
2126 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2127 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2140 r0 = 0.5*(ze1(k)+ze1(k+1)) - 3.2e3
2142 r0 = (0.5*(ze1(k)+ze1(k+1)) - 3.0e3) / 2.e3
2147 p2(1) = agrid(i,j,1)
2148 p2(2) = agrid(i,j,2)
2151 dist = sqrt( r**2 + r0**2 ) / 3.2e3
2154 dist = sqrt( r**2 + r0**2 )
2156 if ( dist<=1. )
then 2157 q(i,j,k,1) = pk0 * pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
2158 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
2163 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
2180 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2183 pe1(z) = ak(z) + bk(z)*p00
2188 ze1(z) = ze1(z+1) + ztop/
real(npz)
2192 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2196 ps(i,j) = pe1(npz+1)
2204 peln(i,z,j) = log(pe1(z))
2205 pk(i,j,z) = exp(kappa*peln(i,z,j))
2218 vort(i,j) = 0.5*(1.+cos(pi*r/r0))
2229 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*pi/ztop )
2232 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(kappa*(peln(i,z+1,j)-peln(i,z,j)))
2233 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2235 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2236 q(i,j,z,1) = q(i,j,z,1) + vort(i,j)*zmid
2249 call gw_1d(npz, p00, ak, bk, ptop, ztop, ppt)
2252 pe1(z) = ak(z) + bk(z)*p00
2257 ze1(z) = ze1(z+1) + ztop/
real(npz)
2261 if ( is_master() )
write(*,*)
'Model top (pa)=', ptop
2265 ps(i,j) = pe1(npz+1)
2273 peln(i,z,j) = log(pe1(z))
2274 pk(i,j,z) = exp(kappa*peln(i,z,j))
2287 vort(i,j) = 0.5*(1.+cos(pi*r/r0))
2297 zmid = sin( 0.5*(ze1(z)+ze1(z+1))*pi/ztop )
2300 pkz(i,j,z) = (pk(i,j,z+1)-pk(i,j,z))/(kappa*(peln(i,z+1,j)-peln(i,z,j)))
2301 delp(i,j,z) = pe(i,z+1,j)-pe(i,z,j)
2303 pt(i,j,z) = ( ppt(z) + pturb*vort(i,j)*zmid ) * pkz(i,j,z)
2311 n2 = grav**2 / (cp_air*pt0)
2320 phis(i,j) = grav*2.e3*exp( -(r/1500.e3)**2 )
2322 (sin(agrid(i,j,2))**2-1.) - n2/(grav*grav*kappa)*phis(i,j))
2330 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
2336 p1(:) = grid(i ,j ,1:2)
2337 p2(:) = grid(i,j+1 ,1:2)
2341 utmp =
ubar * cos(p3(2))
2350 p1(:) = grid(i, j,1:2)
2351 p2(:) = grid(i+1,j,1:2)
2355 utmp =
ubar * cos(p3(2))
2381 p1(:) = grid(i ,j ,1:2)
2382 p2(:) = grid(i,j+1 ,1:2)
2386 utmp =
ubar * cos(p3(2))
2393 p1(:) = grid(i, j,1:2)
2394 p2(:) = grid(i+1,j,1:2)
2398 utmp =
ubar * cos(p3(2))
2419 p1(1) = (0.5-0.125) * pi
2426 p2(:) = agrid(i,j,1:2)
2428 if ( r < pi*
radius )
then 2429 p4(:) = p2(:) - p1(:)
2430 if ( abs(p4(1)) > 1.e-12 )
then 2431 zeta = asin( p4(2) / sqrt(p4(1)**2 + p4(2)**2) )
2435 if ( p4(1) <= 0. ) zeta = pi - zeta
2437 v1 = r/uu1 * cos( zeta )
2438 v2 = r/uu2 * sin( zeta )
2439 phis(i,j) = ftop / ( 1. + v1**2 + v2**2 )
2446 if ( hybrid_z )
then 2450 elseif( npz.eq.31 .or. npz.eq.41 .or. npz.eq.51 )
then 2454 if ( is_master() )
write(*,*)
'Using const DZ' 2456 dz1(1) = ztop /
real(npz) 2461 dz1(1) = max( 1.0e3, 3.*dz1(2) )
2467 ze1(k) = ze1(k+1) + dz1(k)
2471 call set_hybrid_z( is, ie, js, je, ng, npz, ztop, dz1, rgrav, &
2474 call mpp_error(fatal,
'This test case is only currently setup for hybrid_z')
2480 delz(i,j,k) = ze0(i,j,k+1) - ze0(i,j,k)
2490 s0 = grav*grav / (cp_air*n2)
2494 pe1(k) = p00*( (1.-s0/t00) + s0/t00*exp(-n2*ze1(k)/grav) )**(1./kappa)
2498 if ( is_master() )
write(*,*)
'Lee vortex testcase: model top (mb)=', ptop/100.
2504 bk(k) = (pe1(k) - pe1(1)) / (pe1(npz+1)-pe1(1))
2505 ak(k) = pe1(1)*(1.-bk(k))
2514 pk(i,j,k) = pk0 - (1.-exp(-n2/grav*ze0(i,j,k))) * (grav*grav)/(n2*cp_air*pt0)
2515 pe(i,k,j) = pk(i,j,k) ** (1./kappa)
2516 peln(i,k,j) = log(pe(i,k,j))
2524 peln(i,1,j) = log(pe(i,1,j))
2525 pk(i,j,1) = pe(i,1,j) ** kappa
2526 ps(i,j) = pe(i,npz+1,j)
2533 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2534 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
2535 pt(i,j,k) = pkz(i,j,k)*grav*delz(i,j,k) / ( cp_air*(pk(i,j,k)-pk(i,j,k+1)) )
2546 if (.not.hydrostatic) w(:,:,:)= 0.0
2559 dz = 12000./
real(npz)
2561 allocate(
zz0(npz+1))
2562 allocate(
pz0(npz+1))
2570 if (is_master()) print*,
'TRACER ADVECTION TEST CASE' 2571 if (is_master()) print*,
'INITIAL LEVELS' 2577 if (is_master())
write(*,*) k,
pz0(k),
zz0(k)
2590 ptmp = 0.5*(
pz0(k) +
pz0(k+1))
2597 delp(i,j,k) =
pz0(k+1)-
pz0(k)
2602 ptop = 100000.*exp(-12000.*grav/t00/rdgas)
2610 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
2613 call mpp_update_domains( psi, domain )
2616 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
2617 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
2625 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
2626 if (dist==0) vc(i,j,k) = 0.
2632 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
2633 if (dist==0) uc(i,j,k) = 0.
2640 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
2641 if (dist==0) v(i,j,k) = 0.
2647 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
2648 if (dist==0) u(i,j,k) = 0.
2654 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
2655 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
2657 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
2658 if (dist==0) ua(i,j,k) = 0.
2659 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
2660 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
2662 va(i,j,k) = (psi2 - psi1) / (dist)
2663 if (dist==0) va(i,j,k) = 0.
2670 uc(:,:,k) = uc(:,:,1)
2671 vc(:,:,k) = vc(:,:,1)
2672 ua(:,:,k) = ua(:,:,1)
2673 va(:,:,k) = va(:,:,1)
2676 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
2677 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
2685 call mpp_error(fatal,
'Value of tracer_test not implemented ')
2703 if (.not.hydrostatic) w(:,:,:)= 0.0
2707 dz = 12000./
real(npz)
2712 exponent = rdgas*gamma/grav
2713 px = ((t00-9000.*gamma)/t00)**(1./exponent)
2717 height = 12000. - dz*
real(k-1)
2718 if (height >= 9000. )
then 2719 ak(k) = p00*((t00-height*gamma)/t00)**(1./exponent)
2722 ak(k) = (((t00-height*gamma)/t00)**(1./exponent)-1.)/(px - 1.)*px*p00
2723 bk(k) = (((t00-height*gamma)/t00)**(1./exponent)-px)/(1.-px)
2725 if (is_master())
write(*,*) k, ak(k), bk(k), height, ak(k)+bk(k)*p00
2731 p1(1) = 3.*pi/2. ; p1(2) = 0.
2738 p2(:) = agrid(i,j,1:2)
2741 phis(i,j) = grav*0.5*2000.*(1. + cos(pi*r/r0))*cos(pi*r/zetam)**2.
2742 pe(i,npz+1,j) = p00*(1.-gamma/t00*phis(i,j)/grav)**(1./exponent)
2747 ps(i,j) = pe(i,npz+1,j)
2754 pe(i,k,j) = ak(k) + bk(k)*ps(i,j)
2755 gz(i,j,k) = t00/gamma*(1. - (pe(i,k,j)/p00)**exponent)
2767 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
2772 pt(i,j,k) = -grav*t00*p00/(rdgas*gamma + grav)/delp(i,j,k) * &
2773 ( (pe(i,k,j)/p00)**(exponent+1.) - (pe(i,k+1,j)/p00)**(exponent+1.) )
2789 zvir = rvgas/rdgas - 1.
2795 pk(i,j,1) = ptop**kappa
2797 peln(i,1,j) = log(ptop)
2804 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
2805 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
2806 peln(i,k+1,j) = log(pe(i,k+1,j))
2807 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
2815 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
2825 pp0(1) = 262.0/180.*pi
2826 pp0(2) = 35.0/180.*pi
2833 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
2840 ze1(k) = ze1(k+1) - delz(is,js,k)
2844 if (is_master())
then 2846 write(6,*)
'Toy supercell winds, piecewise approximation' 2848 write(6,*)
'Toy supercell winds, tanh approximation' 2853 zm = 0.5*(ze1(k)+ze1(k+1))
2858 if ( zm .le. 2.e3 )
then 2859 utmp = 8.*(1.-cos(pi*zm/4.e3))
2860 vtmp = 8.*sin(pi*zm/4.e3)
2861 elseif (zm .le. 6.e3 )
then 2862 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
2872 utmp = 15.0*(1.+tanh(zm/2000. - 1.5))
2873 vtmp = 8.5*tanh(zm/1000.)
2888 if( is_master() )
then 2889 write(6,*) k, utmp, vtmp
2894 p1(:) = grid(i ,j ,1:2)
2895 p2(:) = grid(i,j+1 ,1:2)
2906 p1(:) = grid(i, j,1:2)
2907 p2(:) = grid(i+1,j,1:2)
2918 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
2919 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
2920 .true., hydrostatic, nwat, domain, adiabatic)
2927 zm = 0.5*(ze1(k)+ze1(k+1))
2928 ptmp = ( (zm-zc)/zc ) **2
2929 if ( ptmp < 1. )
then 2933 if ( dist < 1. )
then 2934 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
2943 call mpp_error(fatal,
' test_case 32 not yet implemented')
2969 p0(1) = 60./180. * pi
2984 p2(2) = agrid(i,j,2)
2985 pt1 = exp( dum*(us0*sin(p2(2)))**2 )
2987 pt2 = exp( dum*(us0*sin(p2(2)))**2 )
2989 pt3 = exp( dum*(us0*sin(p2(2)))**2 )
2991 pt4 = exp( dum*(us0*sin(p2(2)))**2 )
2993 pt5 = exp( dum*(us0*sin(p2(2)))**2 )
2995 pt6 = exp( dum*(us0*sin(p2(2)))**2 )
2996 p2(2) = grid(i+1,j,2)
2997 pt7 = exp( dum*(us0*sin(p2(2)))**2 )
2998 p2(2) = grid(i+1,j+1,2)
2999 pt8 = exp( dum*(us0*sin(p2(2)))**2 )
3000 p2(2) = grid(i,j+1,2)
3001 pt9 = exp( dum*(us0*sin(p2(2)))**2 )
3002 ptmp = t00*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
3004 ptmp = t00*exp( dum*(us0*sin(agrid(i,j,2)))**2 )
3020 p2(1:2) = agrid(i,j,1:2)
3022 pt1 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3025 pt2 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3028 pt3 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3031 pt4 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3034 pt5 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3035 p2(1:2) = grid(i,j,1:2)
3037 pt6 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3038 p2(1:2) = grid(i+1,j,1:2)
3040 pt7 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3041 p2(1:2) = grid(i+1,j+1,1:2)
3043 pt8 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3044 p2(1:2) = grid(i,j+1,1:2)
3046 pt9 = cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3047 phis(i,j) = grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
3049 p2(1:2) = agrid(i,j,1:2)
3051 phis(i,j) = grav*h0*cos(p2(2))*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3068 pt1 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3071 pt2 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3074 pt3 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3077 pt4 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3080 pt5 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3082 pt6 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3084 pt7 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3086 pt8 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3088 pt9 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3089 phis(i,j) = grav*h0*(0.25*pt1+0.125*(pt2+pt3+pt4+pt5)+0.0625*(pt6+pt7+pt8+pt9))
3092 pt1 = exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3093 phis(i,j) = grav*h0*exp(-(r/5.e3)**2)*cos(pi*r/4.e3)**2
3102 ps(i,j) = p00*exp( -0.5*(us0*sin(agrid(i,j,2)))**2/(rdgas*t00)-phis(i,j)/(rdgas*pt(i,j,1)) )
3104 peln(i,1,j) = log(ptop)
3105 pk(i,j,1) = ptop**kappa
3112 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
3113 peln(i,k,j) = log(pe(i,k,j))
3114 pk(i,j,k) = exp( kappa*peln(i,k,j) )
3122 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3123 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(peln(i,k,j)-peln(i,k+1,j))
3131 ze1(npz+1) = phis(i,j)/grav
3133 ze1(k) = ze1(k+1) - delz(i,j,k)
3136 w(i,j,k) = 0.5*(ze1(k)+ze1(k+1))
3140 call mpp_update_domains( w, domain )
3145 p1(:) = grid(i ,j, 1:2)
3146 p2(:) = grid(i,j+1, 1:2)
3151 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i-1,j,k)+w(i,j,k)) )
3157 p1(:) = grid(i, j, 1:2)
3158 p2(:) = grid(i+1,j, 1:2)
3162 utmp = us0*cos(p3(2))*sqrt( 1. + cs_m3*(w(i,j-1,k)+w(i,j,k)) )
3171 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3172 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
3173 .true., hydrostatic, nwat, domain, adiabatic)
3185 zvir = rvgas/rdgas - 1.
3198 ze1(k) = ze1(k+1) + ztop/
real(npz)
3207 call var_dz(npz, ztop, ze1)
3210 zs1(k) = 0.5*(ze1(k)+ze1(k+1))
3218 ts1(k) = cp_air*ts1(k)*(1.+zvir*qs1(k))
3224 call balanced_k(npz, is, ie, js, je, ng, pe1(npz+1), ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
3225 delz, zvir, ptop, ak, bk, agrid)
3228 ps(i,j) = pe(i,npz+1,j)
3235 peln(i,k,j) = log(pe(i,k,j))
3236 pk(i,j,k) = exp( kappa*peln(i,k,j) )
3244 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
3245 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3254 delz(i,j,k) = rdgas/grav*pt(i,j,k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
3261 delz(i,j,k) = ze1(k+1) - ze1(k)
3262 pt(i,j,k) = delz(i,j,k)*grav/(rdgas*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j)))
3271 p1(:) = grid(i ,j ,1:2)
3272 p2(:) = grid(i,j+1 ,1:2)
3276 v(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e2,ex)
3281 p1(:) = grid(i, j,1:2)
3282 p2(:) = grid(i+1,j,1:2)
3286 u(i,j,k) = uz1(k)*cos(p3(2))*
inner_prod(e1,ex)
3303 zm = 0.5*(ze1(k)+ze1(k+1))
3304 ptmp = ( (zm-zc)/zc ) **2
3305 if ( ptmp < 1. )
then 3310 if ( dist < 1. )
then 3311 pt(i,j,k) = pt(i,j,k) + (pkz(i,j,k)/pk0)*pturb*cos(0.5*pi*dist)**2
3332 zvir = rvgas/rdgas - 1.
3339 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3347 peln(i,1,j) = log(pe(i,1,j))
3348 pk(i,j,1) = exp(kappa*peln(i,1,j))
3352 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3353 peln(i,k,j) = log(pe(i,k,j))
3354 pk(i,j,k) = exp(kappa*peln(i,k,j))
3366 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3368 if ( dist .le. r0 )
then 3380 if (.not.hydrostatic)
then 3384 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))
3406 zvir = rvgas/rdgas - 1.
3413 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
3421 peln(i,1,j) = log(pe(i,1,j))
3422 pk(i,j,1) = exp(kappa*peln(i,1,j))
3426 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3427 peln(i,k,j) = log(pe(i,k,j))
3428 pk(i,j,k) = exp(kappa*peln(i,k,j))
3449 p1(:) = grid(i ,j ,1:2)
3450 p2(:) = grid(i,j+1 ,1:2)
3453 utmp =
ubar*exp(-(r/r0)**2)
3461 p1(:) = grid(i, j,1:2)
3462 p2(:) = grid(i+1,j,1:2)
3465 utmp =
ubar*exp(-(r/r0)**2)
3474 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
3476 pt(i,j,k) = pt0/p00**kappa
3478 pt(i,j,k) = pt(i,j,k)*pkz(i,j,k)
3492 q(i,j,k,:) = agrid(i,j,1)*0.180/pi
3498 ncnst, npz, q, agrid(is:ie,js:je,1), agrid(is:ie,js:je,2), 9., 9.)
3501 if ( .not. hydrostatic )
then 3505 delz(i,j,k) = rdgas*pt(i,j,k)/grav*log(pe(i,k,j)/pe(i,k+1,j))
3522 p0(1) = 180. * pi / 180.
3523 p0(2) = 10. * pi / 180.
3538 p2(:) = agrid(i,j,1:2)
3540 ps(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3545 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3551 delp(i,j,z) = ak(z+1)-ak(z) + ps(i,j)*(bk(z+1)-bk(z))
3563 pe(i,k,j) = pe(i,k-1,j) + delp(i,j,k-1)
3571 p2(:) = 0.5*(grid(i,j,1:2)+grid(i,j+1,1:2))
3573 ps_v(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3578 p2(:) = 0.5*(grid(i,j,1:2)+grid(i+1,j,1:2))
3580 ps_u(i,j) = p00 - dp*exp(-(r/rp)**1.5)
3591 pe_v(i,k,j) = ak(k) + ps_v(i,j)*bk(k)
3601 pe_u(i,k,j) = ak(k) + ps_u(i,j)*bk(k)
3610 zvir = rvgas/rdgas - 1.
3613 p0 = (/ pi, pi/18. /)
3620 t00 = ts0*(1.+zvir*q00)
3621 exponent = rdgas*gamma/grav
3625 cor = 2.*omega*sin(p0(2))
3630 p1(:) = grid(i ,j ,1:2)
3631 p2(:) = grid(i,j+1 ,1:2)
3636 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3637 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3638 d = max(1.e-15,sqrt(d1**2+d2**2))
3643 ptmp = 0.5*(pe_v(i,k,j)+pe_v(i,k+1,j))
3644 height = (t00/gamma)*(1.-(ptmp/ps_v(i,j))**exponent)
3645 if (height > ztrop)
then 3648 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3649 - exppr*(r/rp)**exppr*rdgas*(t00-gamma*height) &
3650 /(exppz*height*rdgas*(t00-gamma*height)/(grav*zp**exppz) &
3651 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3663 p1(:) = grid(i, j,1:2)
3664 p2(:) = grid(i+1,j,1:2)
3669 d1 = sin(p0(2))*cos(p3(2)) - cos(p0(2))*sin(p3(2))*cos(p3(1)-p0(1))
3670 d2 = cos(p0(2))*sin(p3(1)-p0(1))
3671 d = max(1.e-15,sqrt(d1**2+d2**2))
3676 ptmp = 0.5*(pe_u(i,k,j)+pe_u(i,k+1,j))
3677 height = (t00/gamma)*(1.-(ptmp/ps_u(i,j))**exponent)
3678 if (height > ztrop)
then 3681 utmp = 1.d0/d*(-cor*r/2.d0+sqrt((cor*r/2.d0)**(2.d0) &
3682 - exppr*(r/rp)**exppr*rdgas*(t00-gamma*height) &
3683 /(exppz*height*rdgas*(t00-gamma*height)/(grav*zp**exppz) &
3684 +(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz)))))
3696 ttrop = t00 - gamma*ztrop
3705 ptmp = 0.5*(pe(i,k,j)+pe(i,k+1,j))
3706 height = (t00/gamma)*(1.-(ptmp/ps(i,j))**exponent)
3707 if (height > ztrop)
then 3711 q(i,j,k,1) = q00*exp(-height/zq1)*exp(-(height/zq2)**exppz)
3712 p2(:) = agrid(i,j,1:2)
3714 pt(i,j,k) = (t00-gamma*height)/(1.d0+zvir*q(i,j,k,1))/(1.d0+exppz*rdgas*(t00-gamma*height)*height &
3715 /(grav*zp**exppz*(1.d0-p00/dp*exp((r/rp)**exppr)*exp((height/zp)**exppz))))
3724 ps(i,j) = pe(i,npz+1,j)
3728 if (.not.hydrostatic)
then 3732 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))
3739 call dtoa(u , v , ua, va, dx,dy,dxa,dya,dxc,dyc,npx, npy, ng, bd)
3741 call prt_maxmin(
'PS', ps(is:ie,js:je), is, ie, js, je, 0, 1, 0.01)
3759 call dcmip16_tc (delp, pt, u, v, q, w, delz, &
3760 is, ie, js, je, isd, ied, jsd, jed, npz, ncnst, &
3761 ak, bk, ptop, pk, peln, pe, pkz, gz, phis, &
3762 ps, grid, agrid, hydrostatic, nwat, adiabatic)
3766 call mpp_error(fatal,
" test_case not defined" )
3770 call mpp_update_domains( phis, domain )
3772 ftop =
g_sum(domain, phis(is:ie,js:je), is, ie, js, je, ng, area, 1)
3773 if(is_master())
write(*,*)
'mean terrain height (m)=', ftop/grav
3777 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
3778 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., mountain, &
3779 moist_phys, hydrostatic, nwat, domain, adiabatic, .not.hydrostatic)
3782 #ifdef COLUMN_TRACER 3783 if( ncnst>1 ) q(:,:,:,2:ncnst) = 0.0
3791 p1(:) = grid(i ,j ,1:2)
3792 p2(:) = grid(i,j+1 ,1:2)
3798 if (-(r/r0)**2.0 > -40.0) q(i,j,z,1) = exp(-(r/r0)**2.0)
3834 nullify(cubed_sphere)
3839 nullify(have_south_pole)
3840 nullify(have_north_pole)
3849 subroutine get_vorticity(isc, iec, jsc, jec ,isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
3850 integer isd, ied, jsd, jed, npz
3851 integer isc, iec, jsc, jec
3852 real,
intent(in) :: u(isd:ied, jsd:jed+1, npz), v(isd:ied+1, jsd:jed, npz)
3853 real,
intent(out) :: vort(isc:iec, jsc:jec, npz)
3854 real,
intent(IN) :: dx(isd:ied,jsd:jed+1)
3855 real,
intent(IN) :: dy(isd:ied+1,jsd:jed)
3856 real,
intent(IN) :: rarea(isd:ied,jsd:jed)
3858 real :: utmp(isc:iec, jsc:jec+1), vtmp(isc:iec+1, jsc:jec)
3864 utmp(i,j) = u(i,j,k)*dx(i,j)
3869 vtmp(i,j) = v(i,j,k)*dy(i,j)
3875 vort(i,j,k) = rarea(i,j)*(utmp(i,j)-utmp(i,j+1)-vtmp(i,j)+vtmp(i+1,j))
3882 subroutine checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, &
3883 nq, km, q, lon, lat, nx, ny, rn)
3893 integer,
intent(in):: nq
3894 integer,
intent(in):: km
3895 integer,
intent(in):: i0, i1
3896 integer,
intent(in):: j0, j1
3897 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3898 real,
intent(in):: nx
3899 real,
intent(in):: ny
3900 real,
intent(in),
optional:: rn
3901 real(kind=R_GRID),
intent(in),
dimension(i0:i1,j0:j1):: lon, lat
3902 real,
intent(out):: q(ifirst:ilast,jfirst:jlast,km,nq)
3904 real:: qt(i0:i1,j0:j1)
3912 qtmp = sin(nx*lon(i,j))*sin(ny*lat(i,j))
3913 if ( qtmp < 0. )
then 3921 if (
present(rn) )
then 3929 call random_number(ftmp)
3930 q(i,j,k,iq) = qt(i,j) + rn*ftmp
3942 q(i,j,k,iq) = qt(i,j)
3952 km, q, delp, ncnst, lon, lat, bd)
3959 integer,
intent(in):: km
3960 integer,
intent(in):: i0, i1
3961 integer,
intent(in):: j0, j1
3962 integer,
intent(in):: ifirst, ilast, jfirst, jlast
3963 integer,
intent(in):: ncnst
3964 real(kind=R_GRID),
intent(in),
dimension(ifirst:ilast,jfirst:jlast):: lon, lat
3965 real,
intent(inout):: q(ifirst:ilast,jfirst:jlast,km,ncnst)
3966 real,
intent(in):: delp(ifirst:ilast,jfirst:jlast,km)
3968 real:: D, k1, r, ll, sinthc, costhc, mm
3974 real,
parameter :: qcly = 4.e-6
3975 real,
parameter :: lc = 5.*pi/3.
3976 real,
parameter :: thc = pi/9.
3977 real,
parameter :: k2 = 1.
3982 cl = get_tracer_index(model_atmos,
'Cl')
3983 cl2 = get_tracer_index(model_atmos,
'Cl2')
3987 k1 = max(0., sin(lat(i,j))*sinthc + cos(lat(i,j))*costhc*cos(lon(i,j) - lc))
3989 d = sqrt(r*r + 2.*r*qcly)
3991 q(i,j,1,cl2) = 0.5*(qcly - q(i,j,1,cl))
3998 q(i,j,k,cl) = q(i,j,1,cl)
3999 q(i,j,k,cl2) = q(i,j,1,cl2)
4006 if (is_master())
then 4011 qcly0 =
qcly0 + (q(i,j,k,cl) + 2.*q(i,j,k,cl2))*delp(i,j,k)
4012 mm = mm + delp(i,j,k)
4017 if (is_master()) print*,
' qcly0 = ',
qcly0 4028 real,
intent(in):: ubar
4029 real,
intent(in):: r0
4030 real,
intent(in):: p1(2)
4031 real,
intent(inout):: u(bd%isd:bd%ied, bd%jsd:bd%jed+1)
4032 real,
intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed)
4033 real(kind=R_GRID),
intent(IN) :: grid(bd%isd:bd%ied+1,bd%jsd:bd%jed+1,2)
4035 real(kind=R_GRID):: p2(2), p3(2), p4(2)
4036 real(kind=R_GRID):: e1(3), e2(3), ex(3), ey(3)
4037 real:: vr, r, d2, cos_p, x1, y1
4041 integer :: is, ie, js, je
4042 integer :: isd, ied, jsd, jed
4058 p2(1) = p2(1) - p1(1)
4059 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
4067 x1 = cos(p2(2))*sin(p2(1))
4068 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
4069 d2 = max(1.e-25, sqrt(x1**2 + y1**2))
4072 p3(1) = grid(i,j, 1) - p1(1)
4073 p3(2) = grid(i,j, 2)
4074 p4(1) = grid(i+1,j,1) - p1(1)
4075 p4(2) = grid(i+1,j,2)
4087 p2(1) = p2(1) - p1(1)
4088 cos_p = sin(p2(2))*sin(p1(2)) + cos(p2(2))*cos(p1(2))*cos(p2(1))
4095 x1 = cos(p2(2))*sin(p2(1))
4096 y1 = sin(p2(2))*cos(p1(2)) - cos(p2(2))*sin(p1(2))*cos(p2(1))
4097 d2 = max(1.e-25, sqrt(x1**2 + y1**2))
4100 p3(1) = grid(i,j, 1) - p1(1)
4101 p3(2) = grid(i,j, 2)
4102 p4(1) = grid(i,j+1,1) - p1(1)
4103 p4(2) = grid(i,j+1,2)
4113 real function gh_jet(npy, lat_in)
4114 integer,
intent(in):: npy
4115 real,
intent(in):: lat_in
4116 real lat, lon, dp, uu
4123 dp = pi /
real(jm-1)
4133 lat = -pi/2. + (
real(j-1)-0.5)*dp
4135 ft = 2.*omega*sin(lat)
4160 real function u_jet(lat)
4162 real umax, en, ph0, ph1
4167 en = exp( -4./(ph1-ph0)**2 )
4169 if ( lat>ph0 .and. lat<ph1 )
then 4170 u_jet = (umax/en)*exp( 1./( (lat-ph0)*(lat-ph1) ) )
4176 subroutine get_case9_b(B, agrid, isd, ied, jsd, jed)
4177 integer,
intent(IN) :: isd, ied, jsd, jed
4178 real,
intent(OUT) :: B(isd:ied,jsd:jed)
4179 real,
intent(IN) :: agrid(isd:ied,jsd:jed,2)
4187 if (sin(agrid(i,j,2)) > 0.)
then 4188 myc = sin(agrid(i,j,1))
4189 yy = (cos(agrid(i,j,2))/sin(agrid(i,j,2)))**2
4190 myb =
gh0*yy*exp(1.-yy)
4208 integer,
intent(IN) :: isd,ied,jsd,jed
4209 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4210 real ,
intent(IN) :: time_since_start
4216 tday = time_since_start/86400.0
4217 if (tday >= 20.)
then 4218 aoft(2) = 0.5*(1.-cos(0.25*pi*(tday-20)))
4219 if (tday == 24)
aoft(2) = 1.0
4220 elseif (tday <= 4.)
then 4221 aoft(2) = 0.5*(1.-cos(0.25*pi*tday))
4222 elseif (tday <= 16.)
then 4225 aoft(2) = 0.5*(1.+cos(0.25*pi*(tday-16.)))
4230 phis(i,j) = amean*
case9_b(i,j)
4243 integer,
intent(IN) :: isd,ied,jsd,jed
4244 real ,
intent(INOUT) :: phis(isd:ied ,jsd:jed )
4261 subroutine case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain, bd)
4264 real,
intent(INOUT) :: delp(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
4265 real,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz)
4266 real,
intent(INOUT) :: vc(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz)
4267 real,
intent(INOUT) :: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,npz)
4268 real,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,npz)
4269 real,
intent(INOUT) :: ua(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
4270 real,
intent(INOUT) :: va(bd%isd:bd%ied,bd%jsd:bd%jed,npz)
4271 real,
intent(INOUT) :: pe(bd%is-1:bd%ie+1, npz+1,bd%js-1:bd%je+1)
4272 real,
intent(IN) :: time, dt
4273 real,
intent(INOUT) :: ptop
4274 integer,
intent(IN) :: npx, npy, npz
4276 type(domain2d),
intent(INOUT) :: domain
4283 real :: s, l, dt2, V0, phase
4284 real :: ull, vll, lonp
4285 real :: p0(2), elon(3), elat(3)
4287 real :: psi(bd%isd:bd%ied,bd%jsd:bd%jed)
4288 real :: psi_b(bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
4289 real :: dist, psi1, psi2
4294 real(kind=R_GRID) :: e1(3), e2(3), ex(3), ey(3), pt(2), p1(2), p2(2), p3(2), rperiod, timefac, t00
4296 integer :: wind_field = 1
4298 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
4299 real,
pointer,
dimension(:,:) :: dx, dxa, dy, dya, dxc, dyc
4301 integer :: is, ie, js, je
4302 integer :: isd, ied, jsd, jed, ng
4314 agrid => gridstruct%agrid_64
4315 grid => gridstruct%grid_64
4318 dxa => gridstruct%dxa
4319 dxc => gridstruct%dxc
4321 dya => gridstruct%dya
4322 dyc => gridstruct%dyc
4324 period =
real( 12*24*3600 ) 4329 phase = pi*time/period
4339 select case (wind_field)
4342 omega0 = 23000.*pi/period
4345 ptop = 100000.*exp(-12000.*grav/t00/rdgas)
4350 s = min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*pi)))
4351 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(agrid(i,j,1)-period*(time+dt2))*cos(agrid(i,j,2))* &
4352 cos(period*(time+dt2))*sin(s*0.5*pi)
4360 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4371 cos(agrid(i,j,1))*cos(agrid(i,j,2))*sin(
alpha) ) )
4374 call mpp_update_domains( psi, domain )
4377 psi_b(i,j) = (-1.0 *
ubar *
radius *( sin(grid(i,j,2)) *cos(
alpha) - &
4378 cos(grid(i,j,1))*cos(grid(i,j,2))*sin(
alpha) ) )
4387 vc(i,j,k) = (psi_b(i+1,j)-psi_b(i,j))/dist
4388 if (dist==0) vc(i,j,k) = 0.
4394 uc(i,j,k) = -1.0*(psi_b(i,j+1)-psi_b(i,j))/dist
4395 if (dist==0) uc(i,j,k) = 0.
4402 v(i,j,k) = (psi(i,j)-psi(i-1,j))/dist
4403 if (dist==0) v(i,j,k) = 0.
4409 u(i,j,k) = -1.0*(psi(i,j)-psi(i,j-1))/dist
4410 if (dist==0) u(i,j,k) = 0.
4416 psi1 = 0.5*(psi(i,j)+psi(i,j-1))
4417 psi2 = 0.5*(psi(i,j)+psi(i,j+1))
4419 ua(i,j,k) = -1.0 * (psi2 - psi1) / (dist)
4420 if (dist==0) ua(i,j,k) = 0.
4421 psi1 = 0.5*(psi(i,j)+psi(i-1,j))
4422 psi2 = 0.5*(psi(i,j)+psi(i+1,j))
4424 va(i,j,k) = (psi2 - psi1) / (dist)
4425 if (dist==0) va(i,j,k) = 0.
4431 omega0 = 23000.*pi/period
4436 s = min(1.,2.*sqrt(sin((pe(i,k,j)-ptop)/(pe(i,npz+1,j)-ptop)*pi)))
4437 pe(i,k,j) = pe(i,k,j) + dt*omega0*sin(agrid(i,j,1)-period*(time+dt2))*cos(agrid(i,j,2))* &
4438 cos(period*(time+dt2))*sin(s*0.5*pi)
4446 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4456 p1(:) = grid(i ,j ,1:2)
4457 p2(:) = grid(i,j+1 ,1:2)
4461 l = p3(1) - 2.*pi*time/period
4462 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(pi*time/period) + 2.*pi*
radius/period*cos(p3(2))
4463 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(pi*time/period)
4469 p1(:) = grid(i, j,1:2)
4470 p2(:) = grid(i+1,j,1:2)
4474 l = p3(1) - 2.*pi*time/period
4475 utmp =
ubar * sin(l)**2 * sin(2.*p3(2)) * cos(pi*time/period) + 2.*pi*
radius/period*cos(p3(2))
4476 vtmp =
ubar * sin(2.*l) * cos(p3(2)) * cos(pi*time/period)
4499 call dtoa( u(:,:,1), v(:,:,1),ua(:,:,1),va(:,:,1),dx,dy,dxa,dya,dxc,dyc,npx,npy,ng,bd)
4500 call mpp_update_domains( ua, va, domain, gridtype=agrid_param)
4501 call atoc(ua(:,:,1),va(:,:,1),uc(:,:,1),vc(:,:,1),dx,dy,dxa,dya,npx,npy,ng, gridstruct%bounded_domain, domain, bd)
4506 ua(i,j,k) = ua(i,j,1)
4511 va(i,j,k) = va(i,j,1)
4519 vc(i,j,k) = vc(i,j,1)
4524 uc(i,j,k) = uc(i,j,1)
4537 pe(i,k,j) = pe(i,k,j) + dt*omega0*grav*pe(i,k,j)/rdgas/300./k_cell* &
4538 (-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)))* &
4539 sin(pi*
zz0(k)/12000.)*cos(phase)
4547 delp(i,j,k) = pe(i,k+1,j) - pe(i,k,j)
4558 utmp =
ubar*cos(agrid(i,j,2))
4559 vtmp = -
radius * omega0 * pi / k_cell / 12000. * &
4560 cos(agrid(i,j,2)) * sin(k_cell * agrid(i,j,2)) * &
4561 sin(pi*
zz0(k)/12000.)*cos(phase)
4570 uc(:,:,k) = uc(:,:,1)
4571 vc(:,:,k) = vc(:,:,1)
4572 ua(:,:,k) = ua(:,:,1)
4573 va(:,:,k) = va(:,:,1)
4576 call mpp_update_domains( uc, vc, domain, gridtype=cgrid_ne_param)
4577 call fill_corners(uc, vc, npx, npy, npz, vector=.true., cgrid=.true.)
5062 real ,
intent(IN) :: p1(2), p2(2)
5063 real ,
intent(IN) :: dist
5064 real ,
intent(IN) :: heading
5065 real ,
intent(OUT) :: p3(2)
5071 p3(2) = asin( (cos(heading)*cos(p1(2))*sin(pha)) + (sin(p1(2))*cos(pha)) )
5072 dp = atan2( sin(heading)*sin(pha)*cos(p1(2)) , cos(pha) - sin(p1(2))*sin(p3(2)) )
5073 p3(1) = mod( (p1(1)-pi)-dp+pi , 2.*pi )
6173 subroutine init_double_periodic(u,v,w,pt,delp,q,phis, ps,pe,peln,pk,pkz, uc,vc, ua,va, ak, bk, &
6174 gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, &
6175 mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd)
6179 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6180 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6181 real ,
intent(INOUT) :: w(bd%isd: ,bd%jsd: ,1:)
6182 real ,
intent(INOUT) :: pt(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6183 real ,
intent(INOUT) :: delp(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6184 real ,
intent(INOUT) :: q(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz, ncnst)
6186 real ,
intent(INOUT) :: phis(bd%isd:bd%ied ,bd%jsd:bd%jed )
6188 real ,
intent(INOUT) :: ps(bd%isd:bd%ied ,bd%jsd:bd%jed )
6189 real ,
intent(INOUT) :: pe(bd%is-1:bd%ie+1,npz+1,bd%js-1:bd%je+1)
6190 real ,
intent(INOUT) :: pk(bd%is:bd%ie ,bd%js:bd%je ,npz+1)
6191 real ,
intent(INOUT) :: peln(bd%is :bd%ie ,npz+1 ,bd%js:bd%je)
6192 real ,
intent(INOUT) :: pkz(bd%is:bd%ie ,bd%js:bd%je ,npz )
6193 real ,
intent(INOUT) :: uc(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
6194 real ,
intent(INOUT) :: vc(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
6195 real ,
intent(INOUT) :: ua(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6196 real ,
intent(INOUT) :: va(bd%isd:bd%ied ,bd%jsd:bd%jed ,npz)
6197 real ,
intent(inout) :: delz(bd%is:,bd%js:,1:)
6198 real ,
intent(inout) :: ze0(bd%is:,bd%js:,1:)
6200 real ,
intent(inout) :: ak(npz+1)
6201 real ,
intent(inout) :: bk(npz+1)
6203 integer,
intent(IN) :: npx, npy, npz
6204 integer,
intent(IN) :: ng, ncnst, nwat
6205 integer,
intent(IN) :: ndims
6206 integer,
intent(IN) :: nregions
6208 real,
intent(IN) :: dry_mass
6209 logical,
intent(IN) :: mountain
6210 logical,
intent(IN) :: moist_phys
6211 logical,
intent(IN) :: hydrostatic, hybrid_z
6212 integer,
intent(INOUT) :: ks
6213 integer,
intent(INOUT),
target :: tile_in
6214 real,
intent(INOUT) :: ptop
6216 type(domain2d),
intent(IN),
target :: domain_in
6221 real,
dimension(bd%is:bd%ie):: pm, qs
6222 real,
dimension(1:npz):: pk1, ts1, qs1
6224 real :: dist, r0, f0_const, prf, rgrav
6225 real :: ptmp, ze, zc, zm, utmp, vtmp
6226 real :: t00, p00, xmax, xc, xx, yy, pk0, pturb, ztop
6230 integer :: i, j, k, m, icenter, jcenter
6232 real,
pointer,
dimension(:,:,:) :: agrid, grid
6233 real(kind=R_GRID),
pointer,
dimension(:,:) :: area
6234 real,
pointer,
dimension(:,:) :: rarea, fC, f0
6235 real,
pointer,
dimension(:,:,:) :: ee1, ee2, en1, en2
6236 real,
pointer,
dimension(:,:,:,:) :: ew, es
6237 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
6239 logical,
pointer :: cubed_sphere, latlon
6241 type(domain2d),
pointer :: domain
6242 integer,
pointer :: tile
6244 logical,
pointer :: have_south_pole, have_north_pole
6246 integer,
pointer :: ntiles_g
6247 real,
pointer :: acapN, acapS, globalarea
6249 real(kind=R_GRID),
pointer :: dx_const, dy_const
6251 integer :: is, ie, js, je
6252 integer :: isd, ied, jsd, jed
6263 agrid => gridstruct%agrid
6264 grid => gridstruct%grid
6266 area => gridstruct%area_64
6270 dxa => gridstruct%dxa
6271 dya => gridstruct%dya
6272 rdxa => gridstruct%rdxa
6273 rdya => gridstruct%rdya
6274 dxc => gridstruct%dxc
6275 dyc => gridstruct%dyc
6281 dx_const => flagstruct%dx_const
6282 dy_const => flagstruct%dy_const
6287 have_south_pole => gridstruct%have_south_pole
6288 have_north_pole => gridstruct%have_north_pole
6290 ntiles_g => gridstruct%ntiles_g
6291 acapn => gridstruct%acapN
6292 acaps => gridstruct%acapS
6293 globalarea => gridstruct%globalarea
6295 f0_const = 2.*omega*sin(flagstruct%deglat/180.*pi)
6316 if (j>0 .and. j<5)
then 6318 if (i>0 .and. i<5)
then 6324 call mpp_update_domains( delp, domain )
6331 r0 = 5.*sqrt(dx_const**2 + dy_const**2)
6336 dist=(i-icenter)*dx_const*(i-icenter)*dx_const &
6337 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6338 dist=min(r0,sqrt(dist))
6339 phis(i,j)=1500.*(1. - (dist/r0))
6360 call hydro_eq(npz, is, ie, js, je, ps, phis, dry_mass, &
6361 delp, ak, bk, pt, delz, area, ng, .false., hydrostatic, hybrid_z, domain)
6365 r0 = 100.*sqrt(dx_const**2 + dy_const**2)
6371 dist = (i-icenter)*dx_const*(i-icenter)*dx_const &
6372 +(j-jcenter)*dy_const*(j-jcenter)*dy_const
6373 dist = min(r0, sqrt(dist))
6375 prf = ak(k) + ps(i,j)*bk(k)
6376 if ( prf > 100.e2 )
then 6377 pt(i,j,k) = pt(i,j,k) + 0.01*(1. - (dist/r0)) * prf/ps(i,j)
6383 if ( hydrostatic )
then 6384 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6385 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
6386 moist_phys, .true., nwat , domain, flagstruct%adiabatic)
6389 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6390 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
6391 moist_phys, hydrostatic, nwat, domain, flagstruct%adiabatic, .true. )
6398 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6401 call qsmith((ie-is+1)*(je-js+1), npz, &
6402 ie-is+1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6404 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6407 q(i,j,k,1) = max(2.e-6, 0.8*pm(i)/ps(i,j)*qs(i) )
6423 if ( .not. hydrostatic ) w(:,:,:) = 0.
6435 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6444 ptmp = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6450 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6451 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
6452 moist_phys, .false., nwat, domain, flagstruct%adiabatic)
6455 r0 = 5.*max(dx_const, dy_const)
6464 zm = ze - 0.5*delz(i,j,k)
6465 ze = ze - delz(i,j,k)
6466 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + &
6469 if ( dist <= r0 )
then 6470 pt(i,j,k) = pt(i,j,k) + 5.*(1.-dist/r0)
6493 ze1(k) = ze1(k+1) + ztop/
real(npz)
6507 delz(i,j,k) = ze1(k+1) - ze1(k)
6508 pk(i,j,k) = pk(i,j,k+1) + grav*delz(i,j,k)/(cp_air*t00)*pk0
6509 pe(i,k,j) = pk(i,j,k)**(1./kappa)
6515 if ( is_master() )
write(*,*)
'Density curent testcase: model top (mb)=', ptop/100.
6520 peln(i,k,j) = log(pe(i,k,j))
6529 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6530 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6541 zm = (0.5*(ze1(k)+ze1(k+1))-3.e3) / 2.e3
6545 xx = (dx_const * (0.5+
real(i-1)) - xc) / 4.e3
6546 yy = (dy_const * (0.5+
real(j-1)) - xc) / 4.e3
6547 dist = sqrt( xx**2 + yy**2 + zm**2 )
6548 if ( dist<=1. )
then 6549 pt(i,j,k) = pt(i,j,k) - pturb/pkz(i,j,k)*(cos(pi*dist)+1.)/2.
6552 pt(i,j,k) = pt(i,j,k) * pkz(i,j,k)
6561 zvir = rvgas/rdgas - 1.
6567 pk(i,j,1) = ptop**kappa
6569 peln(i,1,j) = log(ptop)
6576 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6577 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6578 peln(i,k+1,j) = log(pe(i,k+1,j))
6579 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
6587 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6603 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6610 ze1(k) = ze1(k+1) - delz(is,js,k)
6614 zm = 0.5*(ze1(k)+ze1(k+1))
6615 utmp = us0*tanh(zm/3.e3)
6623 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6624 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
6625 .true., hydrostatic, nwat, domain, flagstruct%adiabatic)
6631 icenter = (npx-1)/3 + 1
6632 jcenter = (npy-1)/2 + 1
6634 zm = 0.5*(ze1(k)+ze1(k+1))
6635 ptmp = ( (zm-zc)/zc ) **2
6636 if ( ptmp < 1. )
then 6639 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6640 if ( dist < 1. )
then 6641 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6653 zvir = rvgas/rdgas - 1.
6659 pk(i,j,1) = ptop**kappa
6661 peln(i,1,j) = log(ptop)
6668 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
6669 pe(i,k+1,j) = ak(k+1) + ps(i,j)*bk(k+1)
6670 peln(i,k+1,j) = log(pe(i,k+1,j))
6671 pk(i,j,k+1) = exp( kappa*peln(i,k+1,j) )
6679 pk1(k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6693 delz(i,j,k) = rdgas/grav*ts1(k)*(1.+zvir*qs1(k))*(peln(i,k,j)-peln(i,k+1,j))
6700 ze1(k) = ze1(k+1) - delz(is,js,k)
6706 zm = 0.5*(ze1(k)+ze1(k+1))
6707 if ( zm .le. 2.e3 )
then 6708 utmp = 8.*(1.-cos(pi*zm/4.e3))
6709 vtmp = 8.*sin(pi*zm/4.e3)
6710 elseif (zm .le. 6.e3 )
then 6711 utmp = 8. + (us0-8.)*(zm-2.e3)/4.e3
6720 u(i,j,k) = utmp - 8.
6726 v(i,j,k) = vtmp - 4.
6732 call p_var(npz, is, ie, js, je, ptop,
ptop_min, delp, delz, pt, ps, &
6733 pe, peln, pk, pkz, kappa, q, ng, ncnst, area, dry_mass, .false., .false., &
6734 .true., hydrostatic, nwat, domain, flagstruct%adiabatic)
6740 icenter = (npx-1)/2 + 1
6741 jcenter = (npy-1)/2 + 1
6743 zm = 0.5*(ze1(k)+ze1(k+1))
6744 ptmp = ( (zm-zc)/zc ) **2
6745 if ( ptmp < 1. )
then 6748 dist = ptmp+((i-icenter)*dx_const/r0)**2+((j-jcenter)*dy_const/r0)**2
6749 if ( dist < 1. )
then 6750 pt(i,j,k) = pt(i,j,k) + pturb*(1.-sqrt(dist))
6772 if (.not.hybrid_z)
call mpp_error(fatal,
'hybrid_z must be .TRUE.')
6777 call mpp_error(fatal,
'npz must be == 101 ')
6782 call set_hybrid_z(is, ie, js, je, ng, npz, ztop, dz1, rgrav, &
6790 peln(i,npz+1,j) = log(p00)
6797 peln(i,k,j) = peln(i,k+1,j) + grav*delz(i,j,k)/(rdgas*t00)
6798 pe(i,k,j) = exp(peln(i,k,j))
6799 pk(i,j,k) = pe(i,k,j)**kappa
6808 if ( is_master() )
write(*,*)
'LES testcase: computed model top (mb)=', ptop/100.
6813 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
6814 delp(i,j,k) = pe(i,k+1,j)-pe(i,k,j)
6822 pm(i) = delp(i,j,k)/(peln(i,k+1,j)-peln(i,k,j))
6825 call qsmith((ie-is+1)*(je-js+1), npz, &
6826 ie-is+1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6828 call qsmith(ie-is+1, 1, 1, pt(is:ie,j,k), pm, q(is:ie,j,k,1), qs)
6831 if ( pm(i) > 100.e2 )
then 6832 q(i,j,k,1) = 0.9*qs(i)
6849 zm = 0.5*(ze0(i,j,k)+ze0(i,j,k+1))
6850 dist = ((i-icenter)*dx_const)**2 + ((j-jcenter)*dy_const)**2 + (zm-zc)**2
6852 if ( dist <= r0 )
then 6853 pt(i,j,k) = pt(i,j,k) + 2.0*(1.-dist/r0)
6891 nullify(have_south_pole)
6892 nullify(have_north_pole)
6903 character(*),
intent(IN) :: nml_filename
6904 integer :: ierr, f_unit, unit, ios
6906 #include<file_version.h> 6916 #ifdef INTERNAL_FILE_NML 6918 read (input_nml_file,test_case_nml,iostat=ios)
6919 ierr = check_nml_error(ios,
'test_case_nml')
6921 f_unit = open_namelist_file(nml_filename)
6925 read (f_unit,test_case_nml,iostat=ios)
6926 ierr = check_nml_error(ios,
'test_case_nml')
6927 call close_file(f_unit)
6929 write(unit, nml=test_case_nml)
6938 integer,
intent(in):: km
6939 real,
intent(in):: p00
6940 real,
intent(inout),
dimension(km+1):: pe
6941 real,
intent(in),
dimension(km+1):: ze
6944 real,
intent(out),
dimension(km):: pt, qz
6946 integer,
parameter:: nx = 5
6947 real,
parameter:: qst = 1.0e-6
6948 real,
parameter:: qv0 = 1.4e-2
6949 real,
parameter:: ztr = 12.e3
6950 real,
parameter:: ttr = 213.
6951 real,
parameter:: ptr = 343.
6952 real,
parameter:: pt0 = 300.
6953 real,
dimension(km):: zs, rh, temp, dp, dp0
6954 real,
dimension(km+1):: peln, pk
6955 real:: qs, zvir, fac_z, pk0, temp1, pm
6958 zvir = rvgas/rdgas - 1.
6960 if ( (is_master()) )
then 6961 write(*,*)
'Computing sounding for HIWPP super-cell test using p00=', p00
6968 zs(k) = 0.5*(ze(k)+ze(k+1))
6970 if ( zs(k) .gt. ztr )
then 6972 pt(k) = ptr*exp(grav*(zs(k)-ztr)/(cp_air*ttr))
6975 fac_z = (zs(k)/ztr)**1.25
6976 pt(k) = pt0 + (ptr-pt0)* fac_z
6977 rh(k) = 1. - 0.75 * fac_z
6979 qz(k) = qv0 - (qv0-qst)*fac_z
6981 if ( is_master() )
write(*,*) zs(k), pt(k), qz(k)
6986 #ifdef USE_MOIST_P00 6993 peln(km+1) = log(p00)
6998 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k)*(1.+zvir*qz(k)))
6999 peln(k) = log(pk(k)) / kappa
7000 pe(k) = exp(peln(k))
7003 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
7004 temp(k) = pt(k)*pm**kappa
7006 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
7007 qz(k) = min( qv0, rh(k)*qs )
7008 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
7015 peln(km+1) = log(p00)
7019 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k))
7020 peln(k) = log(pk(k)) / kappa
7021 pe(k) = exp(peln(k))
7024 dp0(k) = pe(k+1) - pe(k)
7025 pm = dp0(k)/(peln(k+1)-peln(k))
7026 temp(k) = pt(k)*pm**kappa
7028 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
7029 qz(k) = min( qv0, rh(k)*qs )
7035 dp(k) = dp0(k)*(1. + qz(k))
7036 pe(k+1) = pe(k) + dp(k)
7039 pk(km+1) = pe(km+1)**kappa
7040 peln(km+1) = log(pe(km+1))
7044 pk(k) = pk(k+1) - grav*(ze(k)-ze(k+1))/(cp_air*pt(k)*(1.+zvir*qz(k)))
7045 peln(k) = log(pk(k)) / kappa
7046 pe(k) = exp(peln(k))
7049 pm = (pe(k+1)-pe(k))/(peln(k+1)-peln(k))
7050 temp(k) = pt(k)*pm**kappa
7052 qs = 380./pm*exp(17.27*(temp(k)-273.)/(temp(k)-36.))
7053 qz(k) = min( qv0, rh(k)*qs )
7054 if ( n==nx .and. is_master() )
write(*,*) 0.01*pm, temp(k), qz(k), qs
7059 if ( is_master() )
then 7060 write(*,*)
'Super_K: computed ptop (mb)=', 0.01*pe(1),
' PS=', 0.01*pe(km+1)
7061 call prt_m1(
'1D Sounding T0', temp, 1, km, 1, 1, 0, 1, 1.)
7066 subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, &
7067 delz, zvir, ptop, ak, bk, agrid)
7068 integer,
intent(in):: is, ie, js, je, ng, km
7069 real,
intent(in),
dimension(km ):: ts1, qs1, uz1, dudz
7070 real,
intent(in),
dimension(km+1):: ze1
7071 real,
intent(in):: zvir, ps0
7072 real,
intent(inout):: ptop
7073 real(kind=R_GRID),
intent(in):: agrid(is-ng:ie+ng,js-ng:je+ng,2)
7074 real,
intent(inout),
dimension(km+1):: ak, bk
7075 real,
intent(inout),
dimension(is:ie,js:je,km):: pt
7076 real,
intent(inout),
dimension(is:,js:,1:) :: delz
7077 real,
intent(out),
dimension(is:ie,js:je,km+1):: pk
7079 real,
intent(inout),
dimension(is-1:ie+1,km+1,js-1:je+1):: pe
7081 integer,
parameter:: nt=5
7082 integer,
parameter:: nlat=1001
7083 real,
dimension(nlat,km):: pt2, pky, dzc
7084 real,
dimension(nlat,km+1):: pk2, pe2, peln2, pte
7085 real,
dimension(km+1):: pe1
7086 real:: lat(nlat), latc(nlat-1)
7087 real:: fac_y, dlat, dz0, pk0, tmp1, tmp2, tmp3, pint
7088 integer::i,j,k,n, jj, k1
7092 dz0 = ze1(km) - ze1(km+1)
7095 dlat = 0.5*pi/
real(nlat-1)
7097 lat(j) = dlat*
real(j-1)
7099 dzc(j,k) = ze1(k) - ze1(k+1)
7103 latc(j) = 0.5*(lat(j)+lat(j+1))
7112 if ( is_master() )
then 7114 call prt_m1(
'Super_K PT0', pt2, 1, nlat, 1, km, 0, 1, tmp1)
7121 call ppme(pt2, pte, dzc, nlat, km)
7124 tmp1 = 0.5*(pte(j-1,k ) + pte(j,k ))
7125 tmp3 = 0.5*(pte(j-1,k+1) + pte(j,k+1))
7126 pt2(j,k) = pt2(j-1,k) + dlat/(2.*grav)*sin(2.*latc(j-1))*uz1(k)* &
7127 ( uz1(k)*(tmp1-tmp3)/dzc(j,k) - (pt2(j-1,k)+pt2(j,k))*dudz(k) )
7130 if ( is_master() )
then 7131 call prt_m1(
'Super_K PT', pt2, 1, nlat, 1, km, 0, 1, pk0/cp_air)
7137 pk2(1,km+1) = ps0**kappa
7139 pk2(j,km+1) = pk2(j-1,km+1) - dlat*uz1(km)*uz1(km)*sin(2.*latc(j-1)) &
7140 / (pt2(j-1,km) + pt2(j,km))
7145 pk2(j,k) = pk2(j,k+1) - grav*dzc(j,k)/pt2(j,k)
7151 peln2(j,k) = log(pk2(j,k)) / kappa
7152 pe2(j,k) = exp(peln2(j,k))
7158 pky(j,k) = (pk2(j,k+1)-pk2(j,k))/(kappa*(peln2(j,k+1)-peln2(j,k)))
7159 pt2(j,k) = pt2(j,k)*pky(j,k)/(cp_air*(1.+zvir*qs1(k)))
7167 if ( is_master() )
then 7168 write(*,*)
'SuperK ptop at EQ=', 0.01*pe1(1),
'new ptop=', 0.01*ptop
7169 call prt_m1(
'Super_K pe', pe2, 1, nlat, 1, km+1, 0, 1, 0.01)
7170 call prt_m1(
'Super_K Temp', pt2, 1, nlat, 1, km, 0, 1, 1.)
7177 if (abs(agrid(i,j,2))>=lat(jj) .and. abs(agrid(i,j,2))<=lat(jj+1) )
then 7179 fac_y = (abs(agrid(i,j,2))-lat(jj)) / dlat
7181 pt(i, j,k) = pt2(jj, k) + fac_y*(pt2(jj+1, k)-pt2(jj,k))
7184 pe(i,k,j) = pe2(jj,k) + fac_y*(pe2(jj+1,k)-pe2(jj,k))
7207 bk(k) = (pe1(k) - pint) / (pe1(km+1)-pint)
7208 ak(k) = pe1(k) - bk(k) * pe1(km+1)
7209 if ( is_master() )
write(*,*) k, ak(k), bk(k)
7222 subroutine superk_u(km, zz, um, dudz)
7223 integer,
intent(in):: km
7224 real,
intent(in):: zz(km)
7225 real,
intent(out):: um(km), dudz(km)
7227 real,
parameter:: zs = 5.e3
7228 real,
parameter:: us = 30.
7235 if ( zz(k) .gt. zs+1.e3 )
then 7238 elseif ( abs(zz(k)-zs) .le. 1.e3 )
then 7239 um(k) = us*(-4./5. + 3.*zz(k)/zs - 5./4.*(zz(k)/zs)**2)
7240 dudz(k) = us/zs*(3. - 5./2.*zz(k)/zs)
7249 um(k) = us*tanh( zz(k)/zs ) - uc
7250 dudz(k) = (us/zs)/cosh(zz(k)/zs)**2
7258 use gfdl_cloud_microphys_mod
, only: wqsat_moist, qsmith_init, qs_blend
7261 integer,
intent(in):: km
7262 real,
intent(in):: ps
7263 real,
intent(in),
dimension(km):: pk1
7264 real,
intent(out),
dimension(km):: tp, qp
7266 integer,
parameter:: ns = 401
7267 integer,
parameter:: nx = 3
7268 real,
dimension(ns):: zs, pt, qs, us, rh, pp, pk, dpk, dqdt
7269 real,
parameter:: Tmin = 175.
7270 real,
parameter:: p00 = 1.0e5
7271 real,
parameter:: qst = 3.0e-6
7272 real,
parameter:: qv0 = 1.4e-2
7273 real,
parameter:: ztr = 12.e3
7274 real,
parameter:: ttr = 213.
7275 real,
parameter:: ptr = 343.
7276 real,
parameter:: pt0 = 300.
7277 real:: dz0, zvir, fac_z, pk0, temp1, p2
7286 zvir = rvgas/rdgas - 1.
7290 if ( (is_master()) )
then 7291 write(*,*)
'Computing sounding for super-cell test' 7302 zs(k) = zs(k+1) + dz0
7307 if ( zs(k) .gt. ztr )
then 7309 pt(k) = ptr*exp(grav*(zs(k)-ztr)/(cp_air*ttr))
7312 fac_z = (zs(k)/ztr)**1.25
7313 pt(k) = pt0 + (ptr-pt0)* fac_z
7314 rh(k) = 1. - 0.75 * fac_z
7316 qs(k) = qv0 - (qv0-qst)*fac_z
7326 temp1 = 0.5*(pt(k)*(1.+zvir*qs(k)) + pt(k+1)*(1.+zvir*qs(k+1)))
7327 dpk(k) = grav*(zs(k)-zs(k+1))/(cp_air*temp1)
7331 pk(k) = pk(k+1) - dpk(k)
7337 if ( pk(k) > 0. )
then 7338 pp(k) = exp(log(pk(k))/kappa)
7340 qs(k) = 380./pp(k)*exp(17.27*(temp1-273.)/(temp1-36.))
7341 qs(k) = min( qv0, rh(k)*qs(k) )
7342 if ( (is_master()) )
write(*,*) 0.01*pp(k), qs(k)
7345 #ifdef USE_MIXED_TABLE 7346 qs(k) = min(qv0, rh(k)*qs_blend(temp1, pp(k), qs(k)))
7348 qs(k) = min(qv0, rh(k)*wqsat_moist(temp1, qs(k), pp(k)))
7353 if ( (is_master()) )
write(*,*) n, k, pk(k)
7354 call mpp_error(fatal,
'Super-Cell case: pk < 0')
7361 if ( pk1(k) .le. pk(1) )
then 7362 tp(k) = pt(1)*pk(1)/pk1(k)
7364 elseif ( pk1(k) .ge. pk(ns) )
then 7369 if( (pk1(k).le.pk(kk+1)) .and. (pk1(k).ge.pk(kk)) )
then 7370 fac_z = (pk1(k)-pk(kk))/(pk(kk+1)-pk(kk))
7371 tp(k) = pt(kk) + (pt(kk+1)-pt(kk))*fac_z
7372 qp(k) = qs(kk) + (qs(kk+1)-qs(kk))*fac_z
7380 tp(k) = tp(k)*pk1(k)
7381 tp(k) = max(tmin, tp(k))
7390 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7391 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7392 hydrostatic, nwat, adiabatic, do_pert, domain, bd)
7396 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7397 real,
intent(IN) :: ptop
7398 real,
intent(IN),
dimension(npz+1) :: ak, bk
7399 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7400 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w
7401 real,
intent(OUT),
dimension(is:,js:,1:) :: delz
7402 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7403 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7404 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7405 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7406 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7407 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7408 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7409 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7410 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7411 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7412 logical,
intent(IN) :: hydrostatic,adiabatic,do_pert
7413 type(domain2d),
intent(INOUT) :: domain
7415 real,
parameter :: p0 = 1.e5
7416 real,
parameter :: u0 = 35.
7417 real,
parameter :: b = 2.
7418 real,
parameter :: KK = 3.
7419 real,
parameter :: Te = 310.
7420 real,
parameter :: Tp = 240.
7421 real,
parameter :: T0 = 0.5*(te + tp)
7422 real,
parameter :: up = 1.
7423 real,
parameter :: zp = 1.5e4
7424 real(kind=R_GRID),
parameter :: lamp = pi/9.
7425 real(kind=R_GRID),
parameter :: phip = 2.*lamp
7426 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7427 real,
parameter :: Rp =
radius/10.
7428 real,
parameter :: lapse = 5.e-3
7429 real,
parameter :: dT = 4.8e5
7430 real,
parameter :: phiW = 2.*pi/9.
7431 real,
parameter :: pW = 34000.
7432 real,
parameter :: q0 = .018
7433 real,
parameter :: qt = 1.e-12
7434 real,
parameter :: ptrop = 1.e4
7436 real,
parameter :: zconv = 1.e-6
7437 real,
parameter :: rdgrav = rdgas/grav
7440 real,
parameter :: rrdgrav = grav/rdgas
7442 integer :: i,j,k,iter, sphum, cl, cl2, n
7443 real :: p,z,z0,ziter,piter,titer,uu,vv,pl,pt_u,pt_v
7444 real(kind=R_GRID),
dimension(2) :: pa
7445 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7446 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2
7447 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7448 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2
7449 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7456 zvir = rvgas/rdgas - 1.
7469 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7480 peln(i,1,j) = log(ptop)
7481 pk(i,j,1) = ptop**kappa
7485 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7488 pk(i,j,k) = exp(kappa*log(pe(i,k,j)))
7489 peln(i,k,j) = log(pe(i,k,j))
7497 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
7519 z = ziter + (piter - p)*rdgrav*titer/piter
7525 if (abs(z - ziter) < zconv)
exit 7536 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7541 call mpp_update_domains(pt, domain)
7542 call mpp_update_domains(gz, domain)
7549 peln_u(i,j) = log(p0)
7564 p = ak(k) + ps_u(i,j)*bk(k)
7573 z = ziter + (piter - p)*rdgrav*titer/piter
7574 if (abs(z - ziter) < zconv)
exit 7577 pt_u = rrdgrav * ( z - gz_u(i,j) ) / (peln_u(i,j) - pl)
7584 u(i,j,k) = u1(i,j)*uu
7597 peln_v(i,j) = log(p0)
7612 p = ak(k) + ps_v(i,j)*bk(k)
7621 z = ziter + (piter - p)*rdgrav*titer/piter
7622 if (abs(z - ziter) < zconv)
exit 7625 pt_v = rrdgrav * ( z - gz_v(i,j) ) / (peln_v(i,j) - pl)
7631 v(i,j,k) = v1(i,j)*uu
7640 if (.not. hydrostatic)
then 7646 delz(i,j,k) = rdgrav * (peln(i,k+1,j) - peln(i,k,j)) * pt(i,j,k)
7662 sphum = get_tracer_index(model_atmos,
'sphum')
7666 p = delp(i,j,k)/(peln(i,k+1,j) - peln(i,k,j))
7667 q(i,j,k,sphum) =
dcmip16_bc_sphum(p,ps(i,j),agrid(i,j,2),agrid(i,j,1))
7672 cl = get_tracer_index(model_atmos,
'cl')
7673 cl2 = get_tracer_index(model_atmos,
'cl2')
7674 if (cl > 0 .and. cl2 > 0)
then 7676 q, delp,nq,agrid(isd,jsd,1),agrid(isd,jsd,2),bd)
7677 call mpp_update_domains(q,domain)
7680 if (.not. adiabatic)
then 7685 pt(i,j,k) = pt(i,j,k) / ( 1. + zvir*q(i,j,k,sphum))
7696 real,
intent(IN) :: z
7697 real(kind=R_GRID),
intent(IN) :: lat
7698 real :: IT, T1, T2, Tr, zsc
7700 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7701 zsc = z*grav/(b*rdgas*t0)
7702 tr = ( 1. - 2.*zsc**2.) * exp(-zsc**2. )
7704 t1 = (1./t0)*exp(lapse*z/t0) + (t0 - tp)/(t0*tp) * tr
7705 t2 = 0.5* ( kk + 2.) * (te - tp)/(te*tp) * tr
7713 real,
intent(IN) :: z
7714 real(kind=R_GRID),
intent(IN) :: lat
7715 real :: IT, Ti1, Ti2, Tir
7717 it = exp(kk * log(cos(lat))) - kk/(kk+2.)*exp((kk+2.)*log(cos(lat)))
7718 tir = z*exp(-(z*grav/(b*rdgas*t0))*(z*grav/(b*rdgas*t0)) )
7720 ti1 = 1./lapse* (exp(lapse*z/t0) - 1.) + tir*(t0-tp)/(t0*tp)
7721 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7729 real,
intent(IN) :: z, T
7730 real(kind=R_GRID),
intent(IN) :: lat
7731 real :: Tir, Ti2, UU, ur
7733 tir = z*exp(-(z*grav/(b*rdgas*t0))*(z*grav/(b*rdgas*t0)) )
7734 ti2 = 0.5*(kk+2.)*(te-tp)/(te*tp) * tir
7736 uu = grav*kk/
radius * ti2 * ( cos(lat)**(int(kk)-1) - cos(lat)**(int(kk)+1) ) * t
7737 ur = - omega *
radius * cos(lat) + sqrt( (omega*
radius*cos(lat))**2 +
radius*cos(lat)*uu)
7745 real,
intent(IN) :: z
7746 real(kind=R_GRID),
intent(IN) :: lat, lon
7748 real(kind=R_GRID) :: dst, pphere(2)
7751 zz = max(1. - 3.*zrat*zrat + 2.*zrat*zrat*zrat, 0.)
7753 pphere = (/ lon, lat /)
7762 real,
intent(IN) :: p, ps
7763 real(kind=R_GRID),
intent(IN) :: lat, lon
7778 is,ie,js,je,isd,ied,jsd,jed,npz,nq,ak,bk,ptop, &
7779 pk,peln,pe,pkz,gz,phis,ps,grid,agrid, &
7780 hydrostatic, nwat, adiabatic)
7782 integer,
intent(IN) :: is,ie,js,je,isd,ied,jsd,jed,npz,nq, nwat
7783 real,
intent(IN) :: ptop
7784 real,
intent(IN),
dimension(npz+1) :: ak, bk
7785 real,
intent(INOUT),
dimension(isd:ied,jsd:jed,npz,nq) :: q
7786 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz) :: delp, pt, w
7787 real,
intent(OUT),
dimension(is:,js:,1:) :: delz
7788 real,
intent(OUT),
dimension(isd:ied,jsd:jed+1,npz) :: u
7789 real,
intent(OUT),
dimension(isd:ied+1,jsd:jed,npz) :: v
7790 real,
intent(OUT),
dimension(is:ie,js:je,npz+1) :: pk
7791 real,
intent(OUT),
dimension(is:ie,npz+1,js:je) :: peln
7792 real,
intent(OUT),
dimension(is-1:ie+1,npz+1,js-1:je+1) :: pe
7793 real,
intent(OUT),
dimension(is:ie,js:je,npz) :: pkz
7794 real,
intent(OUT),
dimension(isd:ied,jsd:jed) :: phis,ps
7795 real(kind=R_GRID),
intent(IN),
dimension(isd:ied,jsd:jed,2) :: agrid
7796 real(kind=R_GRID),
intent(IN),
dimension(isd:ied+1,jsd:jed+1,2) :: grid
7797 real,
intent(OUT),
dimension(isd:ied,jsd:jed,npz+1) :: gz
7798 logical,
intent(IN) :: hydrostatic,adiabatic
7800 real,
parameter :: zt = 15000
7801 real,
parameter :: q0 = 0.021
7802 real,
parameter :: qt = 1.e-11
7803 real,
parameter :: T0 = 302.15
7804 real,
parameter :: Tv0 = 302.15*(1.+0.608*q0)
7805 real,
parameter :: Ts = 302.15
7806 real,
parameter :: zq1 = 3000.
7807 real,
parameter :: zq2 = 8000.
7808 real,
parameter :: lapse = 7.e-3
7809 real,
parameter :: Tvt = tv0 - lapse*zt
7810 real,
parameter :: pb = 101500.
7811 real,
parameter :: ptt = pb*(tvt/tv0)**(grav/rdgas/lapse)
7812 real(kind=R_GRID),
parameter :: lamp = pi
7813 real(kind=R_GRID),
parameter :: phip = pi/18.
7814 real(kind=R_GRID),
parameter :: ppcenter(2) = (/ lamp, phip /)
7815 real,
parameter :: dp = 1115.
7816 real,
parameter :: rp = 282000.
7817 real,
parameter :: zp = 7000.
7818 real,
parameter :: fc = 2.*omega*sin(phip)
7820 real,
parameter :: zconv = 1.e-6
7821 real,
parameter :: rdgrav = rdgas/grav
7822 real,
parameter :: rrdgrav = grav/rdgas
7823 real,
parameter :: zvir = rvgas/rdgas - 1.
7825 integer :: i,j,k,iter, sphum, cl, cl2, n
7826 real :: p,z,z0,ziter,piter,titer,uu,vv,pl, r
7827 real(kind=R_GRID),
dimension(2) :: pa
7828 real(kind=R_GRID),
dimension(3) :: e1,e2,ex,ey
7829 real,
dimension(is:ie,js:je) :: rc
7830 real,
dimension(is:ie,js:je+1) :: gz_u,p_u,peln_u,ps_u,u1,u2, rc_u
7831 real(kind=R_GRID),
dimension(is:ie,js:je+1) :: lat_u,lon_u
7832 real,
dimension(is:ie+1,js:je) :: gz_v,p_v,peln_v,ps_v,v1,v2, rc_v
7833 real(kind=R_GRID),
dimension(is:ie+1,js:je) :: lat_v,lon_v
7851 ps(i,j) = pb - dp*exp( -sqrt((rc(i,j)/rp)**3) )
7859 delp(i,j,k) = ak(k+1)-ak(k) + ps(i,j)*(bk(k+1)-bk(k))
7870 peln(i,1,j) = log(ptop)
7871 pk(i,j,1) = ptop**kappa
7875 pe(i,k,j) = ak(k) + ps(i,j)*bk(k)
7878 pk(i,j,k) = exp(kappa*log(pe(i,k,j)))
7879 peln(i,k,j) = log(pe(i,k,j))
7887 pkz(i,j,k) = (pk(i,j,k+1)-pk(i,j,k))/(kappa*(peln(i,k+1,j)-peln(i,k,j)))
7909 z = ziter + (piter - p)*rdgrav*titer/piter
7915 if (abs(z - ziter) < zconv)
exit 7926 pt(i,j,k) = rrdgrav * ( gz(i,j,k) - gz(i,j,k+1) ) / ( peln(i,k+1,j) - peln(i,k,j))
7944 p_u(i,j) = pb - dp*exp( -sqrt((rc_u(i,j)/rp)**3) )
7945 peln_u(i,j) = log(p_u(i,j))
7946 ps_u(i,j) = p_u(i,j)
7953 p = ak(k) + ps_u(i,j)*bk(k)
7962 z = ziter + (piter - p)*rdgrav*titer/piter
7963 if (abs(z - ziter) < zconv)
exit 7967 u(i,j,k) = u1(i,j)*uu + u2(i,j)*vv
7987 p_v(i,j) = pb - dp*exp( - sqrt((rc_v(i,j)/rp)**3) )
7988 peln_v(i,j) = log(p_v(i,j))
7989 ps_v(i,j) = p_v(i,j)
7996 p = ak(k) + ps_v(i,j)*bk(k)
8005 z = ziter + (piter - p)*rdgrav*titer/piter
8006 if (abs(z - ziter) < zconv)
exit 8010 v(i,j,k) = v1(i,j)*uu + v2(i,j)*vv
8028 if (.not. adiabatic)
then 8029 sphum = get_tracer_index(model_atmos,
'sphum')
8033 z = 0.5*(gz(i,j,k) + gz(i,j,k+1))
8036 pt(i,j,k) = pt(i,j,k) / ( 1. + zvir*q(i,j,k,sphum))
8043 if (.not. hydrostatic)
then 8048 delz(i,j,k) = gz(i,j,k) - gz(i,j,k+1)
8059 real,
intent(IN) :: z, r
8060 real :: Tv, term1, term2
8068 term1 = grav*zp*zp* ( 1. - pb/dp * exp( sqrt(r/rp)**3 + (z/zp)**2 ) )
8069 term2 = 2*rdgas*tv*z
8077 real,
intent(IN) :: z, r
8080 dcmip16_tc_pressure = pb*exp(grav/(rdgas*lapse) * log( (tv0-lapse*z)/tv0) ) -dp* exp(-sqrt((r/rp)**3) - (z/zp)**2) * &
8081 exp( grav/(rdgas*lapse) * log( (tv0-lapse*z)/tv0) )
8090 real,
intent(IN) :: z, r
8091 real(kind=R_GRID),
intent(IN) :: lon, lat
8092 real,
intent(OUT) :: uu, vv
8093 real :: rfac, Tvrd, vt, fr5, d1, d2, d
8094 real(kind=R_GRID) :: dst, pphere(2)
8102 rfac = sqrt(r/rp)**3
8105 tvrd = (tv0 - lapse*z)*rdgas
8107 vt = -fr5 + sqrt( fr5**2 - (1.5 * rfac * tvrd) / &
8108 ( 1. + 2*tvrd*z/(grav*zp**2) - pb/dp*exp( rfac + (z/zp)**2) ) )
8110 d1 = sin(phip)*cos(lat) - cos(phip)*sin(lat)*cos(lon - lamp)
8111 d2 = cos(phip)*sin(lon - lamp)
8112 d = max(1.e-25,sqrt(d1*d1 + d2*d2))
8121 real,
intent(IN) :: z
8753 subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, bounded_domain, domain, bd)
8756 integer,
intent(IN) :: npx, npy, ng
8757 real ,
intent(IN) :: uin(bd%isd:bd%ied ,bd%jsd:bd%jed )
8758 real ,
intent(IN) :: vin(bd%isd:bd%ied ,bd%jsd:bd%jed )
8759 real ,
intent(OUT) :: uout(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
8760 real ,
intent(OUT) :: vout(bd%isd:bd%ied+1,bd%jsd:bd%jed )
8761 logical,
intent(IN) :: bounded_domain
8762 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed) :: dxa, dya
8763 real ,
intent(IN),
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed) :: dxc
8764 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1) :: dyc
8765 type(domain2d),
intent(INOUT) :: domain
8769 real :: tmp1i(bd%isd:bd%ied+1)
8770 real :: tmp2i(bd%isd:bd%ied)
8771 real :: tmp3i(bd%isd:bd%ied)
8772 real :: tmp1j(bd%jsd:bd%jed+1)
8773 real :: tmp2j(bd%jsd:bd%jed)
8774 real :: tmp3j(bd%jsd:bd%jed)
8776 integer :: jsd, jed, isd, ied
8784 tmp2i(:) = vin(:,j)*dxa(:,j)
8787 vout(:,j) = tmp1i(:)/dxc(:,j)
8791 tmp2j(:) = uin(i,:)*dya(i,:)
8794 uout(i,:) = tmp1j(:)/dyc(i,:)
8797 if (.not. bounded_domain)
call fill_corners(uout, vout, npx, npy, vector=.true., dgrid=.true.)
8808 subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng, bd)
8811 integer,
intent(IN) :: npx, npy, ng
8812 real ,
intent(IN) :: uin(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
8813 real ,
intent(IN) :: vin(bd%isd:bd%ied+1,bd%jsd:bd%jed )
8814 real ,
intent(OUT) :: uout(bd%isd:bd%ied ,bd%jsd:bd%jed )
8815 real ,
intent(OUT) :: vout(bd%isd:bd%ied ,bd%jsd:bd%jed )
8816 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1) :: dx, dyc
8817 real ,
intent(IN),
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed) :: dy, dxc
8818 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed) :: dxa, dya
8822 real :: tmp1i(bd%isd:bd%ied+1)
8823 real :: tmp2i(bd%isd:bd%ied+1)
8824 real :: tmp3i(bd%isd:bd%ied+1)
8825 real :: tmp1j(bd%jsd:bd%jed+1)
8826 real :: tmp2j(bd%jsd:bd%jed+1)
8827 real :: tmp3j(bd%jsd:bd%jed+1)
8829 integer :: is, ie, js, je
8830 integer :: isd, ied, jsd, jed
8846 uout(i,j) = 0.5*(uin(i,j)*dx(i,j)+uin(i,j+1)*dx(i,j+1))/dxa(i,j)
8847 vout(i,j) = 0.5*(vin(i,j)*dy(i,j)+vin(i+1,j)*dy(i+1,j))/dya(i,j)
8853 tmp2j(:) = uin(i,:)*dyc(i,:)
8856 uout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dya(i,jsd:jed)
8860 tmp2i(:) = vin(:,j)*dxc(:,j)
8863 vout(isd:ied,j) = tmp1i(isd+1:ied+1)/dxa(isd:ied,j)
8877 subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, bounded_domain, domain, bd, noComm)
8880 integer,
intent(IN) :: npx, npy, ng
8881 real ,
intent(IN) :: uin(bd%isd:bd%ied ,bd%jsd:bd%jed )
8882 real ,
intent(IN) :: vin(bd%isd:bd%ied ,bd%jsd:bd%jed )
8883 real ,
intent(OUT) :: uout(bd%isd:bd%ied+1,bd%jsd:bd%jed )
8884 real ,
intent(OUT) :: vout(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
8885 logical,
intent(IN) :: bounded_domain
8886 logical,
OPTIONAL,
intent(IN) :: noComm
8887 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1) :: dx
8888 real ,
intent(IN),
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed) :: dy
8889 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed) :: dxa, dya
8890 type(domain2d),
intent(INOUT) :: domain
8895 real :: tmp1i(bd%isd:bd%ied+1)
8896 real :: tmp2i(bd%isd:bd%ied)
8897 real :: tmp3i(bd%isd:bd%ied)
8898 real :: tmp1j(bd%jsd:bd%jed+1)
8899 real :: tmp2j(bd%jsd:bd%jed)
8900 real :: tmp3j(bd%jsd:bd%jed)
8902 integer :: is, ie, js, je
8903 integer :: isd, ied, jsd, jed
8915 #if !defined(ALT_INTERP) 8920 uout(i,j) = ( uin(i,j)*dxa(i,j) + uin(i-1,j)*dxa(i-1,j) ) &
8921 / ( dxa(i,j) + dxa(i-1,j) )
8926 vout(i,j) = ( vin(i,j)*dya(i,j) + vin(i,j-1)*dya(i,j-1) ) &
8927 / ( dya(i,j) + dya(i,j-1) )
8939 vout(i,:) = tmp1j(:)
8946 tmp2i(:) = uin(:,j)*dya(:,j)
8949 uout(:,j) = tmp1i(:)/dy(:,j)
8953 tmp2j(:) = vin(i,:)*dxa(i,:)
8956 vout(i,:) = tmp1j(:)/dx(i,:)
8959 if (cubed_sphere .and. .not. bounded_domain)
then 8960 csfac = cos(30.0*pi/180.0)
8962 if ( (is==1) .and. (js==1) )
then 8965 uout(i,j)=uout(i,j)*csfac
8966 uout(i,j-1)=uout(i,j-1)*csfac
8967 vout(i,j)=vout(i,j)*csfac
8968 vout(i-1,j)=vout(i-1,j)*csfac
8970 if ( (is==1) .and. (je==npy-1) )
then 8973 uout(i,j)=uout(i,j)*csfac
8974 uout(i,j+1)=uout(i,j+1)*csfac
8975 vout(i,j+1)=vout(i,j+1)*csfac
8976 vout(i-1,j+1)=vout(i-1,j+1)*csfac
8978 if ( (ie==npx-1) .and. (je==npy-1) )
then 8981 uout(i+1,j)=uout(i+1,j)*csfac
8982 uout(i+1,j+1)=uout(i+1,j+1)*csfac
8983 vout(i,j+1)=vout(i,j+1)*csfac
8984 vout(i+1,j+1)=vout(i+1,j+1)*csfac
8986 if ( (ie==npx-1) .and. (js==1) )
then 8989 uout(i+1,j)=uout(i+1,j)*csfac
8990 uout(i+1,j-1)=uout(i+1,j-1)*csfac
8991 vout(i,j)=vout(i,j)*csfac
8992 vout(i+1,j)=vout(i+1,j)*csfac
8998 if (
present(nocomm))
then 8999 if (.not. nocomm)
call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
9001 call mpp_update_domains( uout,vout, domain, gridtype=cgrid_ne_param, complete=.true.)
9003 if (.not. bounded_domain)
call fill_corners(uout, vout, npx, npy, vector=.true., cgrid=.true.)
9015 subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng, bd)
9019 integer,
intent(IN) :: npx, npy, ng
9020 real ,
intent(IN) :: uin(bd%isd:bd%ied+1,bd%jsd:bd%jed )
9021 real ,
intent(IN) :: vin(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
9022 real ,
intent(OUT) :: uout(bd%isd:bd%ied ,bd%jsd:bd%jed )
9023 real ,
intent(OUT) :: vout(bd%isd:bd%ied ,bd%jsd:bd%jed )
9024 real ,
intent(IN),
dimension(bd%isd:bd%ied+1,bd%jsd:bd%jed) :: dxc, dy
9025 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed+1) :: dyc, dx
9026 real ,
intent(IN),
dimension(bd%isd:bd%ied,bd%jsd:bd%jed) :: dxa, dya
9030 real :: tmp1i(bd%isd:bd%ied+1)
9031 real :: tmp2i(bd%isd:bd%ied+1)
9032 real :: tmp3i(bd%isd:bd%ied+1)
9033 real :: tmp1j(bd%jsd:bd%jed+1)
9034 real :: tmp2j(bd%jsd:bd%jed+1)
9035 real :: tmp3j(bd%jsd:bd%jed+1)
9037 integer :: is, ie, js, je
9038 integer :: isd, ied, jsd, jed
9061 tmp2j(:) = vin(i,:)*dx(i,:)
9064 vout(i,jsd:jed) = tmp1j(jsd+1:jed+1)/dxa(i,jsd:jed)
9068 tmp2i(:) = uin(:,j)*dy(:,j)
9071 uout(isd:ied,j) = tmp1i(isd+1:ied+1)/dya(isd:ied,j)
9084 subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
9087 integer,
intent(IN) :: ndims
9088 real ,
intent(INOUT) :: myU
9089 real ,
intent(INOUT) :: myV
9090 real(kind=R_GRID) ,
intent(IN) :: p1(ndims)
9091 real(kind=R_GRID) ,
intent(IN) :: p2(ndims)
9092 real(kind=R_GRID) ,
intent(IN) :: p3(ndims)
9093 real(kind=R_GRID) ,
intent(IN) :: p4(ndims)
9094 real(kind=R_GRID) ,
intent(IN) :: t1(ndims)
9095 integer,
intent(IN) :: dir
9097 real(kind=R_GRID) :: ee1(3), ee2(3), ee3(3), elon(3), elat(3)
9099 real :: g11, g12, g21, g22
9105 elon(1) = -sin(t1(1) - pi)
9106 elon(2) = cos(t1(1) - pi)
9108 elat(1) = -sin(t1(2))*cos(t1(1) - pi)
9109 elat(2) = -sin(t1(2))*sin(t1(1) - pi)
9110 elat(3) = cos(t1(2))
9118 newu = myu*g11 + myv*g12
9119 newv = myu*g21 + myv*g22
9121 newu = ( myu*g22 - myv*g12)/(g11*g22 - g21*g12)
9122 newv = (-myu*g21 + myv*g11)/(g11*g22 - g21*g12)
9130 use mpp_parameter_mod
, only: dgrid_ne
9132 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1)
9133 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed )
9134 integer,
intent(IN) :: npx, npy
9135 type(domain2d),
intent(INOUT) :: domain
9137 call mpp_update_domains( u, v, domain, gridtype=dgrid_ne, complete=.true.)
9149 use mpp_parameter_mod
, only: dgrid_ne
9151 real ,
intent(INOUT) :: u(bd%isd:bd%ied ,bd%jsd:bd%jed+1,npz)
9152 real ,
intent(INOUT) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed ,npz)
9153 integer,
intent(IN) :: npx, npy, npz
9154 type(domain2d),
intent(INOUT) :: domain
9157 call mpp_update_domains( u, v, domain, gridtype=dgrid_ne, complete=.true.)
9169 real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
result (gsum)
9171 integer,
intent(IN) :: npx, npy
9172 integer,
intent(IN) :: ifirst, ilast
9173 integer,
intent(IN) :: jfirst, jlast
9174 integer,
intent(IN) :: isd, ied
9175 integer,
intent(IN) :: jsd, jed, tile
9176 real ,
intent(IN) :: p(ifirst:ilast,jfirst:jlast)
9182 real,
allocatable :: p_R8(:,:,:)
9184 real,
pointer,
dimension(:,:,:) :: agrid, grid
9185 real,
pointer,
dimension(:,:) :: area, rarea, fC, f0
9186 real,
pointer,
dimension(:,:) :: dx,dy, dxa,dya, rdxa, rdya, dxc,dyc
9188 logical,
pointer :: cubed_sphere, latlon
9190 logical,
pointer :: have_south_pole, have_north_pole
9192 integer,
pointer :: ntiles_g
9193 real,
pointer :: acapN, acapS, globalarea
9195 grid => gridstruct%grid
9196 agrid=> gridstruct%agrid
9198 area => gridstruct%area
9199 rarea => gridstruct%rarea
9206 dxa => gridstruct%dxa
9207 dya => gridstruct%dya
9208 rdxa => gridstruct%rdxa
9209 rdya => gridstruct%rdya
9210 dxc => gridstruct%dxc
9211 dyc => gridstruct%dyc
9213 cubed_sphere => gridstruct%cubed_sphere
9214 latlon => gridstruct%latlon
9216 have_south_pole => gridstruct%have_south_pole
9217 have_north_pole => gridstruct%have_north_pole
9219 ntiles_g => gridstruct%ntiles_g
9220 acapn => gridstruct%acapN
9221 acaps => gridstruct%acapS
9222 globalarea => gridstruct%globalarea
9224 allocate(p_r8(npx-1,npy-1,ntiles_g))
9231 gsum = gsum + p(1,1)*acaps
9232 gsum = gsum + p(1,npy-1)*acapn
9235 gsum = gsum + p(i,j)*cos(agrid(i,j,2))
9243 p_r8(i,j,n) = p(i,j)*area(i,j)
9247 call mp_gather(p_r8, ifirst,ilast, jfirst,jlast, npx-1, npy-1, ntiles_g)
9248 if (is_master())
then 9252 gsum = gsum + p_r8(i,j,n)
9256 gsum = gsum/globalarea
9258 call mpp_broadcast(gsum, mpp_root_pe())
9268 real(kind=R_GRID),
intent(in):: p1(2), p2(2), p3(2)
9269 real(kind=R_GRID),
intent(out):: uvect(3)
9272 real(kind=R_GRID) :: xyz1(3), xyz2(3), xyz3(3)
9279 uvect(n) = xyz3(n)-xyz1(n)
9292 integer,
intent(in):: np
9293 real(kind=R_GRID),
intent(inout):: e(3,np)
9299 pdot = sqrt(e(1,n)**2+e(2,n)**2+e(3,n)**2)
9301 e(k,n) = e(k,n) / pdot
9311 subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, &
9312 kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
9315 integer,
intent(in):: im, jm, km, nq
9316 integer,
intent(in):: ifirst, ilast
9317 integer,
intent(in):: jfirst, jlast
9318 integer,
intent(in):: kfirst, klast
9319 integer,
intent(in):: ng_e
9320 integer,
intent(in):: ng_w
9321 integer,
intent(in):: ng_s
9322 integer,
intent(in):: ng_n
9323 real,
intent(inout):: q_ghst(ifirst-ng_w:ilast+ng_e,jfirst-ng_s:jlast+ng_n,kfirst:klast,nq)
9324 real,
optional,
intent(in):: q(ifirst:ilast,jfirst:jlast,kfirst:klast,nq)
9338 if (
present(q))
then 9339 q_ghst(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq) = &
9340 q(ifirst:ilast,jfirst:jlast,kfirst:klast,1:nq)
9346 do j=jfirst-ng_s,jlast+ng_n
9348 q_ghst(ifirst-i,j,k,n) = q_ghst(ilast-i+1,j,k,n)
9351 q_ghst(ilast+i,j,k,n) = q_ghst(ifirst+i-1,j,k,n)
9374 integer,
intent(in):: ifirst,ilast
9375 real,
intent(out) :: qout(ifirst:)
9376 real,
intent(in) :: qin(ifirst:)
9377 real,
intent(in) :: dx(ifirst:)
9378 integer,
intent(in):: order
9381 real :: dm(ifirst:ilast),qmax,qmin
9382 real :: r3, da1, da2, a6da, a6, al, ar
9383 real :: qLa, qLb1, qLb2
9392 qout(i) = 0.5 * (qin(i-1) + qin(i))
9394 elseif (order==2)
then 9397 qout(i) = (dx(i-1)*qin(i-1) + dx(i)*qin(i))/(dx(i-1)+dx(i))
9399 elseif (order==3)
then 9402 do i=ifirst+1,ilast-1
9403 dm(i) = 0.25*(qin(i+1) - qin(i-1))
9408 do i=ifirst+1,ilast-1
9409 qmax = max(qin(i-1),qin(i),qin(i+1)) - qin(i)
9410 qmin = qin(i) - min(qin(i-1),qin(i),qin(i+1))
9411 dm(i) = sign(min(abs(dm(i)),qmin,qmax),dm(i))
9414 do i=ifirst+1,ilast-1
9415 qout(i) = 0.5*(qin(i-1)+qin(i)) + r3*(dm(i-1) - dm(i))
9422 qout(ifirst+1) = 0.5 * (qin(ifirst) + qin(ifirst+1))
9423 qout(ilast) = 0.5 * (qin(ilast-1) + qin(ilast))
9425 elseif (order==4)
then 9428 do i=ifirst+1,ilast-1
9429 dm(i) = ( (2.*dx(i-1) + dx(i) ) / &
9430 ( dx(i+1) + dx(i) ) ) * ( qin(i+1) - qin(i) ) + &
9431 ( (dx(i) + 2.*dx(i+1)) / &
9432 (dx(i-1) + dx(i) ) ) * ( qin(i) - qin(i-1) )
9433 dm(i) = ( dx(i) / ( dx(i-1) + dx(i) + dx(i+1) ) ) * dm(i)
9434 if ( (qin(i+1)-qin(i))*(qin(i)-qin(i-1)) > 0.)
then 9435 dm(i) = sign( min( abs(dm(i)), 2.*abs(qin(i)-qin(i-1)), 2.*abs(qin(i+1)-qin(i)) ) , dm(i) )
9441 do i=ifirst+2,ilast-1
9442 qla = ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) - &
9443 ( (dx(i+1) + dx(i)) / (2.*dx(i) + dx(i-1)) )
9444 qla = ( (2.*dx(i) * dx(i-1)) / (dx(i-1) + dx(i)) ) * qla * &
9446 qlb1 = dx(i-1) * ( (dx(i-2) + dx(i-1)) / (2.*dx(i-1) + dx(i)) ) * &
9448 qlb2 = dx(i) * ( (dx(i) + dx(i+1)) / (dx(i-1) + 2.*dx(i)) ) * &
9451 qout(i) = 1. / ( dx(i-2) + dx(i-1) + dx(i) + dx(i+1) )
9452 qout(i) = qout(i) * ( qla - qlb1 + qlb2 )
9453 qout(i) = qin(i-1) + ( dx(i-1) / ( dx(i-1) + dx(i) ) ) * (qin(i) - qin(i-1)) + qout(i)
9456 elseif (order==5)
then 9459 do i=ifirst+1,ilast-1
9460 x = float(i-(ifirst+1))*float(ilast-ifirst+1-1)/float(ilast-ifirst-1)
9461 qout(i) = qin(ifirst+nint(x)) + (x - nint(x)) * (qin(ifirst+nint(x+1)) - qin(ifirst+nint(x)))
9482 subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, &
9483 ng_d, ng_s, jfirst, jlast)
9490 integer,
intent(in):: ng_s, ng_d
9491 real,
intent(in):: coslon(im,jm), sinlon(im,jm)
9492 real,
intent(in):: cosl5(im,jm),sinl5(im,jm)
9493 real,
intent(in):: u(im,jfirst-ng_d:jlast+ng_s)
9496 real,
intent(inout):: v(im,jfirst-ng_d:jlast+ng_d)
9513 real uanp(im), uasp(im), vanp(im), vasp(im)
9514 real un, vn, us, vs, r2im
9517 r2im = 0.5d0/dble(im)
9522 if ( jfirst-ng_d <= 1 )
then 9524 uasp(i) = u(i, 2) + u(i,3)
9528 vasp(i) = v(i, 2) + v(i+1,2)
9530 vasp(im) = v(im,2) + v(1,2)
9536 us = us + (uasp(i+imh)-uasp(i))*sinlon(i,1) &
9537 + (vasp(i)-vasp(i+imh))*coslon(i,1)
9538 vs = vs + (uasp(i+imh)-uasp(i))*coslon(i,1) &
9539 + (vasp(i+imh)-vasp(i))*sinlon(i,1)
9547 v(i, 1) = us*cosl5(i,1) - vs*sinl5(i,1)
9548 v(i+imh,1) = -v(i,1)
9553 if ( jlast+ng_d >= jm )
then 9556 uanp(i) = u(i,jm-1) + u(i,jm)
9560 vanp(i) = v(i,jm-1) + v(i+1,jm-1)
9562 vanp(im) = v(im,jm-1) + v(1,jm-1)
9569 un = un + (uanp(i+imh)-uanp(i))*sinlon(i,jm) &
9570 + (vanp(i+imh)-vanp(i))*coslon(i,jm)
9571 vn = vn + (uanp(i)-uanp(i+imh))*coslon(i,jm) &
9572 + (vanp(i+imh)-vanp(i))*sinlon(i,jm)
9580 v(i, jm) = -un*cosl5(i,jm) - vn*sinl5(i,jm)
9581 v(i+imh,jm) = -v(i,jm)
9586 end subroutine vpol5 9588 subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
9590 character(len=*),
intent(in):: qname
9591 integer,
intent(in):: is, ie, js, je
9592 integer,
intent(in):: n_g, km
9593 real,
intent(in):: q(is-n_g:ie+n_g, js-n_g:je+n_g, km)
9594 real,
intent(in):: fac
9605 if( q(i,j,k) < qmin )
then 9607 elseif( q(i,j,k) > qmax )
then 9614 write(*,*) qname,
' max = ', qmax*fac,
' min = ', qmin*fac
9618 subroutine var_dz(km, ztop, ze)
9619 integer,
intent(in):: km
9620 real,
intent(in):: ztop
9621 real,
intent(out),
dimension(km+1):: ze
9623 real,
dimension(km):: dz, s_fac
9634 s_fac(k) = 1.05 * s_fac(k+1)
9636 s_fac(4) = 1.1*s_fac(5)
9637 s_fac(3) = 1.2*s_fac(4)
9638 s_fac(2) = 1.3*s_fac(3)
9639 s_fac(1) = 1.5*s_fac(2)
9643 sum1 = sum1 + s_fac(k)
9649 dz(k) = s_fac(k) * dz0
9654 ze(k) = ze(k+1) + dz(k)
9659 dz(k) = dz(k) * (ztop/ze(1))
9663 ze(k) = ze(k+1) + dz(k)
9667 call sm1_edge(1, 1, 1, 1, km, 1, 1, ze, 1)
9669 if ( is_master() )
then 9670 write(*,*)
'var_dz: model top (km)=', ztop*0.001
9672 dz(k) = ze(k) - ze(k+1)
9673 write(*,*) k, 0.5*(ze(k)+ze(k+1)),
'dz=', dz(k)
9679 subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
9680 integer,
intent(in):: is, ie, js, je, km
9681 integer,
intent(in):: ntimes, i, j
9682 real,
intent(inout):: ze(is:ie,js:je,km+1)
9684 real,
parameter:: df = 0.25
9687 integer k, n, k1, k2
9691 dz(k) = ze(i,j,k+1) - ze(i,j,k)
9700 flux(k) = df*(dz(k) - dz(k-1))
9704 dz(k) = dz(k) - flux(k) + flux(k+1)
9709 ze(i,j,k) = ze(i,j,k+1) - dz(k)
real, dimension(:,:), allocatable case9_b
subroutine, public cart_to_latlon(np, q, xs, ys)
subroutine, public mid_pt_sphere(p1, p2, pm)
subroutine dcmip16_tc(delp, pt, u, v, q, w, delz, is, ie, js, je, isd, ied, jsd, jed, npz, nq, ak, bk, ptop, pk, peln, pe, pkz, gz, phis, ps, grid, agrid, hydrostatic, nwat, adiabatic)
subroutine mp_ghost_ew(im, jm, km, nq, ifirst, ilast, jfirst, jlast, kfirst, klast, ng_w, ng_e, ng_s, ng_n, q_ghst, q)
subroutine interp_left_edge_1d(qout, qin, dx, ifirst, ilast, order)
subroutine dtoa(uin, vin, uout, vout, dx, dy, dxa, dya, dxc, dyc, npx, npy, ng, bd)
subroutine, public init_double_periodic(u, v, w, pt, delp, q, phis, ps, pe, peln, pk, pkz, uc, vc, ua, va, ak, bk, gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, ks, ptop, domain_in, tile_in, bd)
subroutine ctoa(uin, vin, uout, vout, dx, dy, dxc, dyc, dxa, dya, npx, npy, ng, bd)
subroutine, public init_case(u, v, w, pt, delp, q, phis, ps, pe, peln, pk, pkz, uc, vc, ua, va, ak, bk, gridstruct, flagstruct, npx, npy, npz, ng, ncnst, nwat, ndims, nregions, dry_mass, mountain, moist_phys, hydrostatic, hybrid_z, delz, ze0, adiabatic, ks, npx_global, ptop, domain_in, tile_in, bd)
The module 'fv_mp_mod' is a single program multiple data (SPMD) parallel decompostion/communication m...
integer, parameter interporder
subroutine get_vorticity(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, u, v, vort, dx, dy, rarea)
The type 'fv_grid_type' is made up of grid-dependent information from fv_grid_tools and fv_grid_utils...
real(kind=r_grid), parameter radius
subroutine, public hybrid_z_dz(km, dz, ztop, s_rate)
real tvort_orig
enstrophy (integral of total vorticity)
subroutine sm1_edge(is, ie, js, je, km, i, j, ze, ntimes)
subroutine init_winds(UBar, u, v, ua, va, uc, vc, defOnGrid, npx, npy, ng, ndims, nregions, bounded_domain, gridstruct, domain, tile, bd)
real, dimension(:), allocatable lats_table
real, parameter pi_shift
3.0*pi/4.
subroutine, public eqv_pot(theta_e, pt, delp, delz, peln, pkz, q, is, ie, js, je, ng, npz, hydrostatic, moist)
The subroutine 'eqv_pot' calculates the equivalent potential temperature.
subroutine, public set_hybrid_z(is, ie, js, je, ng, km, ztop, dz, rgrav, hs, ze, dz3)
subroutine, public case51_forcing(delp, uc, vc, u, v, ua, va, pe, time, dt, gridstruct, npx, npy, npz, ptop, domain, bd)
subroutine superk_u(km, zz, um, dudz)
real function dcmip16_bc_uwind(z, T, lat)
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, adiabatic, make_nh)
the subroutine 'p_var' computes auxiliary pressure variables for a hydrostatic state.
subroutine normalize_vect(np, e)
subroutine dcmip16_tc_uwind_pert(z, r, lon, lat, uu, vv)
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, bd)
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
The function 'g_sum' is the fast version of 'globalsum'.
subroutine var_dz(km, ztop, ze)
real function, public inner_prod(v1, v2)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine, public get_unit_vect2(e1, e2, uc)
integer, parameter, public r_grid
subroutine, public gw_1d(km, p0, ak, bk, ptop, ztop, pt1)
The module 'fv_sg' performs FV sub-grid mixing.
integer, parameter initwindscase2
subroutine terminator_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, km, q, delp, ncnst, lon, lat, bd)
subroutine supercell_sounding(km, ps, pk1, tp, qp)
real function gh_jet(npy, lat_in)
subroutine, public surfdrv(npx, npy, grid, agrid, area, dx, dy, dxa, dya, dxc, dyc, sin_sg, phis, stretch_fac, nested, bounded_domain, npx_global, domain, grid_number, bd)
subroutine, public make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
subroutine, public hydro_eq(km, is, ie, js, je, ps, hs, drym, delp, ak, bk, pt, delz, area, ng, mountain, hydrostatic, hybrid_z, domain)
The subroutine 'hydro_eq' computes a hydrostatically balanced and isothermal basic state from input h...
integer, parameter initwindscase6
subroutine rankine_vortex(ubar, r0, p1, u, v, grid, bd)
subroutine atod(uin, vin, uout, vout, dxa, dya, dxc, dyc, npx, npy, ng, bounded_domain, domain, bd)
real, parameter, public ptop_min
integer, parameter, public f_p
subroutine, public ppme(p, qe, delp, im, km)
subroutine, public compute_dz_l32(km, ztop, dz)
subroutine get_pt_on_great_circle(p1, p2, dist, heading, p3)
real function globalsum(p, npx, npy, ifirst, ilast, jfirst, jlast, isd, ied, jsd, jed, gridstruct, tile)
real function dcmip16_tc_temperature(z, r)
subroutine, public qsmith(im, km, k1, t, p, q, qs, dqdt)
subroutine, public case9_forcing1(phis, time_since_start, isd, ied, jsd, jed)
integer, parameter initwindscase9
real tmass_orig
total mass
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, bounded_domain, c2l_ord, bd)
The module 'fv_arrays' contains the 'fv_atmos_type' and associated datatypes.
real function, public great_circle_dist(q1, q2, radius)
integer, parameter initwindscase0
The module 'fv_eta' contains routine to set up the reference (Eulerian) pressure coordinate.
real function dcmip16_bc_sphum(p, ps, lat, lon)
subroutine, public read_namelist_test_case_nml(nml_filename)
real, dimension(:), allocatable, public pz0
subroutine get_case9_b(B, agrid, isd, ied, jsd, jed)
subroutine, public case9_forcing2(phis, isd, ied, jsd, jed)
subroutine mp_update_dwinds_2d(u, v, npx, npy, domain, bd)
real function dcmip16_tc_sphum(z)
subroutine, public project_sphere_v(np, f, e)
real, dimension(:,:,:), allocatable ua0
Validating U-Wind.
@ The module 'fv_diagnostics' contains routines to compute diagnosic fields.
subroutine balanced_k(km, is, ie, js, je, ng, ps0, ze1, ts1, qs1, uz1, dudz, pe, pk, pt, delz, zvir, ptop, ak, bk, agrid)
The module 'fv_grid_utils' contains routines for setting up and computing grid-related quantities...
subroutine get_unit_vector(p1, p2, p3, uvect)
subroutine mp_update_dwinds_3d(u, v, npx, npy, npz, domain, bd)
subroutine, public prt_maxmin(qname, q, is, ie, js, je, n_g, km, fac)
subroutine rotate_winds(myU, myV, p1, p2, p3, p4, t1, ndims, dir)
real, dimension(:), allocatable, public zz0
subroutine, public checker_tracers(i0, i1, j0, j1, ifirst, ilast, jfirst, jlast, nq, km, q, lon, lat, nx, ny, rn)
real function dcmip16_bc_uwind_pert(z, lat, lon)
subroutine prt_m1(qname, q, is, ie, js, je, n_g, km, fac)
real, dimension(:,:,:), allocatable va0
Validating V-Windfms_io_exit, get_tile_string, &.
real function dcmip16_bc_temperature(z, lat)
real, dimension(:,:,:), allocatable phi0
Validating Field.
real function dcmip16_tc_pressure(z, r)
real(kind=r_grid), parameter one
real function dcmip16_bc_pressure(z, lat)
subroutine superk_sounding(km, pe, p00, ze, pt, qz)
subroutine vpol5(u, v, im, jm, coslon, sinlon, cosl5, sinl5, ng_d, ng_s, jfirst, jlast)
integer, parameter initwindscase1
real, dimension(:), allocatable gh_table
subroutine atoc(uin, vin, uout, vout, dx, dy, dxa, dya, npx, npy, ng, bounded_domain, domain, bd, noComm)
integer, parameter initwindscase5
subroutine, public compute_dz_l101(km, ztop, dz)
subroutine, public latlon2xyz(p, e, id)
The subroutine 'latlon2xyz' maps (lon, lat) to (x,y,z)