orog_mask_tools 1.14.0
Loading...
Searching...
No Matches
orog_utils.F90
Go to the documentation of this file.
1
4
8 module orog_utils
9
10 implicit none
11
12 private
13
14 real, parameter :: earth_radius = 6371200.
15 real, parameter :: pi=3.1415926535897931
16 real, parameter :: rad2deg = 180./3.14159265358979
17 real, parameter :: deg2rad = 3.14159265358979/180.
18
20 public :: find_poles
21 public :: get_index
22 public :: get_lat_angle
23 public :: get_lon_angle
24 public :: get_xnsum
25 public :: get_xnsum2
26 public :: get_xnsum3
27 public :: inside_a_polygon
28 public :: latlon2xyz
29 public :: minmax
30 public :: remove_isolated_pts
31 public :: timef
32 public :: transpose_orog
33 public :: transpose_mask
34
35 contains
36
49 subroutine minmax(im,jm,a,title,imax,jmax)
50
51 implicit none
52
53 character(len=8), intent(in) :: title
54
55 integer, intent(in) :: im, jm
56 integer, intent(out), optional :: imax, jmax
57
58 real, intent(in) :: a(im,jm)
59
60 integer :: i, j
61
62 real :: rmin,rmax
63
64 rmin=huge(a)
65 rmax=-rmin
66
67 if (present(imax) .and. present(jmax)) then
68 imax=0
69 jmax=0
70 endif
71
72 do j=1,jm
73 do i=1,im
74 if(a(i,j) >= rmax) then
75 rmax=a(i,j)
76 if (present(imax) .and. present(jmax)) then
77 imax = i
78 jmax = j
79 endif
80 endif
81 if(a(i,j) <= rmin)rmin=a(i,j)
82 enddo
83 enddo
84
85 write(6,150) title,rmin,rmax
86150 format(' - ',a8,' MIN=',e13.4,2x,'MAX=',e13.4)
87
88 end subroutine minmax
89
100 subroutine latlon2xyz(siz,lon, lat, x, y, z)
101
102 implicit none
103
104 integer, intent(in) :: siz
105 real, intent(in) :: lon(siz), lat(siz)
106 real, intent(out) :: x(siz), y(siz), z(siz)
107
108 integer :: n
109
110 do n = 1, siz
111 x(n) = cos(lat(n))*cos(lon(n))
112 y(n) = cos(lat(n))*sin(lon(n))
113 z(n) = sin(lat(n))
114 enddo
115
116 end subroutine latlon2xyz
117
126
127 function get_lat_angle(dy)
128
129 implicit none
130
131 real, intent(in) :: dy
132
133 real :: get_lat_angle
134
136
137 end function get_lat_angle
138
150
151 function get_lon_angle(dx,lat_in)
152
153 implicit none
154
155 real, intent(in) :: dx, lat_in
156
157 real :: get_lon_angle, lat
158
159 lat = lat_in
160 if (lat > 89.5) lat = 89.5
161 if (lat < -89.5) lat = -89.5
162
163 get_lon_angle = 2*asin( sin(dx/earth_radius*0.5)/cos(lat*deg2rad) )*rad2deg
164
165 end function get_lon_angle
166
175
176 subroutine transpose_mask(imn, jmn, mask)
177
178 implicit none
179
180 integer, intent(in) :: imn, jmn
181 integer(1), intent(inout) :: mask(imn,jmn)
182
183 integer :: i, j, it, jt
184 integer(1) :: isave
185
186! Transpose from S to N to the NCEP standard N to S.
187
188 do j=1,jmn/2
189 do i=1,imn
190 jt=jmn - j + 1
191 isave = mask(i,j)
192 mask(i,j)=mask(i,jt)
193 mask(i,jt) = isave
194 enddo
195 enddo
196
197! Data begins at dateline. NCEP standard is Greenwich.
198
199 do j=1,jmn
200 do i=1,imn/2
201 it=imn/2 + i
202 isave = mask(i,j)
203 mask(i,j)=mask(it,j)
204 mask(it,j) = isave
205 enddo
206 enddo
207
208 end subroutine transpose_mask
209
218
219 subroutine transpose_orog(imn, jmn, glob)
220
221 implicit none
222
223 integer, intent(in) :: imn, jmn
224 integer(2), intent(inout) :: glob(imn,jmn)
225
226 integer :: i, j, it, jt
227 integer(2) :: i2save
228
229! Transpose from S to N to the NCEP standard N to S.
230
231 do j=1,jmn/2
232 do i=1,imn
233 jt=jmn - j + 1
234 i2save = glob(i,j)
235 glob(i,j)=glob(i,jt)
236 glob(i,jt) = i2save
237 enddo
238 enddo
239
240! Data begins at dateline. NCEP standard is Greenwich.
241
242 do j=1,jmn
243 do i=1,imn/2
244 it=imn/2 + i
245 i2save = glob(i,j)
246 glob(i,j)=glob(it,j)
247 glob(it,j) = i2save
248 enddo
249 enddo
250
251 end subroutine transpose_orog
252
260
261 function spherical_angle(v1, v2, v3)
262
263 implicit none
264
265 real :: spherical_angle
266
267 real, parameter :: epsln30 = 1.e-30
268
269 real, intent(in) :: v1(3), v2(3), v3(3)
270
271 real :: px, py, pz, qx, qy, qz, ddd
272
273! vector product between v1 and v2
274
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)
278
279! vector product between v1 and v3
280
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);
284
285 ddd = (px*px+py*py+pz*pz)*(qx*qx+qy*qy+qz*qz);
286 if ( ddd <= 0.0 ) then
287 spherical_angle = 0.
288 else
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
293 !FIX to correctly handle co-linear points (angle near pi or 0) */
294 if (ddd < 0.) then
296 else
297 spherical_angle = 0.
298 endif
299 else
300 spherical_angle = acos( ddd )
301 endif
302 endif
303
304 end function spherical_angle
305
316
317 function inside_a_polygon(lon1, lat1, npts, lon2, lat2)
318
319 implicit none
320
321 logical inside_a_polygon
322
323 real, parameter :: epsln10 = 1.e-10
324 real, parameter :: epsln8 = 1.e-8
325 real, parameter :: range_check_criteria=0.05
326
327 integer, intent(in) :: npts
328
329 real, intent(in) :: lon1, lat1
330 real, intent(in) :: lon2(npts), lat2(npts)
331
332 integer :: i, ip1
333
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
340
341! first convert to cartesian grid.
342
343 call latlon2xyz(npts,lon2, lat2, x2, y2, z2);
344 lon1_1d(1) = lon1
345 lat1_1d(1) = lat1
346 call latlon2xyz(1,lon1_1d, lat1_1d, x1, y1, z1);
347 inside_a_polygon = .false.
348 max_x2 = maxval(x2)
349 if( x1(1) > max_x2+range_check_criteria ) return
350 min_x2 = minval(x2)
351 if( x1(1)+range_check_criteria < min_x2 ) return
352 max_y2 = maxval(y2)
353 if( y1(1) > max_y2+range_check_criteria ) return
354 min_y2 = minval(y2)
355 if( y1(1)+range_check_criteria < min_y2 ) return
356 max_z2 = maxval(z2)
357 if( z1(1) > max_z2+range_check_criteria ) return
358 min_z2 = minval(z2)
359 if( z1(1)+range_check_criteria < min_z2 ) return
360
361 pnt0(1) = x1(1)
362 pnt0(2) = y1(1)
363 pnt0(3) = z1(1)
364
365 anglesum = 0
366
367 do i = 1, npts
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 ! same as the corner point
371 inside_a_polygon = .true.
372 return
373 endif
374 ip1 = i+1
375 if(ip1>npts) ip1 = 1
376 pnt1(1) = x2(i)
377 pnt1(2) = y2(i)
378 pnt1(3) = z2(i)
379 pnt2(1) = x2(ip1)
380 pnt2(2) = y2(ip1)
381 pnt2(3) = z2(ip1)
382 angle = spherical_angle(pnt0, pnt2, pnt1);
383 anglesum = anglesum + angle
384 enddo
385
386 if(abs(anglesum-2*pi) < epsln8) then
387 inside_a_polygon = .true.
388 else
389 inside_a_polygon = .false.
390 endif
391
392 end function inside_a_polygon
393
409 subroutine find_poles(geolat, nx, ny, i_north_pole, j_north_pole, &
410 i_south_pole, j_south_pole)
411
412 implicit none
413
414 integer, intent(in) :: nx, ny
415
416 real, intent(in) :: geolat(nx+1,ny+1)
417
418 integer, intent(out) :: i_north_pole, j_north_pole
419 integer, intent(out) :: i_south_pole, j_south_pole
420
421 integer :: i, j
422
423 real :: maxlat, minlat
424
425 print*,'- CHECK IF THE TILE CONTAINS A POLE.'
426
427!--- figure out pole location.
428
429 maxlat = -90
430 minlat = 90
431 i_north_pole = 0
432 j_north_pole = 0
433 i_south_pole = 0
434 j_south_pole = 0
435 do j = 1, ny+1; do i = 1, nx+1
436 if( geolat(i,j) > maxlat ) then
437 i_north_pole=i
438 j_north_pole=j
439 maxlat = geolat(i,j)
440 endif
441 if( geolat(i,j) < minlat ) then
442 i_south_pole=i
443 j_south_pole=j
444 minlat = geolat(i,j)
445 endif
446 enddo ; enddo
447
448!--- only when maxlat is close to 90. the point is north pole
449
450 if(maxlat < 89.9 ) then
451 i_north_pole = 0
452 j_north_pole = 0
453 endif
454 if(minlat > -89.9 ) then
455 i_south_pole = 0
456 j_south_pole = 0
457 endif
458
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
464
465 end subroutine find_poles
466
483
484 subroutine find_nearest_pole_points(i_north_pole, j_north_pole, &
485 i_south_pole, j_south_pole, im, jm, is_north_pole, &
486 is_south_pole)
487
488 implicit none
489
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
493
494 logical, intent(out) :: is_north_pole(im,jm)
495 logical, intent(out) :: is_south_pole(im,jm)
496
497 integer :: i, j
498
499 print*,'- FIND NEAREST POLE POINTS.'
500
501 is_north_pole=.false.
502 is_south_pole=.false.
503
504 if(i_south_pole >0 .and. j_south_pole > 0) then
505 if(mod(i_south_pole,2)==0) then ! stretched grid
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
511 endif
512 enddo; enddo
513 else
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
520 endif
521 enddo; enddo
522 endif
523 endif
524
525 if(i_north_pole >0 .and. j_north_pole > 0) then
526 if(mod(i_north_pole,2)==0) then ! stretched grid
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
532 endif
533 enddo; enddo
534 else
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
541 endif
542 enddo; enddo
543 endif
544 endif
545
546 end subroutine find_nearest_pole_points
547
559
560 subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol)
561
562 implicit none
563
564 integer, intent(in) :: im, jm
565
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)
572
573 integer :: i, j, jn, js, k
574 integer :: iw, ie, wgta, is, ise
575 integer :: in, ine, inw, isw
576
577 real :: slma, oroa, vara, var4a
578 real, allocatable :: oaa(:), ola(:)
579
580! REMOVE ISOLATED POINTS
581
582 print*,"- REMOVE ISOLATED POINTS."
583
584 allocate (oaa(4),ola(4))
585
586 iso_loop : DO j=2,jm-1
587 jn=j-1
588 js=j+1
589 i_loop : DO i=2,im-1
590! Check the points to the 'west' and 'east'.
591 iw=i-1
592 ie=i+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)
597 DO k=1,4
598 oaa(k)=oa(iw,j,k)+oa(ie,j,k)
599! --- (*j*) fix typo:
600 ola(k)=ol(iw,j,k)+ol(ie,j,k)
601 ENDDO
602 wgta=2
603! Check the points to the 'northwest', 'north' and 'northeast'
604 in=i
605 inw=i-1
606 ine=i+1
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)
611 DO k=1,4
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)
614 ENDDO
615 wgta=wgta+3
616! Check the points to the 'southwest', 'south' and 'southeast'
617 is=i
618 isw=i-1
619 ise=i+1
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)
624 DO k=1,4
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)
627 ENDDO
628 wgta=wgta+3
629! Take the average of the surrounding the points.
630 oroa=oroa/wgta
631 vara=vara/wgta
632 var4a=var4a/wgta
633 DO k=1,4
634 oaa(k)=oaa(k)/wgta
635 ola(k)=ola(k)/wgta
636 ENDDO
637! Isolated water point.
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
641 slm(i,j)=1.
642 oro(i,j)=oroa
643 var(i,j)=vara
644 var4(i,j)=var4a
645 DO k=1,4
646 oa(i,j,k)=oaa(k)
647 ol(i,j,k)=ola(k)
648 ENDDO
649! Isolated land point.
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
653 slm(i,j)=0.
654 oro(i,j)=oroa
655 var(i,j)=vara
656 var4(i,j)=var4a
657 DO k=1,4
658 oa(i,j,k)=oaa(k)
659 ol(i,j,k)=ola(k)
660 ENDDO
661 ENDIF
662 ENDDO i_loop
663 ENDDO iso_loop
664
665 deallocate (oaa,ola)
666
667 end subroutine remove_isolated_pts
668
689 subroutine get_index(imn,jmn,npts,lonO,latO,delxn, &
690 jst,jen,ilist,numx)
691
692 implicit none
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
700
701 integer :: i2, ii, ist, ien
702 real :: minlat,maxlat,minlon,maxlon
703
704 minlat = minval(lato)
705 maxlat = maxval(lato)
706 minlon = minval(lono)
707 maxlon = maxval(lono)
708 ist = minlon/delxn+1
709 ien = maxlon/delxn+1
710 jst = (minlat+90)/delxn+1
711 jen = (maxlat+90)/delxn
712!--- add a few points to both ends of j-direction
713 jst = jst - 5
714 if(jst<1) jst = 1
715 jen = jen + 5
716 if(jen>jmn) jen = jmn
717
718!--- when around the pole, just search through all the points.
719 if((jst == 1 .OR. jen == jmn) .and. &
720 (ien-ist+1 > imn/2) )then
721 numx = imn
722 do i2 = 1, imn
723 ilist(i2) = i2
724 enddo
725 else if( ien-ist+1 > imn/2 ) then ! cross longitude = 0
726!--- find the minimum that greater than IMN/2
727!--- and maximum that less than IMN/2
728 ist = 0
729 ien = imn+1
730 do i2 = 1, npts
731 ii = lono(i2)/delxn+1
732 if(ii <0 .or. ii>imn) then
733 print*,"- II=",ii,imn,lono(i2),delxn
734 endif
735 if( ii < imn/2 ) then
736 ist = max(ist,ii)
737 else if( ii > imn/2 ) then
738 ien = min(ien,ii)
739 endif
740 enddo
741 if(ist<1 .OR. ist>imn) then
742 print*, .or."FATAL ERROR: ist<1 ist>IMN"
743 call abort()
744 endif
745 if(ien<1 .OR. ien>imn) then
746 print*, .or."FATAL ERROR: iend<1 iend>IMN"
747 call abort()
748 endif
749
750 numx = imn - ien + 1
751 do i2 = 1, numx
752 ilist(i2) = ien + (i2-1)
753 enddo
754 do i2 = 1, ist
755 ilist(numx+i2) = i2
756 enddo
757 numx = numx+ist
758 else
759 numx = ien-ist+1
760 do i2 = 1, numx
761 ilist(i2) = ist + (i2-1)
762 enddo
763 endif
764
765 end subroutine get_index
766
790
791 function get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, &
792 glat,zavg,zslm,delxn)
793
794 implicit none
795
796 real :: get_xnsum
797
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)
802
803 integer :: i, j, ist, ien, jst, jen, i1
804 real :: oro, height
805 real :: xland,xwatr,xl1,xs1,slm,xnsum
806
807!-- Figure out ist,ien,jst,jen
808
809 do j = 1, jmn
810 if( glat(j) .gt. lat1 ) then
811 jst = j
812 exit
813 endif
814 enddo
815
816 do j = 1, jmn
817 if( glat(j) .gt. lat2 ) then
818 jen = j
819 exit
820 endif
821 enddo
822
823 ist = lon1/delxn + 1
824 ien = lon2/delxn
825 if(ist .le.0) ist = ist + imn
826 if(ien < ist) ien = ien + imn
827
828!--- Compute average oro
829
830 oro = 0.0
831 xnsum = 0
832 xland = 0
833 xwatr = 0
834 xl1 = 0
835 xs1 = 0
836 do j = jst,jen
837 do i1 = 1, ien - ist + 1
838 i = ist + i1 -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))
843 xnsum = xnsum + 1.
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))
848 enddo
849 enddo
850
851 if( xnsum > 1.) then
852 slm = float(nint(xland/xnsum))
853 if(slm.ne.0.) then
854 oro= xl1 / xland
855 else
856 oro = xs1 / xwatr
857 endif
858 endif
859
860 get_xnsum = 0
861 do j = jst, jen
862 do i1= 1, ien-ist+1
863 i = ist + i1 -1
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
868 if ( height .gt. oro ) get_xnsum = get_xnsum + 1
869 enddo
870 enddo
871
872 end function get_xnsum
873
901
902 subroutine get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn, &
903 glat,zavg,delxn,xnsum1,xnsum2,hc)
904
905 implicit none
906
907 integer, intent(in) :: imn,jmn
908 integer, intent(in) :: zavg(imn,jmn)
909
910 real, intent(in) :: lon1,lat1,lon2,lat2,delxn
911 real, intent(in) :: glat(jmn)
912 real, intent(out) :: xnsum1,xnsum2,hc
913
914 integer :: i, j, ist, ien, jst, jen, i1
915
916 real :: height, var
917 real :: xw1,xw2,xnsum
918
919!-- Figure out ist,ien,jst,jen
920
921 do j = 1, jmn
922 if( glat(j) .gt. lat1 ) then
923 jst = j
924 exit
925 endif
926 enddo
927
928 do j = 1, jmn
929 if( glat(j) .gt. lat2 ) then
930 jen = j
931 exit
932 endif
933 enddo
934
935 ist = lon1/delxn + 1
936 ien = lon2/delxn
937 if(ist .le.0) ist = ist + imn
938 if(ien < ist) ien = ien + imn
939
940!--- Compute average oro
941
942 xnsum = 0
943 xw1 = 0
944 xw2 = 0
945 do j = jst,jen
946 do i1 = 1, ien - ist + 1
947 i = ist + i1 -1
948 if( i .le. 0) i = i + imn
949 if( i .gt. imn) i = i - imn
950 xnsum = xnsum + 1.
951 height = float(zavg(i,j))
952 if(height.lt.-990.) height = 0.0
953 xw1 = xw1 + height
954 xw2 = xw2 + height ** 2
955 enddo
956 enddo
957
958 var = sqrt(max(xw2/xnsum-(xw1/xnsum)**2,0.))
959 hc = 1116.2 - 0.878 * var
960 xnsum1 = 0
961 xnsum2 = 0
962 do j = jst, jen
963 do i1= 1, ien-ist+1
964 i = ist + i1 -1
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
969 xnsum2 = xnsum2 + 1
970 enddo
971 enddo
972
973 end subroutine get_xnsum2
974
1003
1004 subroutine get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn, &
1005 glat,zavg,delxn,xnsum1,xnsum2,HC)
1006 implicit none
1007
1008 integer, intent(in) :: imn,jmn
1009 integer, intent(in) :: zavg(imn,jmn)
1010
1011 real, intent(in) :: hc, glat(jmn)
1012 real, intent(in) :: lon1,lat1,lon2,lat2,delxn
1013 real, intent(out) :: xnsum1,xnsum2
1014
1015 integer :: i, j, ist, ien, jst, jen, i1
1016
1017 real :: height
1018
1019!-- Figure out ist,ien,jst,jen
1020
1021! if lat1 or lat 2 is 90 degree. set jst = JMN
1022
1023 jst = jmn
1024 jen = jmn
1025 do j = 1, jmn
1026 if( glat(j) .gt. lat1 ) then
1027 jst = j
1028 exit
1029 endif
1030 enddo
1031
1032 do j = 1, jmn
1033 if( glat(j) .gt. lat2 ) then
1034 jen = j
1035 exit
1036 endif
1037 enddo
1038
1039 ist = lon1/delxn + 1
1040 ien = lon2/delxn
1041 if(ist .le.0) ist = ist + imn
1042 if(ien < ist) ien = ien + imn
1043
1044 xnsum1 = 0
1045 xnsum2 = 0
1046 do j = jst, jen
1047 do i1= 1, ien-ist+1
1048 i = ist + i1 -1
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
1053 xnsum2 = xnsum2 + 1
1054 enddo
1055 enddo
1056
1057 end subroutine get_xnsum3
1062
1063 real function timef()
1064
1065 implicit none
1066
1067 character(8) :: date
1068 character(10) :: time
1069 character(5) :: zone
1070 integer,dimension(8) :: values
1071 integer :: total
1072 real :: elapsed
1073
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)))
1077 timef=elapsed
1078
1079 end function timef
1080
1081 end module orog_utils
Module containing utilites used by the orog program.
Definition orog_utils.F90:8
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.
real, parameter pi
pi.
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.