22 #include <fms_platform.h> 35 use fv_mp_nlm_mod,
only: mp_reduce_sum, mp_reduce_min, mp_reduce_max
42 #ifdef NO_QUAD_PRECISION 44 integer,
parameter::
f_p = selected_real_kind(15)
47 integer,
parameter::
f_p = selected_real_kind(20)
52 real(kind=R_GRID) ::
radius=cnst_radius
71 MODULE PROCEDURE fill_ghost_r4
78 subroutine grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order)
81 logical,
intent(in):: non_ortho
82 integer,
intent(in):: npx, npy, npz
83 integer,
intent(in)::
grid_type, c2l_order
93 real(kind=R_GRID) grid3(3,atm%bd%isd:atm%bd%ied+1,atm%bd%jsd:atm%bd%jed+1)
94 real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3), pp(3), ex(3), ey(3), e1(3), e2(3)
95 real(kind=R_GRID) pp1(2), pp2(2), pp3(2)
96 real(kind=R_GRID) sin2, tmp1, tmp2
97 integer i, j, k, n,
ip 99 integer :: is, ie, js, je
100 integer :: isd, ied, jsd, jed
103 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: agrid, grid
104 real(kind=R_GRID),
pointer,
dimension(:,:) :: area, area_c
105 real(kind=R_GRID),
pointer,
dimension(:,:) :: sina, cosa, dx, dy, dxc, dyc, dxa, dya
106 real,
pointer,
dimension(:,:) :: del6_u, del6_v
107 real,
pointer,
dimension(:,:) :: divg_u, divg_v
108 real,
pointer,
dimension(:,:) :: cosa_u, cosa_v, cosa_s
109 real,
pointer,
dimension(:,:) :: sina_u, sina_v
110 real,
pointer,
dimension(:,:) :: rsin_u, rsin_v
111 real,
pointer,
dimension(:,:) :: rsina, rsin2
112 real,
pointer,
dimension(:,:,:) :: sin_sg, cos_sg
113 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, ec1, ec2
114 real(kind=R_GRID),
pointer,
dimension(:,:,:,:) :: ew, es
115 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: en1, en2
117 logical,
pointer :: sw_corner, se_corner, ne_corner, nw_corner
129 agrid => atm%gridstruct%agrid_64
130 grid => atm%gridstruct%grid_64
131 area => atm%gridstruct%area_64
132 area_c => atm%gridstruct%area_c_64
133 dx => atm%gridstruct%dx_64
134 dy => atm%gridstruct%dy_64
135 dxc => atm%gridstruct%dxc_64
136 dyc => atm%gridstruct%dyc_64
137 dxa => atm%gridstruct%dxa_64
138 dya => atm%gridstruct%dya_64
139 sina => atm%gridstruct%sina_64
140 cosa => atm%gridstruct%cosa_64
142 divg_u => atm%gridstruct%divg_u
143 divg_v => atm%gridstruct%divg_v
145 del6_u => atm%gridstruct%del6_u
146 del6_v => atm%gridstruct%del6_v
148 cosa_u => atm%gridstruct%cosa_u
149 cosa_v => atm%gridstruct%cosa_v
150 cosa_s => atm%gridstruct%cosa_s
151 sina_u => atm%gridstruct%sina_u
152 sina_v => atm%gridstruct%sina_v
153 rsin_u => atm%gridstruct%rsin_u
154 rsin_v => atm%gridstruct%rsin_v
155 rsina => atm%gridstruct%rsina
156 rsin2 => atm%gridstruct%rsin2
157 ee1 => atm%gridstruct%ee1
158 ee2 => atm%gridstruct%ee2
159 ec1 => atm%gridstruct%ec1
160 ec2 => atm%gridstruct%ec2
161 ew => atm%gridstruct%ew
162 es => atm%gridstruct%es
163 sin_sg => atm%gridstruct%sin_sg
164 cos_sg => atm%gridstruct%cos_sg
165 en1 => atm%gridstruct%en1
166 en2 => atm%gridstruct%en2
170 sw_corner => atm%gridstruct%sw_corner
171 se_corner => atm%gridstruct%se_corner
172 ne_corner => atm%gridstruct%ne_corner
173 nw_corner => atm%gridstruct%nw_corner
175 if ( atm%flagstruct%do_schmidt .and. abs(atm%flagstruct%stretch_fac-1.) > 1.e-5 )
then 176 atm%gridstruct%stretched_grid = .true.
179 atm%gridstruct%stretched_grid = .false.
190 elseif ( .not. atm%flagstruct%hybrid_z )
then 192 call set_eta(npz, atm%ks, atm%ptop, atm%ak, atm%bk)
193 if ( is_master() )
then 194 write(*,*)
'Grid_init', npz, atm%ks, atm%ptop
195 tmp1 = atm%ak(atm%ks+1)
197 tmp1 =
max(tmp1, (atm%ak(k)-atm%ak(k+1))/
max(1.e-9, (atm%bk(k+1)-atm%bk(k))) )
199 write(*,*)
'Hybrid Sigma-P: minimum allowable surface pressure (hpa)=', tmp1/100.
200 if ( tmp1 > 420.e2 )
write(*,*)
'Warning: the chosen setting in set_eta can cause instability' 206 if (.not.
allocated(sst_ncep))
allocate (sst_ncep(i_sst,j_sst))
207 if (.not.
allocated(sst_anom))
allocate (sst_anom(i_sst,j_sst))
219 if (
grid_type < 3 .and. .not. atm%neststruct%nested)
then 220 if ( is==1 .and. js==1 ) sw_corner = .true.
221 if ( (ie+1)==npx .and. js==1 ) se_corner = .true.
222 if ( (ie+1)==npx .and. (je+1)==npy ) ne_corner = .true.
223 if ( is==1 .and. (je+1)==npy ) nw_corner = .true.
226 if ( sw_corner )
then 229 write(*,*)
'Corner interpolation coefficient=', tmp2/(tmp2-tmp1)
233 if ( .not. atm%neststruct%nested )
then 234 call fill_corners(grid(:,:,1), npx, npy, fill=xdir, bgrid=.true.)
235 call fill_corners(grid(:,:,2), npx, npy, fill=xdir, bgrid=.true.)
248 if (.not. atm%neststruct%nested)
then 258 if ( ( (i<1 .and. j<1 ) .or. (i>npx .and. j<1 ) .or. &
259 (i>npx .and. j>(npy-1)) .or. (i<1 .and. j>(npy-1)) ) .and. .not. atm%neststruct%nested)
then 262 call mid_pt_cart( grid(i,j,1:2), grid(i,j+1,1:2), pp)
263 if (i==1 .and. .not. atm%neststruct%nested)
then 266 elseif(i==npx .and. .not. atm%neststruct%nested)
then 277 call vect_cross(p1, grid3(1,i,j), grid3(1,i,j+1))
286 if ( ( (i<1 .and. j<1 ) .or. (i>(npx-1) .and. j<1 ) .or. &
287 (i>(npx-1) .and. j>npy) .or. (i<1 .and. j>npy) ) .and. .not. atm%neststruct%nested)
then 290 call mid_pt_cart(grid(i,j,1:2), grid(i+1,j,1:2), pp)
291 if (j==1 .and. .not. atm%neststruct%nested)
then 294 elseif (j==npy .and. .not. atm%neststruct%nested)
then 305 call vect_cross(p3, grid3(1,i,j), grid3(1,i+1,j))
322 cos_sg(i,j,6) =
cos_angle( grid3(1,i,j), grid3(1,i+1,j), grid3(1,i,j+1) )
324 cos_sg(i,j,7) = -
cos_angle( grid3(1,i+1,j), grid3(1,i,j), grid3(1,i+1,j+1) )
326 cos_sg(i,j,8) =
cos_angle( grid3(1,i+1,j+1), grid3(1,i+1,j), grid3(1,i,j+1) )
328 cos_sg(i,j,9) = -
cos_angle( grid3(1,i,j+1), grid3(1,i,j), grid3(1,i+1,j+1) )
338 cos_sg(i,j,1) =
cos_angle( p1, p3, grid3(1,i,j+1) )
340 cos_sg(i,j,2) =
cos_angle( p1, grid3(1,i+1,j), p3 )
342 cos_sg(i,j,3) =
cos_angle( p1, p3, grid3(1,i+1,j) )
344 cos_sg(i,j,4) =
cos_angle( p1, grid3(1,i,j+1), p3 )
347 cos_sg(i,j,5) =
inner_prod( ec1(1:3,i,j), ec2(1:3,i,j) )
354 sin_sg(i,j,
ip) =
min(1.0, sqrt(
max(0., 1.-cos_sg(i,j,
ip)**2) ) )
362 if (.not. atm%neststruct%nested)
then 363 if ( sw_corner )
then 365 sin_sg(0,i,3) = sin_sg(i,1,2)
366 sin_sg(i,0,4) = sin_sg(1,i,1)
369 if ( nw_corner )
then 371 sin_sg(0,i,3) = sin_sg(npy-i,npy-1,4)
374 sin_sg(i,npy,2) = sin_sg(1,npx+i,1)
377 if ( se_corner )
then 379 sin_sg(npx,j,1) = sin_sg(npx-j,1,2)
382 sin_sg(i,0,4) = sin_sg(npx-1,npx-i,3)
385 if ( ne_corner )
then 387 sin_sg(npx,i,1) = sin_sg(i,npy-1,4)
388 sin_sg(i,npy,2) = sin_sg(npx-1,i,3)
396 pp1(:) = grid(i ,j ,1:2)
397 pp2(:) = grid(i,j+1 ,1:2)
401 atm%gridstruct%l2c_v(i,j) = cos(pp3(2)) *
inner_prod(e2, ex)
406 pp1(:) = grid(i, j,1:2)
407 pp2(:) = grid(i+1,j,1:2)
411 atm%gridstruct%l2c_u(i,j) = cos(pp3(2)) *
inner_prod(e1, ex)
444 if ( non_ortho )
then 460 if (i==1 .and. .not. atm%neststruct%nested)
then 461 call vect_cross(pp, grid3(1,i, j), grid3(1,i+1,j))
462 elseif(i==npx .and. .not. atm%neststruct%nested)
then 463 call vect_cross(pp, grid3(1,i-1,j), grid3(1,i, j))
465 call vect_cross(pp, grid3(1,i-1,j), grid3(1,i+1,j))
467 call vect_cross(ee1(1:3,i,j), pp, grid3(1:3,i,j))
471 if (j==1 .and. .not. atm%neststruct%nested)
then 472 call vect_cross(pp, grid3(1:3,i,j ), grid3(1:3,i,j+1))
473 elseif(j==npy .and. .not. atm%neststruct%nested)
then 474 call vect_cross(pp, grid3(1:3,i,j-1), grid3(1:3,i,j ))
476 call vect_cross(pp, grid3(1:3,i,j-1), grid3(1:3,i,j+1))
478 call vect_cross(ee2(1:3,i,j), pp, grid3(1:3,i,j))
484 cosa(i,j) = sign(
min(1., abs(tmp1)), tmp1)
485 sina(i,j) = sqrt(
max(0.,1. -cosa(i,j)**2))
487 cosa(i,j) = 0.5*(cos_sg(i-1,j-1,8)+cos_sg(i,j,6))
488 sina(i,j) = 0.5*(sin_sg(i-1,j-1,8)+sin_sg(i,j,6))
500 cosa_u(i,j) = 0.5*(cos_sg(i-1,j,3)+cos_sg(i,j,1))
501 sina_u(i,j) = 0.5*(sin_sg(i-1,j,3)+sin_sg(i,j,1))
508 cosa_v(i,j) = 0.5*(cos_sg(i,j-1,4)+cos_sg(i,j,2))
509 sina_v(i,j) = 0.5*(sin_sg(i,j-1,4)+sin_sg(i,j,2))
517 cosa_s(i,j) = cos_sg(i,j,5)
523 if (.not. atm%neststruct%nested)
then 531 if ( i==npx .and. j==npy .and. .not. atm%neststruct%nested)
then 532 else if ( ( i==1 .or. i==npx .or. j==1 .or. j==npy ) .and. .not. atm%neststruct%nested )
then 543 if ( (i==1 .or. i==npx) .and. .not. atm%neststruct%nested )
then 545 rsin_u(i,j) = 1. / sign(
max(
tiny_number,abs(sina_u(i,j))), sina_u(i,j))
552 if ( (j==1 .or. j==npy) .and. .not. atm%neststruct%nested )
then 554 rsin_v(i,j) = 1. / sign(
max(
tiny_number,abs(sina_v(i,j))), sina_v(i,j))
561 if (.not. atm%neststruct%nested)
then 571 if ( sw_corner )
then 573 sin_sg(0,i,3) = sin_sg(i,1,2)
574 sin_sg(i,0,4) = sin_sg(1,i,1)
575 cos_sg(0,i,3) = cos_sg(i,1,2)
576 cos_sg(i,0,4) = cos_sg(1,i,1)
585 if ( nw_corner )
then 587 sin_sg(0,i,3) = sin_sg(npy-i,npy-1,4)
588 cos_sg(0,i,3) = cos_sg(npy-i,npy-1,4)
593 sin_sg(i,npy,2) = sin_sg(1,npy-i,1)
594 cos_sg(i,npy,2) = cos_sg(1,npy-i,1)
600 if ( se_corner )
then 602 sin_sg(npx,j,1) = sin_sg(npx-j,1,2)
603 cos_sg(npx,j,1) = cos_sg(npx-j,1,2)
608 sin_sg(i,0,4) = sin_sg(npx-1,npx-i,3)
609 cos_sg(i,0,4) = cos_sg(npx-1,npx-i,3)
615 if ( ne_corner )
then 617 sin_sg(npx,npy+i,1) = sin_sg(npx+i,npy-1,4)
618 sin_sg(npx+i,npy,2) = sin_sg(npx-1,npy+i,3)
619 cos_sg(npx,npy+i,1) = cos_sg(npx+i,npy-1,4)
622 cos_sg(npx+i,npy,2) = cos_sg(npx-1,npy+i,3)
650 if (.not. atm%neststruct%nested)
then 654 call vect_cross(ew(1,i,j,1), grid3(1,i,j+1), grid3(1,i,j))
657 if ( (ie+1)==npx )
then 659 call vect_cross(ew(1,i,j,1), grid3(1,i,j+1), grid3(1,i,j))
667 call vect_cross(es(1,i,j,2), grid3(1,i,j),grid3(1,i+1,j))
671 if ( (je+1)==npy )
then 674 call vect_cross(es(1,i,j,2), grid3(1,i,j),grid3(1,i+1,j))
685 call vect_cross(en1(1:3,i,j), grid3(1,i,j), grid3(1,i+1,j))
691 call vect_cross(en2(1:3,i,j), grid3(1,i,j+1), grid3(1,i,j))
701 if ((j==1 .OR. j==npy) .and. .not. atm%neststruct%nested)
then 703 divg_u(i,j) = 0.5*(sin_sg(i,j,2)+sin_sg(i,j-1,4))*dyc(i,j)/dx(i,j)
704 del6_u(i,j) = 0.5*(sin_sg(i,j,2)+sin_sg(i,j-1,4))*dx(i,j)/dyc(i,j)
708 divg_u(i,j) = sina_v(i,j)*dyc(i,j)/dx(i,j)
709 del6_u(i,j) = sina_v(i,j)*dx(i,j)/dyc(i,j)
715 divg_v(i,j) = sina_u(i,j)*dxc(i,j)/dy(i,j)
716 del6_v(i,j) = sina_u(i,j)*dy(i,j)/dxc(i,j)
718 if (is == 1 .and. .not. atm%neststruct%nested)
then 719 divg_v(is,j) = 0.5*(sin_sg(1,j,1)+sin_sg(0,j,3))*dxc(is,j)/dy(is,j)
720 del6_v(is,j) = 0.5*(sin_sg(1,j,1)+sin_sg(0,j,3))*dy(is,j)/dxc(is,j)
722 if (ie+1 == npx .and. .not. atm%neststruct%nested)
then 723 divg_v(ie+1,j) = 0.5*(sin_sg(npx,j,1)+sin_sg(npx-1,j,3))*dxc(ie+1,j)/dy(ie+1,j)
724 del6_v(ie+1,j) = 0.5*(sin_sg(npx,j,1)+sin_sg(npx-1,j,3))*dy(ie+1,j)/dxc(ie+1,j)
731 call global_mx(area, ng, atm%gridstruct%da_min, atm%gridstruct%da_max, atm%bd)
732 if( is_master() )
write(*,*)
'da_max/da_min=', atm%gridstruct%da_max/atm%gridstruct%da_min
734 call global_mx_c(area_c(is:ie,js:je), is, ie, js, je, atm%gridstruct%da_min_c, atm%gridstruct%da_max_c)
736 if( is_master() )
write(*,*)
'da_max_c/da_min_c=', atm%gridstruct%da_max_c/atm%gridstruct%da_min_c
742 if (
grid_type < 3 .and. .not. atm%neststruct%nested)
then 744 gridtype=cgrid_ne_param, complete=.true.)
746 gridtype=cgrid_ne_param, complete=.true.)
747 call edge_factors (atm%gridstruct%edge_s, atm%gridstruct%edge_n, atm%gridstruct%edge_w, &
748 atm%gridstruct%edge_e, non_ortho, grid, agrid, npx, npy, atm%bd)
749 call efactor_a2c_v(atm%gridstruct%edge_vect_s, atm%gridstruct%edge_vect_n, &
750 atm%gridstruct%edge_vect_w, atm%gridstruct%edge_vect_e, &
751 non_ortho, grid, agrid, npx, npy, atm%neststruct%nested, atm%bd)
769 atm%gridstruct%grid = atm%gridstruct%grid_64
770 atm%gridstruct%agrid = atm%gridstruct%agrid_64
771 atm%gridstruct%area = atm%gridstruct%area_64
772 atm%gridstruct%area_c = atm%gridstruct%area_c_64
773 atm%gridstruct%dx = atm%gridstruct%dx_64
774 atm%gridstruct%dy = atm%gridstruct%dy_64
775 atm%gridstruct%dxa = atm%gridstruct%dxa_64
776 atm%gridstruct%dya = atm%gridstruct%dya_64
777 atm%gridstruct%dxc = atm%gridstruct%dxc_64
778 atm%gridstruct%dyc = atm%gridstruct%dyc_64
779 atm%gridstruct%cosa = atm%gridstruct%cosa_64
780 atm%gridstruct%sina = atm%gridstruct%sina_64
786 deallocate ( atm%gridstruct%area_c_64 )
789 deallocate ( atm%gridstruct%dxa_64 )
790 deallocate ( atm%gridstruct%dya_64 )
791 deallocate ( atm%gridstruct%dxc_64 )
792 deallocate ( atm%gridstruct%dyc_64 )
793 deallocate ( atm%gridstruct%cosa_64 )
794 deallocate ( atm%gridstruct%sina_64 )
845 if (
allocated(sst_ncep))
deallocate( sst_ncep )
846 if (
allocated(sst_anom))
deallocate( sst_anom )
857 real(kind=R_GRID),
intent(in):: c
858 real(kind=R_GRID),
intent(in):: lon_p, lat_p
859 integer,
intent(in):: n
860 integer,
intent(in):: i1, i2, j1, j2
862 real(kind=R_GRID),
intent(inout),
dimension(i1:i2,j1:j2):: lon, lat
864 real(f_p):: lat_t, sin_p, cos_p, sin_lat, cos_lat, sin_o, p2, two_pi
865 real(f_p):: c2p1, c2m1
871 if( is_master() .and. n==1 )
then 872 write(*,*) n,
'Schmidt transformation: stretching factor=', c,
' center=', lon_p, lat_p
883 if ( abs(c2m1) > 1.d-7 )
then 884 sin_lat = sin(lat(i,j))
885 lat_t = asin( (c2m1+c2p1*sin_lat)/(c2p1+c2m1*sin_lat) )
891 sin_o = -(sin_p*sin_lat + cos_p*cos_lat*cos(lon(i,j)))
892 if ( (1.-abs(sin_o)) < 1.d-7 )
then 894 lat(i,j) = sign( p2, sin_o )
896 lat(i,j) = asin( sin_o )
897 lon(i,j) = lon_p + atan2( -cos_lat*sin(lon(i,j)), &
898 -sin_lat*cos_p+cos_lat*sin_p*cos(lon(i,j)))
899 if ( lon(i,j) < 0.d0 )
then 900 lon(i,j) = lon(i,j) + two_pi
901 elseif( lon(i,j) >= two_pi )
then 902 lon(i,j) = lon(i,j) - two_pi
912 real(kind=R_GRID),
intent(in):: v1(3), v2(3)
913 real (f_p) :: vp1(3), vp2(3), prod16
917 vp1(k) =
real(v1(k),kind=
f_p)
918 vp2(k) =
real(v2(k),kind=
f_p)
920 prod16 = vp1(1)*vp2(1) + vp1(2)*vp2(2) + vp1(3)*vp2(3)
926 subroutine efactor_a2c_v(edge_vect_s, edge_vect_n, edge_vect_w, edge_vect_e, non_ortho, grid, agrid, npx, npy, nested, bd)
931 type(fv_grid_bounds_type),
intent(IN) :: bd
932 real(kind=R_GRID),
intent(INOUT),
dimension(bd%isd:bd%ied) :: edge_vect_s, edge_vect_n
933 real(kind=R_GRID),
intent(INOUT),
dimension(bd%jsd:bd%jed) :: edge_vect_w, edge_vect_e
934 logical,
intent(in):: non_ortho, nested
935 real(kind=R_GRID),
intent(in):: grid(bd%isd:bd%ied+1,bd%jsd:bd%jed+1,2)
936 real(kind=R_GRID),
intent(in):: agrid(bd%isd:bd%ied ,bd%jsd:bd%jed ,2)
937 integer,
intent(in):: npx, npy
939 real(kind=R_GRID) px(2,bd%isd:bd%ied+1), py(2,bd%jsd:bd%jed+1)
940 real(kind=R_GRID) p1(2,bd%isd:bd%ied+1), p2(2,bd%jsd:bd%jed+1)
941 real(kind=R_GRID) d1, d2
945 integer :: is, ie, js, je
946 integer :: isd, ied, jsd, jed
958 if ( .not. non_ortho )
then 969 if ( npx /= npy .and. .not. nested)
call mpp_error(fatal,
'efactor_a2c_v: npx /= npy')
970 if ( (npx/2)*2 == npx )
call mpp_error(fatal,
'efactor_a2c_v: npx/npy is not an odd number')
978 call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j, 1:2), py(1,j))
979 call mid_pt_sphere( grid(i, j,1:2), grid(i,j+1,1:2), p2(1,j))
990 edge_vect_w(j) = d1 / ( d1 + d2 )
994 edge_vect_w(j) = d1 / ( d2 + d1 )
998 edge_vect_w(0) = edge_vect_w(1)
1000 if ( (je+1)==npy )
then 1001 edge_vect_w(npy) = edge_vect_w(je)
1008 if ( (ie+1)==npx )
then 1011 call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j, 1:2), py(1,j))
1012 call mid_pt_sphere( grid(i, j,1:2), grid(i,j+1,1:2), p2(1,j))
1019 edge_vect_e(j) = d1 / ( d1 + d2 )
1023 edge_vect_e(j) = d1 / ( d2 + d1 )
1027 edge_vect_e(0) = edge_vect_e(1)
1029 if ( (je+1)==npy )
then 1030 edge_vect_e(npy) = edge_vect_e(je)
1040 call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i, j,1:2), px(1,i))
1041 call mid_pt_sphere( grid(i,j, 1:2), grid(i+1,j,1:2), p1(1,i))
1051 edge_vect_s(i) = d1 / ( d1 + d2 )
1055 edge_vect_s(i) = d1 / ( d2 + d1 )
1059 edge_vect_s(0) = edge_vect_s(1)
1061 if ( (ie+1)==npx )
then 1062 edge_vect_s(npx) = edge_vect_s(ie)
1070 if ( (je+1)==npy )
then 1074 call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i, j,1:2), px(1,i))
1075 call mid_pt_sphere( grid(i,j, 1:2), grid(i+1,j,1:2), p1(1,i))
1082 edge_vect_n(i) = d1 / ( d1 + d2 )
1086 edge_vect_n(i) = d1 / ( d2 + d1 )
1090 edge_vect_n(0) = edge_vect_n(1)
1092 if ( (ie+1)==npx )
then 1093 edge_vect_n(npx) = edge_vect_n(ie)
1105 subroutine edge_factors(edge_s, edge_n, edge_w, edge_e, non_ortho, grid, agrid, npx, npy, bd)
1110 type(fv_grid_bounds_type),
intent(IN) :: bd
1111 real(kind=R_GRID),
intent(INOUT),
dimension(npx) :: edge_s, edge_n
1112 real(kind=R_GRID),
intent(INOUT),
dimension(npy) :: edge_w, edge_e
1113 logical,
intent(in):: non_ortho
1114 real(kind=R_GRID),
intent(in):: grid(bd%isd:bd%ied+1,bd%jsd:bd%jed+1,2)
1115 real(kind=R_GRID),
intent(in):: agrid(bd%isd:bd%ied ,bd%jsd:bd%jed ,2)
1116 integer,
intent(in):: npx, npy
1118 real(kind=R_GRID) px(2,npx), py(2,npy)
1119 real(kind=R_GRID) d1, d2
1122 integer :: is, ie, js, je
1123 integer :: isd, ied, jsd, jed
1134 if ( .not. non_ortho )
then 1151 do j=
max(1,js-1),
min(npy-1,je+1)
1152 call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j,1:2), py(1,j))
1154 do j=
max(2,js),
min(npy-1,je+1)
1157 edge_w(j) = d2 / ( d1 + d2 )
1165 if ( (ie+1)==npx )
then 1167 do j=
max(1,js-1),
min(npy-1,je+1)
1168 call mid_pt_sphere(agrid(i-1,j,1:2), agrid(i,j,1:2), py(1,j))
1170 do j=
max(2,js),
min(npy-1,je+1)
1173 edge_e(j) = d2 / ( d1 + d2 )
1186 do i=
max(1,is-1),
min(npx-1,ie+1)
1187 call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i,j,1:2), px(1,i))
1189 do i=
max(2,is),
min(npx-1,ie+1)
1192 edge_s(i) = d2 / ( d1 + d2 )
1200 if ( (je+1)==npy )
then 1202 do i=
max(1,is-1),
min(npx-1,ie+1)
1203 call mid_pt_sphere(agrid(i,j-1,1:2), agrid(i,j,1:2), px(1,i))
1205 do i=
max(2,is),
min(npx-1,ie+1)
1208 edge_n(i) = d2 / ( d1 + d2 )
1219 real(kind=R_GRID),
intent(out):: lon(im+1,im+1)
1220 real(kind=R_GRID),
intent(out):: lat(im+1,im+1)
1232 lon(i,j) = lon(i,j) -
pi 1252 integer,
intent(in):: im
1253 real(kind=R_GRID),
intent(out):: lamda(im+1,im+1)
1254 real(kind=R_GRID),
intent(out):: theta(im+1,im+1)
1257 real(kind=R_GRID) pp(3,im+1,im+1)
1258 real(kind=R_GRID) p1(2), p2(2)
1259 real(f_p):: rsq3, alpha, delx, dely
1262 rsq3 = 1.d0/sqrt(3.d0)
1263 alpha = asin( rsq3 )
1269 dely = 2.d0*alpha /
real(im,kind=
f_p)
1273 lamda(1, j) = 0.75d0*
pi 1274 lamda(im+1,j) = 1.25d0*
pi 1275 theta(1, j) = -alpha + dely*
real(j-1,kind=f_p) 1276 theta(im+1,j) = theta(1,j)
1282 call mirror_latlon(lamda(1,1), theta(1,1), lamda(im+1,im+1), theta(im+1,im+1), &
1283 lamda(1,i), theta(1,i), lamda(i,1), theta(i, 1) )
1284 lamda(i,im+1) = lamda(i,1)
1285 theta(i,im+1) = -theta(i,1)
1289 call latlon2xyz2(lamda(1 , 1), theta(1, 1), pp(1, 1, 1))
1290 call latlon2xyz2(lamda(im+1, 1), theta(im+1, 1), pp(1,im+1, 1))
1291 call latlon2xyz2(lamda(1, im+1), theta(1, im+1), pp(1, 1,im+1))
1292 call latlon2xyz2(lamda(im+1,im+1), theta(im+1,im+1), pp(1,im+1,im+1))
1299 call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,j))
1300 pp(2,i,j) = -pp(2,i,j)*rsq3/pp(1,i,j)
1301 pp(3,i,j) = -pp(3,i,j)*rsq3/pp(1,i,j)
1306 call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,1))
1307 pp(2,i,1) = -pp(2,i,1)*rsq3/pp(1,i,1)
1308 pp(3,i,1) = -pp(3,i,1)*rsq3/pp(1,i,1)
1320 pp(2,i,j) = pp(2,i,1)
1322 pp(3,i,j) = pp(3,1,j)
1329 if ( is_master() )
then 1330 p1(1) = lamda(1,1); p1(2) = theta(1,1)
1331 p2(1) = lamda(2,1); p2(2) = theta(2,1)
1350 integer,
intent(IN) :: im, in, nghost
1351 real(kind=R_GRID),
intent(IN),
dimension(2) :: lL, lR, uL, uR
1352 real(kind=R_GRID),
intent(OUT) :: lamda(1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1353 real(kind=R_GRID),
intent(OUT) :: theta(1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1356 real(kind=R_GRID) pp(3,1-nghost:im+1+nghost,1-nghost:in+1+nghost)
1357 real(kind=R_GRID) p1(2), p2(2)
1358 real(f_p):: rsq3, alpha, delx, dely
1359 integer i, j, k, irefl
1361 rsq3 = 1.d0/sqrt(3.d0)
1362 alpha = asin( rsq3 )
1364 lamda(1,1) = ll(1); theta(1,1) = ll(2)
1365 lamda(im+1,1) = lr(1); theta(im+1,1) = lr(2)
1366 lamda(1,in+1) = ul(1); theta(1,in+1) = ul(2)
1367 lamda(im+1,in+1) = ur(1); theta(im+1,in+1) = ur(2)
1371 dely = (ul(2) - ll(2))/in
1373 theta(1,j) = theta(1,j-1) + dely
1374 theta(in+1,j) = theta(in+1,j-1) + dely
1375 lamda(1,j) = lamda(1,1)
1376 lamda(in+1,j) = lamda(in+1,1)
1379 theta(1,j) = theta(1,j+1) - dely
1380 theta(in+1,j) = theta(in+1,j+1) - dely
1381 lamda(1,j) = lamda(1,1)
1382 lamda(in+1,j) = lamda(in+1,1)
1385 lamda(1,:) = lamda(1,1)
1386 lamda(in+1,:) = lamda(in+1,1)
1390 do i=1-nghost,im+1+nghost
1395 (/lamda(1,1),theta(1,1)/), (/lamda(im+1,1),theta(im+1,1)/), p1 )
1397 (/lamda(1,in+1),theta(1,in+1)/), (/lamda(im+1,in+1),theta(im+1,in+1)/), p2 )
1399 lamda(i,1) = p1(1); theta(i,1) = p1(2)
1400 lamda(i,in+1) = p2(1); theta(i,in+1) = p2(2)
1407 do j=1-nghost,in+1+nghost
1408 call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,j))
1409 pp(2,i,j) = -pp(2,i,j)*rsq3/pp(1,i,j)
1410 pp(3,i,j) = -pp(3,i,j)*rsq3/pp(1,i,j)
1414 do i=1-nghost,im+1+nghost
1415 call latlon2xyz2(lamda(i,j), theta(i,j), pp(1,i,1))
1416 pp(2,i,1) = -pp(2,i,1)*rsq3/pp(1,i,1)
1417 pp(3,i,1) = -pp(3,i,1)*rsq3/pp(1,i,1)
1422 do j=1-nghost,in+1+nghost
1423 do i=1-nghost,im+1+nghost
1428 do j=1-nghost,in+1+nghost
1429 do i=1-nghost,im+1+nghost
1431 pp(2,i,j) = pp(2,i,1)
1433 pp(3,i,j) = pp(3,1,j)
1438 pp(:,1-nghost:im+1+nghost,1-nghost:in+1+nghost), &
1439 lamda(1-nghost:im+1+nghost,1-nghost:in+1+nghost), &
1440 theta(1-nghost:im+1+nghost,1-nghost:in+1+nghost))
1444 if ( is_master() )
then 1445 p1(1) = lamda(1,1); p1(2) = theta(1,1)
1446 p2(1) = lamda(2,1); p2(2) = theta(2,1)
1448 p2(1) = lamda(1,2); p2(2) = theta(1,2)
1461 real(kind=R_GRID) lamda(im+1,im+1)
1462 real(kind=R_GRID) theta(im+1,im+1)
1463 real(kind=R_GRID) p(3,im+1,im+1)
1465 real(kind=R_GRID) rsq3, xf, y0, z0, y, x, z, ds
1466 real(kind=R_GRID) dy, dz
1468 real(kind=R_GRID) dp
1472 rsq3 = 1.d0/sqrt(3.d0)
1476 p(2,j,k) =-rsq3*tan(-0.25d0*
pi+(j-1)*dp)
1477 p(3,j,k) = rsq3*tan(-0.25d0*
pi+(k-1)*dp)
1488 real(kind=R_GRID) lamda(im+1,im+1)
1489 real(kind=R_GRID) theta(im+1,im+1)
1490 real(kind=R_GRID) p(3,im+1,im+1)
1492 real(kind=R_GRID) rsq3, xf, y0, z0, y, x, z, ds
1493 real(kind=R_GRID) dy, dz
1498 rsq3 = 1.d0/sqrt(3.d0)
1500 y0 = rsq3; dy = -2.d0*rsq3/im
1501 z0 = -rsq3; dz = 2.d0*rsq3/im
1506 p(2,j,k) = y0 + (j-1)*dy
1507 p(3,j,k) = z0 + (k-1)*dz
1514 subroutine symm_ed(im, lamda, theta)
1517 real(kind=R_GRID) lamda(im+1,im+1)
1518 real(kind=R_GRID) theta(im+1,im+1)
1520 real(kind=R_GRID) avg
1524 lamda(i,j) = lamda(i,1)
1531 avg = 0.5d0*(lamda(i,j)-lamda(ip,j))
1532 lamda(i, j) = avg +
pi 1533 lamda(ip,j) =
pi - avg
1534 avg = 0.5d0*(theta(i,j)+theta(ip,j))
1544 avg = 0.5d0*(lamda(i,j)+lamda(i,jp))
1547 avg = 0.5d0*(theta(i,j)-theta(i,jp))
1556 real(kind=R_GRID),
intent(in):: lon, lat
1557 real(kind=R_GRID),
intent(out):: p3(3)
1558 real(kind=R_GRID) e(2)
1560 e(1) = lon; e(2) = lat
1570 real(kind=R_GRID),
intent(in) :: p(2)
1571 real(kind=R_GRID),
intent(out):: e(3)
1572 integer,
optional,
intent(in):: id
1576 real (f_p):: e1, e2, e3
1582 e1 = cos(q(2)) * cos(q(1))
1583 e2 = cos(q(2)) * sin(q(1))
1610 real(kind=R_GRID),
intent(in) :: p1(3), p2(3), p0(3)
1611 real(kind=R_GRID),
intent(out):: p(3)
1613 real(kind=R_GRID):: x1, y1, z1, x2, y2, z2, x0, y0, z0
1614 real(kind=R_GRID) nb(3)
1615 real(kind=R_GRID) pdot
1619 pdot = sqrt(nb(1)**2+nb(2)**2+nb(3)**2)
1621 nb(k) = nb(k) / pdot
1624 pdot = p0(1)*nb(1) + p0(2)*nb(2) + p0(3)*nb(3)
1626 p(k) = p0(k) - 2.d0*pdot*nb(k)
1632 subroutine mirror_latlon(lon1, lat1, lon2, lat2, lon0, lat0, lon3, lat3)
1637 real(kind=R_GRID),
intent(in):: lon1, lat1, lon2, lat2, lon0, lat0
1638 real(kind=R_GRID),
intent(out):: lon3, lat3
1640 real(kind=R_GRID) p0(3), p1(3), p2(3), nb(3), pp(3), sp(2)
1641 real(kind=R_GRID) pdot
1649 pdot = sqrt(nb(1)**2+nb(2)**2+nb(3)**2)
1651 nb(k) = nb(k) / pdot
1654 pdot = p0(1)*nb(1) + p0(2)*nb(2) + p0(3)*nb(3)
1656 pp(k) = p0(k) - 2.d0*pdot*nb(k)
1668 integer,
intent(in):: np
1669 real(kind=R_GRID),
intent(inout):: q(3,np)
1670 real(kind=R_GRID),
intent(inout):: xs(np), ys(np)
1672 real(kind=R_GRID),
parameter:: esl=1.d-10
1674 real (f_p):: dist, lat, lon
1681 dist = sqrt(p(1)**2 + p(2)**2 + p(3)**2)
1686 if ( (abs(p(1))+abs(p(2))) < esl )
then 1687 lon =
real(0.,kind=
f_p)
1689 lon = atan2( p(2), p(1) )
1692 if ( lon < 0.) lon =
real(2.,kind=
f_p)*
pi + lon
1709 real(kind=R_GRID),
intent(in) :: p1(3), p2(3)
1710 real(kind=R_GRID),
intent(out):: e(3)
1714 e(1) = p1(2)*p2(3) - p1(3)*p2(2)
1715 e(2) = p1(3)*p2(1) - p1(1)*p2(3)
1716 e(3) = p1(1)*p2(2) - p1(2)*p2(1)
1723 type(fv_grid_bounds_type),
intent(IN) :: bd
1724 integer,
intent(in):: npx, npy
1725 real(kind=R_GRID),
intent(in) :: pp(3,bd%isd:bd%ied+1,bd%jsd:bd%jed+1)
1726 real(kind=R_GRID),
intent(out):: u1(3,bd%isd:bd%ied, bd%jsd:bd%jed)
1727 real(kind=R_GRID),
intent(out):: u2(3,bd%isd:bd%ied, bd%jsd:bd%jed)
1730 real(kind=R_GRID) p1(3), p2(3), pc(3), p3(3)
1732 integer :: isd, ied, jsd, jed
1741 if ( (i<1 .and. j<1 ) .or. (i>(npx-1) .and. j<1) .or. &
1742 (i>(npx-1) .and. j>(npy-1)) .or. (i<1 .and. j>(npy-1)))
then 1748 u1(k,i,j) = pp(k,i+1,j)+pp(k,i+1,j+1) - pp(k,i,j)-pp(k,i,j+1)
1749 u2(k,i,j) = pp(k,i,j+1)+pp(k,i+1,j+1) - pp(k,i,j)-pp(k,i+1,j)
1754 call cell_center3(pp(1,i,j), pp(1,i+1,j), pp(1,i,j+1), pp(1,i+1,j+1), pc)
1776 real(kind=R_GRID),
intent(in) :: e1(2), e2(2)
1777 real(kind=R_GRID),
intent(out):: uc(3)
1779 real(kind=R_GRID),
dimension(3):: pc, p1, p2, p3
1793 real(kind=R_GRID),
intent(in) :: p1(3), p2(3)
1794 real(kind=R_GRID),
intent(out):: uc(3)
1796 real(kind=R_GRID),
dimension(3):: pc, p3
1809 real(kind=R_GRID),
intent(inout):: e(3)
1813 pdot = e(1)**2 + e(2)**2 + e(3)**2
1824 real(kind=R_GRID),
intent(in):: beta
1825 real(kind=R_GRID),
intent(in):: p1(2), p2(2)
1826 real(kind=R_GRID),
intent(out):: x_o, y_o
1828 real(kind=R_GRID):: pm(2)
1829 real(kind=R_GRID):: e1(3), e2(3), e3(3)
1830 real(kind=R_GRID):: s1, s2, s3, dd, alpha
1837 s1 = alpha*e1(1) + beta*e2(1)
1838 s2 = alpha*e1(2) + beta*e2(2)
1839 s3 = alpha*e1(3) + beta*e2(3)
1841 dd = sqrt( s1**2 + s2**2 + s3**2 )
1858 real(kind=R_GRID),
intent(in):: beta
1859 real(kind=R_GRID),
intent(in):: p1(2), p2(2)
1860 real(kind=R_GRID),
intent(out):: pb(2)
1862 real(kind=R_GRID):: pm(2)
1863 real(kind=R_GRID):: e1(3), e2(3), eb(3)
1864 real(kind=R_GRID):: dd, alpha, omg
1866 if ( abs(p1(1) - p2(1)) < 1.d-8 .and. abs(p1(2) - p2(2)) < 1.d-8)
then 1867 call mpp_error(warning,
'spherical_linear_interpolation was passed two colocated points.')
1875 dd = sqrt( e1(1)**2 + e1(2)**2 + e1(3)**2 )
1881 dd = sqrt( e2(1)**2 + e2(2)**2 + e2(3)**2 )
1889 omg = acos( e1(1)*e2(1) + e1(2)*e2(2) + e1(3)*e2(3) )
1891 if ( abs(omg) < 1.d-5 )
then 1892 print*,
'spherical_linear_interpolation: ', omg, p1, p2
1893 call mpp_error(fatal,
'spherical_linear_interpolation: interpolation not well defined between antipodal points')
1896 eb(1) = sin( beta*omg )*e2(1) + sin(alpha*omg)*e1(1)
1897 eb(2) = sin( beta*omg )*e2(2) + sin(alpha*omg)*e1(2)
1898 eb(3) = sin( beta*omg )*e2(3) + sin(alpha*omg)*e1(3)
1900 eb(1) = eb(1) / sin(omg)
1901 eb(2) = eb(2) / sin(omg)
1902 eb(3) = eb(3) / sin(omg)
1909 real(kind=R_GRID) ,
intent(IN) :: p1(2), p2(2)
1910 real(kind=R_GRID) ,
intent(OUT) :: pm(2)
1912 real(kind=R_GRID) e1(3), e2(3), e3(3)
1924 real(kind=R_GRID),
intent(IN) :: p1(3), p2(3)
1925 real(kind=R_GRID),
intent(OUT) :: e(3)
1927 real (f_p):: q1(3), q2(3)
1928 real (f_p):: dd, e1, e2, e3
1940 dd = sqrt( e1**2 + e2**2 + e3**2 )
1954 real(kind=R_GRID),
intent(IN) :: p1(2), p2(2)
1955 real(kind=R_GRID),
intent(OUT) :: e3(3)
1957 real(kind=R_GRID) e1(3), e2(3)
1968 real(kind=R_GRID),
intent(IN) :: q1(2), q2(2)
1969 real(kind=R_GRID),
intent(IN),
optional ::
radius 1971 real (f_p):: p1(2), p2(2)
1980 beta = asin( sqrt( sin((p1(2)-p2(2))/2.)**2 + cos(p1(2))*cos(p2(2))* &
1981 sin((p1(1)-p2(1))/2.)**2 ) ) * 2.
1983 if (
present(
radius) )
then 2000 real(kind=R_GRID),
dimension(3),
intent(in) :: v1, v2
2001 real(kind=R_GRID),
intent(IN),
optional ::
radius 2002 real(kind=R_GRID) :: norm
2004 norm = (v1(1)*v1(1)+v1(2)*v1(2)+v1(3)*v1(3)) &
2005 *(v2(1)*v2(1)+v2(2)*v2(2)+v2(3)*v2(3))
2014 if (
present(
radius) )
then 2023 subroutine intersect(a1,a2,b1,b2,radius,x_inter,local_a,local_b)
2043 real(kind=R_GRID),
dimension(3),
intent(in) :: a1, a2, b1, b2
2044 real(kind=R_GRID),
intent(in) :: radius
2045 real(kind=R_GRID),
dimension(3),
intent(out) :: x_inter
2046 logical,
intent(out) :: local_a,local_b
2050 real(kind=R_GRID) :: a2_xy, b1_xy, b2_xy, a2_xz, b1_xz, b2_xz, &
2051 b1_xyz, b2_xyz, length
2055 a2_xy=a2(1)*a1(2)-a2(2)*a1(1)
2056 b1_xy=b1(1)*a1(2)-b1(2)*a1(1)
2057 b2_xy=b2(1)*a1(2)-b2(2)*a1(1)
2059 a2_xz=a2(1)*a1(3)-a2(3)*a1(1)
2060 b1_xz=b1(1)*a1(3)-b1(3)*a1(1)
2061 b2_xz=b2(1)*a1(3)-b2(3)*a1(1)
2063 b1_xyz=b1_xy*a2_xz-b1_xz*a2_xy
2064 b2_xyz=b2_xy*a2_xz-b2_xz*a2_xy
2066 if (b1_xyz==0.0d0)
then 2068 elseif (b2_xyz==0.0d0)
then 2071 x_inter(:)=b2(:)-b1(:)*b2_xyz/b1_xyz
2072 length=sqrt(x_inter(1)*x_inter(1)+x_inter(2)*x_inter(2)+x_inter(3)*x_inter(3))
2073 x_inter(:)=radius/length*x_inter(:)
2085 real(kind=R_GRID),
dimension(3) :: center, dx
2086 real(kind=R_GRID) :: dist1,dist2
2088 center(:)=0.25*(a1(:)+a2(:)+b1(:)+b2(:))
2089 dx(:)=+x_inter(:)-center(:)
2090 dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2091 dx(:)=-x_inter(:)-center(:)
2092 dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2094 if (dist2<dist1) x_inter(:)=-x_inter(:)
2099 real(kind=R_GRID),
dimension(3),
intent(in) :: x1,x2
2100 logical,
intent(out) :: local
2102 real(kind=R_GRID),
dimension(3) :: dx
2103 real(kind=R_GRID) :: dist, dist1, dist2
2106 dist=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2108 dx(:)=x1(:)-x_inter(:)
2109 dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2110 dx(:)=x2(:)-x_inter(:)
2111 dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2113 if (dist1<=dist .and. dist2<=dist)
then 2123 subroutine intersect_cross(a1,a2,b1,b2,radius,x_inter,local_a,local_b)
2137 real(kind=R_GRID),
dimension(3),
intent(in) :: a1, a2, b1, b2
2138 real(kind=R_GRID),
intent(in) :: radius
2139 real(kind=R_GRID),
dimension(3),
intent(out) :: x_inter
2140 logical,
intent(out) :: local_a,local_b
2141 real(kind=R_GRID),
dimension(3) :: v1, v2
2159 v1 = v1/sqrt(v1(1)**2 + v1(2)**2 + v1(3)**2)
2160 v2 = v2/sqrt(v2(1)**2 + v2(2)**2 + v2(3)**2)
2164 x_inter = x_inter/sqrt(x_inter(1)**2 + x_inter(2)**2 + x_inter(3)**2)
2173 real(kind=R_GRID),
dimension(3) :: center, dx
2174 real(kind=R_GRID) :: dist1,dist2
2176 center(:)=0.25*(a1(:)+a2(:)+b1(:)+b2(:))
2177 dx(:)=+x_inter(:)-center(:)
2178 dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2179 dx(:)=-x_inter(:)-center(:)
2180 dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2182 if (dist2<dist1) x_inter(:)=-x_inter(:)
2187 real(kind=R_GRID),
dimension(3),
intent(in) :: x1,x2
2188 logical,
intent(out) :: local
2190 real(kind=R_GRID),
dimension(3) :: dx
2191 real(kind=R_GRID) :: dist, dist1, dist2
2194 dist=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2196 dx(:)=x1(:)-x_inter(:)
2197 dist1=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2198 dx(:)=x2(:)-x_inter(:)
2199 dist2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
2201 if (dist1<=dist .and. dist2<=dist)
then 2214 real(kind=R_GRID),
intent(IN) :: pp(2)
2215 real(kind=R_GRID),
intent(OUT) :: elon(3), elat(3)
2217 real (f_p):: lon, lat
2218 real (f_p):: sin_lon, cos_lon, sin_lat, cos_lat
2232 elat(1) = -sin_lat*cos_lon
2233 elat(2) = -sin_lat*sin_lon
2240 real(kind=R_GRID) function v_prod(v1, v2)
2241 real(kind=R_GRID) v1(3), v2(3)
2243 v_prod = v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)
2249 type(fv_grid_bounds_type),
intent(IN) :: bd
2250 logical,
intent(in):: hydrostatic
2251 real(kind=R_GRID),
intent(in) :: agrid(bd%isd:bd%ied,bd%jsd:bd%jed,2)
2252 integer,
intent(in) :: grid_type
2253 integer,
intent(in) :: ord
2254 type(fv_grid_type),
intent(INOUT),
target :: gridstruct
2257 integer :: is, ie, js, je
2260 real,
pointer,
dimension(:,:) :: a11, a12, a21, a22
2261 real,
pointer,
dimension(:,:) :: z11, z12, z21, z22
2262 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: vlon, vlat
2263 real(kind=R_GRID),
pointer,
dimension(:,:,:) :: ee1, ee2, ec1, ec2
2272 vlon => gridstruct%vlon
2273 vlat => gridstruct%vlat
2274 a11 => gridstruct%a11
2275 a12 => gridstruct%a12
2276 a21 => gridstruct%a21
2277 a22 => gridstruct%a22
2278 z11 => gridstruct%z11
2279 z12 => gridstruct%z12
2280 z21 => gridstruct%z21
2281 z22 => gridstruct%z22
2282 ee1 => gridstruct%ee1
2283 ee2 => gridstruct%ee2
2284 ec1 => gridstruct%ec1
2285 ec2 => gridstruct%ec2
2295 z11(i,j) =
v_prod(ec1(1:3,i,j), vlon(i,j,1:3))
2296 z12(i,j) =
v_prod(ec1(1:3,i,j), vlat(i,j,1:3))
2297 z21(i,j) =
v_prod(ec2(1:3,i,j), vlon(i,j,1:3))
2298 z22(i,j) =
v_prod(ec2(1:3,i,j), vlat(i,j,1:3))
2300 a11(i,j) = 0.5d0*z22(i,j) / gridstruct%sin_sg(i,j,5)
2301 a12(i,j) = -0.5d0*z12(i,j) / gridstruct%sin_sg(i,j,5)
2302 a21(i,j) = -0.5d0*z21(i,j) / gridstruct%sin_sg(i,j,5)
2303 a22(i,j) = 0.5d0*z11(i,j) / gridstruct%sin_sg(i,j,5)
2313 subroutine cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
2315 integer,
intent(in) :: km, npx, npy,
grid_type, c2l_ord
2316 integer,
intent(in) :: mode
2318 real,
intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2319 real,
intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2320 real,
intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2321 real,
intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2322 type(
domain2d),
intent(INOUT) :: domain
2323 logical,
intent(IN) :: nested
2325 if ( c2l_ord == 2 )
then 2328 call c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km,
grid_type, domain, nested, mode, bd)
2334 subroutine c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km, grid_type, domain, nested, mode, bd)
2336 type(fv_grid_bounds_type),
intent(IN) :: bd
2337 integer,
intent(in) :: km, npx, npy, grid_type
2338 integer,
intent(in):: mode
2339 type(fv_grid_type),
intent(IN),
target :: gridstruct
2340 real,
intent(inout):: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2341 real,
intent(inout):: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2342 real,
intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2343 real,
intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2344 type(domain2d),
intent(INOUT) :: domain
2345 logical,
intent(IN) :: nested
2349 real :: a2 = -0.0625
2352 real utmp(bd%is:bd%ie, bd%js:bd%je+1)
2353 real vtmp(bd%is:bd%ie+1,bd%js:bd%je)
2354 real wu(bd%is:bd%ie, bd%js:bd%je+1)
2355 real wv(bd%is:bd%ie+1,bd%js:bd%je)
2358 integer :: is, ie, js, je
2366 if ( mode > 0 )
then 2378 do j=
max(1,js),
min(npy-1,je)
2379 do i=
max(1,is),
min(npx-1,ie)
2380 utmp(i,j) = c2*(u(i,j-1,k)+u(i,j+2,k)) + c1*(u(i,j,k)+u(i,j+1,k))
2381 vtmp(i,j) = c2*(v(i-1,j,k)+v(i+2,j,k)) + c1*(v(i,j,k)+v(i+1,j,k))
2385 do j=
max(2,js),
min(npy-2,je)
2386 do i=
max(2,is),
min(npx-2,ie)
2387 utmp(i,j) = c2*(u(i,j-1,k)+u(i,j+2,k)) + c1*(u(i,j,k)+u(i,j+1,k))
2388 vtmp(i,j) = c2*(v(i-1,j,k)+v(i+2,j,k)) + c1*(v(i,j,k)+v(i+1,j,k))
2394 wv(i,1) = v(i,1,k)*gridstruct%dy(i,1)
2397 vtmp(i,1) = 2.*(wv(i,1) + wv(i+1,1)) / (gridstruct%dy(i,1)+gridstruct%dy(i+1,1))
2398 utmp(i,1) = 2.*(u(i,1,k)*gridstruct%dx(i,1) + u(i,2,k)*gridstruct%dx(i,2)) &
2399 / ( gridstruct%dx(i,1) + gridstruct%dx(i,2))
2405 if ( (je+1)==npy )
then 2408 wv(i,j) = v(i,j,k)*gridstruct%dy(i,j)
2411 vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j)) / (gridstruct%dy(i,j)+gridstruct%dy(i+1,j))
2412 utmp(i,j) = 2.*(u(i,j,k)*gridstruct%dx(i,j) + u(i,j+1,k)*gridstruct%dx(i,j+1)) &
2413 / ( gridstruct%dx(i,j) + gridstruct%dx(i,j+1))
2422 wv(1,j) = v(1,j,k)*gridstruct%dy(1,j)
2423 wv(2,j) = v(2,j,k)*gridstruct%dy(2,j)
2426 wu(i,j) = u(i,j,k)*gridstruct%dx(i,j)
2429 utmp(i,j) = 2.*(wu(i,j) + wu(i,j+1))/(gridstruct%dx(i,j)+gridstruct%dx(i,j+1))
2430 vtmp(i,j) = 2.*(wv(1,j) + wv(2,j ))/(gridstruct%dy(1,j)+gridstruct%dy(2,j))
2436 if ( (ie+1)==npx)
then 2439 wv(i, j) = v(i, j,k)*gridstruct%dy(i, j)
2440 wv(i+1,j) = v(i+1,j,k)*gridstruct%dy(i+1,j)
2443 wu(i,j) = u(i,j,k)*gridstruct%dx(i,j)
2446 utmp(i,j) = 2.*(wu(i,j) + wu(i, j+1))/(gridstruct%dx(i,j)+gridstruct%dx(i,j+1))
2447 vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j ))/(gridstruct%dy(i,j)+gridstruct%dy(i+1,j))
2458 ua(i,j,k) = gridstruct%a11(i,j)*utmp(i,j) + gridstruct%a12(i,j)*vtmp(i,j)
2459 va(i,j,k) = gridstruct%a21(i,j)*utmp(i,j) + gridstruct%a22(i,j)*vtmp(i,j)
2466 ua(i,j,k) = a2*(u(i,j-1,k)+u(i,j+2,k)) + a1*(u(i,j,k)+u(i,j+1,k))
2467 va(i,j,k) = a2*(v(i-1,j,k)+v(i+2,j,k)) + a1*(v(i,j,k)+v(i+1,j,k))
2474 subroutine c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, do_halo)
2477 real,
intent(in) :: u(bd%isd:bd%ied,bd%jsd:bd%jed+1,km)
2478 real,
intent(in) :: v(bd%isd:bd%ied+1,bd%jsd:bd%jed,km)
2480 logical,
intent(in) :: do_halo
2482 real,
intent(out):: ua(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2483 real,
intent(out):: va(bd%isd:bd%ied, bd%jsd:bd%jed,km)
2486 real wu(bd%is-1:bd%ie+1, bd%js-1:bd%je+2)
2487 real wv(bd%is-1:bd%ie+2, bd%js-1:bd%je+1)
2488 real u1(bd%is-1:bd%ie+1), v1(bd%is-1:bd%ie+1)
2490 integer :: is, ie, js, je
2492 real,
dimension(:,:),
pointer :: a11, a12, a21, a22
2493 real,
dimension(:,:),
pointer :: dx, dy, rdxa, rdya
2495 a11 => gridstruct%a11
2496 a12 => gridstruct%a12
2497 a21 => gridstruct%a21
2498 a22 => gridstruct%a22
2502 rdxa => gridstruct%rdxa
2503 rdya => gridstruct%rdya
2523 wu(i,j) = u(i,j,k)*dx(i,j)
2528 wv(i,j) = v(i,j,k)*dy(i,j)
2535 u1(i) = 2.*(wu(i,j) + wu(i,j+1)) / (dx(i,j)+dx(i,j+1))
2536 v1(i) = 2.*(wv(i,j) + wv(i+1,j)) / (dy(i,j)+dy(i+1,j))
2540 ua(i,j,k) = a11(i,j)*u1(i) + a12(i,j)*v1(i)
2541 va(i,j,k) = a21(i,j)*u1(i) + a22(i,j)*v1(i)
2548 ua(i,j,k) = 0.5*(u(i,j,k)+u(i, j+1,k))
2549 va(i,j,k) = 0.5*(v(i,j,k)+v(i+1,j, k))
2558 subroutine expand_cell(q1, q2, q3, q4, a1, a2, a3, a4, fac)
2565 real(kind=R_GRID),
intent(in):: q1(2), q2(2), q3(2), q4(2)
2566 real(kind=R_GRID),
intent(in):: fac
2569 real(kind=R_GRID),
intent(out):: a1(2), a2(2), a3(2), a4(2)
2571 real(kind=R_GRID) qq1(3), qq2(3), qq3(3), qq4(3)
2572 real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3)
2573 real(kind=R_GRID) ec(3)
2574 real(f_p):: dd, d1, d2, d3, d4
2585 ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2587 dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2595 qq1(k) = ec(k) + fac*(p1(k)-ec(k))
2596 qq2(k) = ec(k) + fac*(p2(k)-ec(k))
2597 qq3(k) = ec(k) + fac*(p3(k)-ec(k))
2598 qq4(k) = ec(k) + fac*(p4(k)-ec(k))
2604 d1 = sqrt( qq1(1)**2 + qq1(2)**2 + qq1(3)**2 )
2605 d2 = sqrt( qq2(1)**2 + qq2(2)**2 + qq2(3)**2 )
2606 d3 = sqrt( qq3(1)**2 + qq3(2)**2 + qq3(3)**2 )
2607 d4 = sqrt( qq4(1)**2 + qq4(2)**2 + qq4(3)**2 )
2609 qq1(k) = qq1(k) / d1
2610 qq2(k) = qq2(k) / d2
2611 qq3(k) = qq3(k) / d3
2612 qq4(k) = qq4(k) / d4
2628 real(kind=R_GRID) ,
intent(in ) :: q1(2), q2(2), q3(2), q4(2)
2629 real(kind=R_GRID) ,
intent(out) :: e2(2)
2631 real(kind=R_GRID) p1(3), p2(3), p3(3), p4(3)
2632 real(kind=R_GRID) ec(3)
2633 real(kind=R_GRID) dd
2642 ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2644 dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2657 real(kind=R_GRID) ,
intent(IN) :: p1(3), p2(3), p3(3), p4(3)
2658 real(kind=R_GRID) ,
intent(OUT) :: ec(3)
2664 ec(k) = p1(k) + p2(k) + p3(k) + p4(k)
2666 dd = sqrt( ec(1)**2 + ec(2)**2 + ec(3)**2 )
2676 real(kind=R_GRID) function get_area(p1, p4, p2, p3, radius)
2678 real(kind=R_GRID),
intent(in),
dimension(2):: p1, p2, p3, p4
2679 real(kind=R_GRID),
intent(in),
optional::
radius 2681 real(kind=R_GRID) e1(3), e2(3), e3(3)
2682 real(kind=R_GRID) ang1, ang2, ang3, ang4
2711 if (
present(
radius) )
then 2729 real(kind=R_GRID),
dimension(3),
intent(in) :: v1, v2, point
2731 real(kind=R_GRID) :: angle, side
2743 real(kind=R_GRID),
dimension(2),
intent(in) :: v1, v2, point
2745 real(kind=R_GRID),
dimension(3) :: c1, c2, cpoint
2747 real(kind=R_GRID) :: angle,side
2775 real(kind=R_GRID) p1(3), p2(3), p3(3)
2777 real (f_p):: e1(3), e2(3), e3(3)
2778 real (f_p):: px, py, pz
2779 real (f_p):: qx, qy, qz
2780 real (f_p):: angle, ddd
2793 px = e1(2)*e2(3) - e1(3)*e2(2)
2794 py = e1(3)*e2(1) - e1(1)*e2(3)
2795 pz = e1(1)*e2(2) - e1(2)*e2(1)
2797 qx = e1(2)*e3(3) - e1(3)*e3(2)
2798 qy = e1(3)*e3(1) - e1(1)*e3(3)
2799 qz = e1(1)*e3(2) - e1(2)*e3(1)
2801 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz)
2803 if ( ddd <= 0.0d0 )
then 2806 ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd)
2807 if ( abs(ddd)>1.d0)
then 2808 angle = 2.d0*atan(1.0)
2810 if (ddd < 0.d0)
then 2811 angle = 4.d0*atan(1.0d0)
2825 real(kind=R_GRID) function cos_angle(p1, p2, p3)
2833 real(kind=R_GRID),
intent(in):: p1(3), p2(3), p3(3)
2835 real (f_p):: e1(3), e2(3), e3(3)
2836 real (f_p):: px, py, pz
2837 real (f_p):: qx, qy, qz
2838 real (f_p):: angle, ddd
2851 px = e1(2)*e2(3) - e1(3)*e2(2)
2852 py = e1(3)*e2(1) - e1(1)*e2(3)
2853 pz = e1(1)*e2(2) - e1(2)*e2(1)
2856 qx = e1(2)*e3(3) - e1(3)*e3(2)
2857 qy = e1(3)*e3(1) - e1(1)*e3(3)
2858 qz = e1(1)*e3(2) - e1(2)*e3(1)
2861 ddd = sqrt( (px**2+py**2+pz**2)*(qx**2+qy**2+qz**2) )
2862 if ( ddd > 0.d0 )
then 2863 angle = (px*qx+py*qy+pz*qz) / ddd
2873 real function g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
2875 integer,
intent(IN) :: ifirst, ilast
2876 integer,
intent(IN) :: jfirst, jlast, ngc
2877 integer,
intent(IN) :: mode
2878 logical,
intent(in),
optional :: reproduce
2879 real,
intent(IN) :: p(ifirst:ilast,jfirst:jlast)
2880 real(kind=R_GRID),
intent(IN) :: area(ifirst-ngc:ilast+ngc,jfirst-ngc:jlast+ngc)
2881 type(
domain2d),
intent(IN) :: domain
2884 logical,
SAVE :: g_sum_initialized = .false.
2885 real(kind=R_GRID),
SAVE :: global_area
2886 real :: tmp(ifirst:ilast,jfirst:jlast)
2888 if ( .not. g_sum_initialized )
then 2889 global_area =
mpp_global_sum(domain, area, flags=bitwise_efp_sum)
2890 if ( is_master() )
write(*,*)
'Global Area=',global_area
2891 g_sum_initialized = .true.
2897 if (
present(reproduce) )
then 2899 gsum =
mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast), &
2900 flags=bitwise_efp_sum)
2902 gsum =
mpp_global_sum(domain, p(:,:)*area(ifirst:ilast,jfirst:jlast))
2911 gsum = gsum + p(i,j)*area(i,j)
2914 call mp_reduce_sum(gsum)
2918 g_sum = gsum / global_area
2926 real function global_qsum(p, ifirst, ilast, jfirst, jlast)
2928 integer,
intent(IN) :: ifirst, ilast
2929 integer,
intent(IN) :: jfirst, jlast
2930 real,
intent(IN) :: p(ifirst:ilast,jfirst:jlast)
2937 gsum = gsum + p(i,j)
2940 call mp_reduce_sum(gsum)
2947 subroutine global_mx(q, n_g, qmin, qmax, bd)
2950 integer,
intent(in):: n_g
2951 real(kind=R_GRID),
intent(in)::q(bd%is-n_g:bd%ie+n_g, bd%js-n_g:bd%je+n_g)
2952 real(kind=R_GRID),
intent(out):: qmin, qmax
2955 integer :: is, ie, js, je
2966 qmin =
min(qmin, q(i,j))
2967 qmax =
max(qmax, q(i,j))
2970 call mp_reduce_min(qmin)
2971 call mp_reduce_max(qmax)
2975 subroutine global_mx_c(q, i1, i2, j1, j2, qmin, qmax)
2977 integer,
intent(in):: i1, i2, j1, j2
2978 real(kind=R_GRID),
intent(in) :: q(i1:i2,j1:j2)
2979 real(kind=R_GRID),
intent(out) :: qmin, qmax
2986 qmin =
min(qmin, q(i,j))
2987 qmax =
max(qmax, q(i,j))
2990 call mp_reduce_min(qmin)
2991 call mp_reduce_max(qmax)
2997 subroutine fill_ghost_r4(q, npx, npy, value, bd)
2999 type(fv_grid_bounds_type),
intent(IN) :: bd
3000 real(kind=4),
intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
3001 integer,
intent(in):: npx, npy
3002 real,
intent(in):: value
3005 integer :: is, ie, js, je
3006 integer :: isd, ied, jsd, jed
3019 if ( (i<1 .and. j<1) )
then 3022 if ( i>(npx-1) .and. j<1 )
then 3025 if ( i>(npx-1) .and. j>(npy-1) )
then 3028 if ( i<1 .and. j>(npy-1) )
then 3034 end subroutine fill_ghost_r4
3039 type(fv_grid_bounds_type),
intent(IN) :: bd
3040 real(kind=R_GRID),
intent(inout):: q(bd%isd:bd%ied,bd%jsd:bd%jed)
3041 integer,
intent(in):: npx, npy
3042 real,
intent(in):: value
3045 integer :: is, ie, js, je
3046 integer :: isd, ied, jsd, jed
3059 if ( (i<1 .and. j<1) )
then 3062 if ( i>(npx-1) .and. j<1 )
then 3065 if ( i>(npx-1) .and. j>(npy-1) )
then 3068 if ( i<1 .and. j>(npy-1) )
then 3078 subroutine make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
3080 integer,
intent(in ):: km
3081 integer,
intent(out):: kks
3082 real(kind=R_GRID),
intent(in):: area(bd%isd:bd%ied,bd%jsd:bd%jed)
3083 real,
intent(INOUT) :: ptop
3084 real,
intent(inout):: pe(bd%is-1:bd%ie+1,km+1,bd%js-1:bd%je+1)
3085 real,
intent(out):: ak(km+1), bk(km+1)
3086 type(
domain2d),
intent(IN) :: domain
3089 real,
allocatable:: pem(:,:)
3092 integer :: is, ie, js, je
3099 allocate ( pem(is:ie,js:je) )
3105 pem(i,j) = pe(i,k,j)
3110 p4 =
g_sum(domain, pem, is, ie, js, je, ng, area, 1)
3129 bk(k) = (ph(k) - ph(1)) / (ph(km+1)-ph(1))
3130 ak(k) = ph(1)*(1.-bk(k))
3133 if ( is_master() )
then 3134 write(*,*)
'Make_eta_level ...., ptop=', ptop
3137 write(*,*) ph(k), ak(k), bk(k)
3147 integer,
intent (in) :: n
3149 real(kind=R_GRID),
intent (inout),
dimension (n,n):: a
3150 real(kind=R_GRID),
intent (out),
dimension (n,n):: x
3151 real(kind=R_GRID),
dimension (n,n) :: b
3164 call elgs (a,n,indx)
3169 b(indx(j),k) = b(indx(j),k) - a(indx(j),i)*b(indx(i),k)
3175 x(n,i) = b(indx(n),i)/a(indx(n),n)
3177 x(j,i) = b(indx(j),i)
3179 x(j,i) = x(j,i)-a(indx(j),k)*x(k,i)
3181 x(j,i) = x(j,i)/a(indx(j),j)
3188 subroutine elgs (a,n,indx)
3196 integer,
intent (in) :: n
3197 integer :: i,j,k,itmp
3198 integer,
intent (out),
dimension (n) :: indx
3199 real(kind=R_GRID),
intent (inout),
dimension (n,n) :: a
3201 real(kind=R_GRID) :: c1, pie, pi1, pj
3202 real(kind=R_GRID),
dimension (n) :: c
3213 c1 =
max(c1,abs(a(i,j)))
3223 pie = abs(a(indx(i),j))/c(indx(i))
3236 pj = a(indx(i),j)/a(indx(j),j)
3245 a(indx(i),k) = a(indx(i),k)-pj*a(indx(j),k)
3253 real(kind=R_GRID),
intent(IN) :: pp(2)
3254 real(kind=R_GRID),
intent(OUT) :: elon(3), elat(3)
3256 elon(1) = -sin(pp(1))
3257 elon(2) = cos(pp(1))
3259 elat(1) = -sin(pp(2))*cos(pp(1))
3260 elat(2) = -sin(pp(2))*sin(pp(1))
3262 elat(3) = cos(pp(2))
3273 integer,
intent(in):: np
3274 real(kind=R_GRID),
intent(in):: e(3,np)
3275 real(kind=R_GRID),
intent(inout):: f(3,np)
3281 ap = f(1,i)*e(1,i) + f(2,i)*e(2,i) + f(3,i)*e(3,i)
3282 f(1,i) = f(1,i) - ap*e(1,i)
3283 f(2,i) = f(2,i) - ap*e(2,i)
3284 f(3,i) = f(3,i) - ap*e(3,i)
3290 subroutine init_mq(phis, gridstruct, npx, npy, is, ie, js, je, ng)
3291 integer,
intent(in):: npx, npy, is, ie, js, je, ng
3292 real,
intent(in):: phis(is-ng:ie+ng, js-ng:je+ng)
3296 real zs(is-ng:ie+ng, js-ng:je+ng)
3297 real zb(is-ng:ie+ng, js-ng:je+ng)
3298 real pdx(3,is:ie,js:je+1)
3299 real pdy(3,is:ie+1,js:je)
3302 real,
pointer :: rarea(:,:)
3303 real,
pointer,
dimension(:,:) :: dx, dy
3304 real,
pointer,
dimension(:,:,:) :: en1, en2, agrid, vlon, vlat
3306 rarea => gridstruct%rarea
3309 en1 => gridstruct%en1
3310 en2 => gridstruct%en2
3311 agrid => gridstruct%agrid
3312 vlon => gridstruct%vlon
3313 vlat => gridstruct%vlat
3319 zs(i,j) = phis(i,j) / grav
3325 call a2b_ord4(zs, zb, gridstruct, npx, npy, is, ie, js, je, ng)
3330 pdx(n,i,j) = 0.5*(zb(i,j)+zb(i+1,j))*dx(i,j)*en1(n,i,j)
3337 pdy(n,i,j) = 0.5*(zb(i,j)+zb(i,j+1))*dy(i,j)*en2(n,i,j)
3345 idiag%zxg(i,j) = vlon(i,j,1)*(pdx(1,i,j+1)-pdx(1,i,j)-pdy(1,i,j)+pdy(1,i+1,j)) &
3346 + vlon(i,j,2)*(pdx(2,i,j+1)-pdx(2,i,j)-pdy(2,i,j)+pdy(2,i+1,j)) &
3347 + vlon(i,j,3)*(pdx(3,i,j+1)-pdx(3,i,j)-pdy(3,i,j)+pdy(3,i+1,j))
3350 idiag%zxg(i,j) = idiag%zxg(i,j)*rarea(i,j) *
radius*cos(agrid(i,j,2))
3355 end subroutine init_mq
subroutine cell_center3(p1, p2, p3, p4, ec)
real, parameter, public radius
Radius of the Earth [m].
subroutine, public grid_utils_end
subroutine elgs(a, n, indx)
real(kind=r_grid) function, public cos_angle(p1, p2, p3)
subroutine, public mid_pt_cart(p1, p2, e3)
real, parameter, public omega
Rotation rate of the Earth [1/s].
subroutine, public get_unit_vect2(e1, e2, uc)
real, parameter, public ptop_min
subroutine, public spherical_linear_interpolation(beta, p1, p2, pb)
real(kind=8), parameter, public pi_8
Ratio of circle circumference to diameter [N/A].
subroutine, public cell_center2(q1, q2, q3, q4, e2)
subroutine mirror_latlon(lon1, lat1, lon2, lat2, lon0, lat0, lon3, lat3)
subroutine latlon2xyz2(lon, lat, p3)
subroutine, public set_eta(km, ks, ptop, ak, bk)
integer, parameter, public ip
real(kind=r_grid) function, public v_prod(v1, v2)
integer, parameter, public corner
subroutine get_center_vect(npx, npy, pp, u1, u2, bd)
void mid_pt3_cart(const double *p1, const double *p2, double *e)
subroutine, public global_mx(q, n_g, qmin, qmax, bd)
subroutine mirror_xyz(p1, p2, p0, p)
subroutine, public gnomonic_grids(grid_type, im, lon, lat)
real, parameter tiny_number
void mid_pt_sphere(const double *p1, const double *p2, double *pm)
void vect_cross(const double *p1, const double *p2, double *e)
subroutine edge_factors(edge_s, edge_n, edge_w, edge_e, non_ortho, grid, agrid, npx, npy, bd)
double spherical_angle(const double *v1, const double *v2, const double *v3)
subroutine, public cubed_to_latlon(u, v, ua, va, gridstruct, npx, npy, km, mode, grid_type, domain, nested, c2l_ord, bd)
subroutine gnomonic_dist(im, lamda, theta)
subroutine invert_matrix(n, a, x)
subroutine efactor_a2c_v(edge_vect_s, edge_vect_n, edge_vect_w, edge_vect_e, non_ortho, grid, agrid, npx, npy, nested, bd)
subroutine, public project_sphere_v(np, f, e)
subroutine get_unit_vect3(p1, p2, uc)
real(kind=r_grid) function dist2side(v1, v2, point)
subroutine global_mx_c(q, i1, i2, j1, j2, qmin, qmax)
integer, parameter, public agrid
subroutine gnomonic_ed_limited(im, in, nghost, lL, lR, uL, uR, lamda, theta)
real(kind=r_grid) function, public get_area(p1, p4, p2, p3, radius)
real, parameter, public big_number
integer, parameter, public f_p
subroutine intersect_cross(a1, a2, b1, b2, radius, x_inter, local_a, local_b)
subroutine, public cart_to_latlon(np, q, xs, ys)
subroutine, public c2l_ord2(u, v, ua, va, gridstruct, km, grid_type, bd, do_halo)
subroutine, public expand_cell(q1, q2, q3, q4, a1, a2, a3, a4, fac)
integer, parameter, public ng
subroutine timing_on(blk_name)
real(kind=r_grid) function, public dist2side_latlon(v1, v2, point)
integer, parameter, public r_grid
real function, public great_circle_dist(q1, q2, radius)
subroutine, public get_latlon_vector(pp, elon, elat)
subroutine, public direct_transform(c, i1, i2, j1, j2, lon_p, lat_p, n, lon, lat)
subroutine, public make_eta_level(km, pe, area, kks, ak, bk, ptop, domain, bd)
real function, public g_sum(domain, p, ifirst, ilast, jfirst, jlast, ngc, area, mode, reproduce)
subroutine, public intp_great_circle(beta, p1, p2, x_o, y_o)
void normalize_vect(double *e)
subroutine, public grid_utils_init(Atm, npx, npy, npz, non_ortho, grid_type, c2l_order)
subroutine fill_ghost_r8(q, npx, npy, value, bd)
real function, public inner_prod(v1, v2)
subroutine intersect(a1, a2, b1, b2, radius, x_inter, local_a, local_b)
subroutine symm_ed(im, lamda, theta)
subroutine check_local(x1, x2, local)
subroutine init_cubed_to_latlon(gridstruct, hydrostatic, agrid, grid_type, ord, bd)
void latlon2xyz(int size, const double *lon, const double *lat, double *x, double *y, double *z)
integer, parameter, public scalar_pair
subroutine c2l_ord4(u, v, ua, va, gridstruct, npx, npy, km, grid_type, domain, nested, mode, bd)
real function, public global_qsum(p, ifirst, ilast, jfirst, jlast)
subroutine gnomonic_ed(im, lamda, theta)
subroutine gnomonic_angl(im, lamda, theta)
real(kind=r_grid) function great_circle_dist_cart(v1, v2, radius)
integer, parameter, public cgrid_ne
Derived type containing the data.
logical, public symm_grid
real(fp), parameter, public pi
subroutine timing_off(blk_name)
void unit_vect_latlon(int size, const double *lon, const double *lat, double *vlon, double *vlat)