15 real,
parameter ::
pi=3.1415926535897931
16 real,
parameter ::
rad2deg = 180./3.14159265358979
17 real,
parameter ::
deg2rad = 3.14159265358979/180.
49 subroutine minmax(im,jm,a,title,imax,jmax)
53 character(len=8),
intent(in) :: title
55 integer,
intent(in) :: im, jm
56 integer,
intent(out),
optional :: imax, jmax
58 real,
intent(in) :: a(im,jm)
67 if (
present(imax) .and.
present(jmax))
then
74 if(a(i,j) >= rmax)
then
76 if (
present(imax) .and.
present(jmax))
then
81 if(a(i,j) <= rmin)rmin=a(i,j)
85 write(6,150) title,rmin,rmax
86150
format(
' - ',a8,
' MIN=',e13.4,2x,
'MAX=',e13.4)
104 integer,
intent(in) :: siz
105 real,
intent(in) :: lon(siz), lat(siz)
106 real,
intent(out) :: x(siz), y(siz), z(siz)
111 x(n) = cos(lat(n))*cos(lon(n))
112 y(n) = cos(lat(n))*sin(lon(n))
131 real,
intent(in) :: dy
155 real,
intent(in) :: dx, lat_in
160 if (lat > 89.5) lat = 89.5
161 if (lat < -89.5) lat = -89.5
180 integer,
intent(in) :: imn, jmn
181 integer(1),
intent(inout) :: mask(imn,jmn)
183 integer :: i, j, it, jt
223 integer,
intent(in) :: imn, jmn
224 integer(2),
intent(inout) :: glob(imn,jmn)
226 integer :: i, j, it, jt
267 real,
parameter :: epsln30 = 1.e-30
269 real,
intent(in) :: v1(3), v2(3), v3(3)
271 real :: px, py, pz, qx, qy, qz, ddd
275 px = v1(2)*v2(3) - v1(3)*v2(2)
276 py = v1(3)*v2(1) - v1(1)*v2(3)
277 pz = v1(1)*v2(2) - v1(2)*v2(1)
281 qx = v1(2)*v3(3) - v1(3)*v3(2);
282 qy = v1(3)*v3(1) - v1(1)*v3(3);
283 qz = v1(1)*v3(2) - v1(2)*v3(1);
285 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
286 if ( ddd <= 0.0 )
then
289 ddd = (px*qx+py*qy+pz*qz) / sqrt(ddd);
290 if( abs(ddd-1) < epsln30 ) ddd = 1;
291 if( abs(ddd+1) < epsln30 ) ddd = -1;
292 if ( ddd>1. .or. ddd<-1. )
then
323 real,
parameter :: epsln10 = 1.e-10
324 real,
parameter :: epsln8 = 1.e-8
325 real,
parameter :: range_check_criteria=0.05
327 integer,
intent(in) :: npts
329 real,
intent(in) :: lon1, lat1
330 real,
intent(in) :: lon2(npts), lat2(npts)
334 real :: anglesum, angle
335 real :: x2(npts), y2(npts), z2(npts)
336 real :: lon1_1d(1), lat1_1d(1)
337 real :: x1(1), y1(1), z1(1)
338 real :: pnt0(3),pnt1(3),pnt2(3)
339 real :: max_x2,min_x2,max_y2,min_y2,max_z2,min_z2
346 call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
349 if( x1(1) > max_x2+range_check_criteria )
return
351 if( x1(1)+range_check_criteria < min_x2 )
return
353 if( y1(1) > max_y2+range_check_criteria )
return
355 if( y1(1)+range_check_criteria < min_y2 )
return
357 if( z1(1) > max_z2+range_check_criteria )
return
359 if( z1(1)+range_check_criteria < min_z2 )
return
368 if(abs(x1(1)-x2(i)) < epsln10 .and. &
369 abs(y1(1)-y2(i)) < epsln10 .and. &
370 abs(z1(1)-z2(i)) < epsln10 )
then
383 anglesum = anglesum + angle
386 if(abs(anglesum-2*
pi) < epsln8)
then
409 subroutine find_poles(geolat, nx, ny, i_north_pole, j_north_pole, &
410 i_south_pole, j_south_pole)
414 integer,
intent(in) :: nx, ny
416 real,
intent(in) :: geolat(nx+1,ny+1)
418 integer,
intent(out) :: i_north_pole, j_north_pole
419 integer,
intent(out) :: i_south_pole, j_south_pole
423 real :: maxlat, minlat
425 print*,
'- CHECK IF THE TILE CONTAINS A POLE.'
435 do j = 1, ny+1;
do i = 1, nx+1
436 if( geolat(i,j) > maxlat )
then
441 if( geolat(i,j) < minlat )
then
450 if(maxlat < 89.9 )
then
454 if(minlat > -89.9 )
then
459 print*,
"- MINLAT=", minlat,
"MAXLAT=", maxlat
460 print*,
"- NORTH POLE SUPERGRID INDEX IS ", &
461 i_north_pole, j_north_pole
462 print*,
"- SOUTH POLE SUPERGRID INDEX IS ", &
463 i_south_pole, j_south_pole
485 i_south_pole, j_south_pole, im, jm, is_north_pole, &
490 integer,
intent(in) :: im, jm
491 integer,
intent(in) :: i_north_pole, j_north_pole
492 integer,
intent(in) :: i_south_pole, j_south_pole
494 logical,
intent(out) :: is_north_pole(im,jm)
495 logical,
intent(out) :: is_south_pole(im,jm)
499 print*,
'- FIND NEAREST POLE POINTS.'
501 is_north_pole=.false.
502 is_south_pole=.false.
504 if(i_south_pole >0 .and. j_south_pole > 0)
then
505 if(mod(i_south_pole,2)==0)
then
506 do j = 1, jm;
do i = 1, im
507 if(i==i_south_pole/2 .and. (j==j_south_pole/2 &
508 .or. j==j_south_pole/2+1) )
then
509 is_south_pole(i,j) = .true.
510 print*,
"- SOUTH POLE AT I,J= ", i, j
514 do j = 1, jm;
do i = 1, im
515 if((i==i_south_pole/2 .or. i==i_south_pole/2+1) &
516 .and. (j==j_south_pole/2 .or. &
517 j==j_south_pole/2+1) )
then
518 is_south_pole(i,j) = .true.
519 print*,
"- SOUTH POLE AT I,J= ", i, j
525 if(i_north_pole >0 .and. j_north_pole > 0)
then
526 if(mod(i_north_pole,2)==0)
then
527 do j = 1, jm;
do i = 1, im
528 if(i==i_north_pole/2 .and. (j==j_north_pole/2 .or. &
529 j==j_north_pole/2+1) )
then
530 is_north_pole(i,j) = .true.
531 print*,
"- NORTH POLE AT I,J= ", i, j
535 do j = 1, jm;
do i = 1, im
536 if((i==i_north_pole/2 .or. i==i_north_pole/2+1) &
537 .and. (j==j_north_pole/2 .or. &
538 j==j_north_pole/2+1) )
then
539 is_north_pole(i,j) = .true.
540 print*,
"- NORTH POLE AT I,J= ", i, j
564 integer,
intent(in) :: im, jm
566 real,
intent(inout) :: slm(im,jm)
567 real,
intent(inout) :: oro(im,jm)
568 real,
intent(inout) :: var(im,jm)
569 real,
intent(inout) :: var4(im,jm)
570 real,
intent(inout) :: oa(im,jm,4)
571 real,
intent(inout) :: ol(im,jm,4)
573 integer :: i, j, jn, js, k
574 integer :: iw, ie, wgta, is, ise
575 integer :: in, ine, inw, isw
577 real :: slma, oroa, vara, var4a
578 real,
allocatable :: oaa(:), ola(:)
582 print*,
"- REMOVE ISOLATED POINTS."
584 allocate (oaa(4),ola(4))
586 iso_loop :
DO j=2,jm-1
593 slma=slm(iw,j)+slm(ie,j)
594 oroa=oro(iw,j)+oro(ie,j)
595 vara=var(iw,j)+var(ie,j)
596 var4a=var4(iw,j)+var4(ie,j)
598 oaa(k)=oa(iw,j,k)+oa(ie,j,k)
600 ola(k)=ol(iw,j,k)+ol(ie,j,k)
607 slma=slma+slm(inw,jn)+slm(in,jn)+slm(ine,jn)
608 oroa=oroa+oro(inw,jn)+oro(in,jn)+oro(ine,jn)
609 vara=vara+var(inw,jn)+var(in,jn)+var(ine,jn)
610 var4a=var4a+var4(inw,jn)+var4(in,jn)+var4(ine,jn)
612 oaa(k)=oaa(k)+oa(inw,jn,k)+oa(in,jn,k)+oa(ine,jn,k)
613 ola(k)=ola(k)+ol(inw,jn,k)+ol(in,jn,k)+ol(ine,jn,k)
620 slma=slma+slm(isw,js)+slm(is,js)+slm(ise,js)
621 oroa=oroa+oro(isw,js)+oro(is,js)+oro(ise,js)
622 vara=vara+var(isw,js)+var(is,js)+var(ise,js)
623 var4a=var4a+var4(isw,js)+var4(is,js)+var4(ise,js)
625 oaa(k)=oaa(k)+oa(isw,js,k)+oa(is,js,k)+oa(ise,js,k)
626 ola(k)=ola(k)+ol(isw,js,k)+ol(is,js,k)+ol(ise,js,k)
638 IF(slm(i,j).EQ.0..AND.slma.EQ.wgta)
THEN
639 print
'(" - SEA ",2F8.0," MODIFIED TO LAND",2F8.0, &
640 " AT ",2I8)',oro(i,j),var(i,j),oroa,vara,i,j
650 ELSEIF(slm(i,j).EQ.1..AND.slma.EQ.0.)
THEN
651 print
'(" - LAND",2F8.0," MODIFIED TO SEA ",2F8.0, &
652 " AT ",2I8)',oro(i,j),var(i,j),oroa,vara,i,j
693 integer,
intent(in) :: imn,jmn
694 integer,
intent(in) :: npts
695 real,
intent(in) :: lono(npts), lato(npts)
696 real,
intent(in) :: delxn
697 integer,
intent(out) :: jst,jen
698 integer,
intent(out) :: ilist(imn)
699 integer,
intent(out) :: numx
701 integer :: i2, ii, ist, ien
702 real :: minlat,maxlat,minlon,maxlon
704 minlat = minval(lato)
705 maxlat = maxval(lato)
706 minlon = minval(lono)
707 maxlon = maxval(lono)
710 jst = (minlat+90)/delxn+1
711 jen = (maxlat+90)/delxn
716 if(jen>jmn) jen = jmn
719 if((jst == 1 .OR. jen == jmn) .and. &
720 (ien-ist+1 > imn/2) )
then
725 else if( ien-ist+1 > imn/2 )
then
731 ii = lono(i2)/delxn+1
732 if(ii <0 .or. ii>imn)
then
733 print*,
"- II=",ii,imn,lono(i2),delxn
735 if( ii < imn/2 )
then
737 else if( ii > imn/2 )
then
741 if(ist<1 .OR. ist>imn)
then
742 print*, .or.
"FATAL ERROR: ist<1 ist>IMN"
745 if(ien<1 .OR. ien>imn)
then
746 print*, .or.
"FATAL ERROR: iend<1 iend>IMN"
752 ilist(i2) = ien + (i2-1)
761 ilist(i2) = ist + (i2-1)
792 glat,zavg,zslm,delxn)
798 integer,
intent(in) :: imn,jmn
799 integer,
intent(in) :: zavg(imn,jmn),zslm(imn,jmn)
800 real,
intent(in) :: lon1,lat1,lon2,lat2,delxn
801 real,
intent(in) :: glat(jmn)
803 integer :: i, j, ist, ien, jst, jen, i1
805 real :: xland,xwatr,xl1,xs1,slm,xnsum
810 if( glat(j) .gt. lat1 )
then
817 if( glat(j) .gt. lat2 )
then
825 if(ist .le.0) ist = ist + imn
826 if(ien < ist) ien = ien + imn
837 do i1 = 1, ien - ist + 1
839 if( i .le. 0) i = i + imn
840 if( i .gt. imn) i = i - imn
841 xland = xland + float(zslm(i,j))
842 xwatr = xwatr + float(1-zslm(i,j))
844 height = float(zavg(i,j))
845 if(height.lt.-990.) height = 0.0
846 xl1 = xl1 + height * float(zslm(i,j))
847 xs1 = xs1 + height * float(1-zslm(i,j))
852 slm = float(nint(xland/xnsum))
864 if( i .le. 0) i = i + imn
865 if( i .gt. imn) i = i - imn
866 height = float(zavg(i,j))
867 if(height.lt.-990.) height = 0.0
903 glat,zavg,delxn,xnsum1,xnsum2,hc)
907 integer,
intent(in) :: imn,jmn
908 integer,
intent(in) :: zavg(imn,jmn)
910 real,
intent(in) :: lon1,lat1,lon2,lat2,delxn
911 real,
intent(in) :: glat(jmn)
912 real,
intent(out) :: xnsum1,xnsum2,hc
914 integer :: i, j, ist, ien, jst, jen, i1
917 real :: xw1,xw2,xnsum
922 if( glat(j) .gt. lat1 )
then
929 if( glat(j) .gt. lat2 )
then
937 if(ist .le.0) ist = ist + imn
938 if(ien < ist) ien = ien + imn
946 do i1 = 1, ien - ist + 1
948 if( i .le. 0) i = i + imn
949 if( i .gt. imn) i = i - imn
951 height = float(zavg(i,j))
952 if(height.lt.-990.) height = 0.0
954 xw2 = xw2 + height ** 2
958 var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
959 hc = 1116.2 - 0.878 * var
965 if( i .le. 0) i = i + imn
966 if( i .gt. imn) i = i - imn
967 height = float(zavg(i,j))
968 if ( height .gt. hc ) xnsum1 = xnsum1 + 1
1005 glat,zavg,delxn,xnsum1,xnsum2,HC)
1008 integer,
intent(in) :: imn,jmn
1009 integer,
intent(in) :: zavg(imn,jmn)
1011 real,
intent(in) :: hc, glat(jmn)
1012 real,
intent(in) :: lon1,lat1,lon2,lat2,delxn
1013 real,
intent(out) :: xnsum1,xnsum2
1015 integer :: i, j, ist, ien, jst, jen, i1
1026 if( glat(j) .gt. lat1 )
then
1033 if( glat(j) .gt. lat2 )
then
1039 ist = lon1/delxn + 1
1041 if(ist .le.0) ist = ist + imn
1042 if(ien < ist) ien = ien + imn
1049 if( i .le. 0) i = i + imn
1050 if( i .gt. imn) i = i - imn
1051 height = float(zavg(i,j))
1052 if ( height .gt. hc ) xnsum1 = xnsum1 + 1
1067 character(8) :: date
1068 character(10) :: time
1069 character(5) :: zone
1070 integer,
dimension(8) :: values
1074 call date_and_time(date,time,zone,values)
1075 total=(3600*values(5)) + (60*values(6))+values(7)
1076 elapsed=float(total) + (1.0e-3*float(values(8)))
Module containing utilites used by the orog program.
real function, public get_lon_angle(dx, lat_in)
Convert the 'x' direction distance of a cubed-sphere grid point to the corresponding distance in long...
real, parameter deg2rad
degrees per radians.
subroutine, public transpose_orog(imn, jmn, glob)
Transpose the global orography data by flipping the poles and moving the starting longitude to Greenw...
subroutine, public find_poles(geolat, nx, ny, i_north_pole, j_north_pole, i_south_pole, j_south_pole)
Find the point on the model grid tile closest to the north and south pole.
real function, public timef()
Get the date/time from the system clock.
subroutine, public latlon2xyz(siz, lon, lat, x, y, z)
Convert from latitude and longitude to x,y,z coordinates.
subroutine, public get_xnsum3(lon1, lat1, lon2, lat2, imn, jmn, glat, zavg, delxn, xnsum1, xnsum2, hc)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
subroutine, public minmax(im, jm, a, title, imax, jmax)
Print out the maximum and minimum values of an array and optionally pass back the i/j location of the...
real function spherical_angle(v1, v2, v3)
Compute spherical angle.
subroutine, public remove_isolated_pts(im, jm, slm, oro, var, var4, oa, ol)
Remove isolated model points.
subroutine, public get_xnsum2(lon1, lat1, lon2, lat2, imn, jmn, glat, zavg, delxn, xnsum1, xnsum2, hc)
Count the number of high-resolution orography points that are higher than a critical value inside a m...
real, parameter earth_radius
earth radius in meters.
real function, public get_xnsum(lon1, lat1, lon2, lat2, imn, jmn, glat, zavg, zslm, delxn)
Count the number of high-resolution orography points that are higher than the model grid box average ...
real function, public get_lat_angle(dy)
Convert the 'y' direction distance of a cubed-sphere grid point to the corresponding distance in lati...
subroutine, public find_nearest_pole_points(i_north_pole, j_north_pole, i_south_pole, j_south_pole, im, jm, is_north_pole, is_south_pole)
Find the point on the model grid tile closest to the north and south pole.
logical function, public inside_a_polygon(lon1, lat1, npts, lon2, lat2)
Check if a point is inside a polygon.
real, parameter rad2deg
radians per degrees.
subroutine, public transpose_mask(imn, jmn, mask)
Transpose the global landmask by flipping the poles and moving the starting longitude to Greenwich.
subroutine, public get_index(imn, jmn, npts, lono, lato, delxn, jst, jen, ilist, numx)
Determine the location of a cubed-sphere point within the high-resolution orography data.