327 TYPE(class_gsu),
POINTER :: ptr => null()
346 REAL(4),
POINTER :: xg4(:,:), yg4(:,:)
347 REAL(8),
POINTER :: xg8(:,:), yg8(:,:)
348 TYPE(
t_nns),
POINTER :: nnp
351 REAL(8) :: xmin, ymin, xmax, ymax
352 TYPE(t_bkt),
POINTER :: b(:,:)
353 TYPE(
t_nns),
POINTER :: nnb
360 INTEGER,
POINTER :: i(:)
361 INTEGER,
POINTER :: j(:)
369 INTEGER,
POINTER :: n1(:)
370 INTEGER,
POINTER :: n2(:)
371 INTEGER,
POINTER :: di(:)
372 INTEGER,
POINTER :: dj(:)
377 REAL(8),
PARAMETER ::
pi = 3.14159265358979323846d0
378 REAL(8),
PARAMETER :: pi2 = 2d0*
pi
379 REAL(8),
PARAMETER :: pi3h = 3d0*
pi/2d0
380 REAL(8),
PARAMETER :: pio2 =
pi/2d0
381 REAL(8),
PARAMETER :: pio4 =
pi/4d0
382 REAL(8),
PARAMETER :: d2r =
pi/180d0
383 REAL(8),
PARAMETER :: r2d = 1d0/d2r
384 REAL(8),
PARAMETER :: d360 = 360d0
385 REAL(8),
PARAMETER :: d270 = 270d0
386 REAL(8),
PARAMETER :: d180 = 180d0
387 REAL(8),
PARAMETER :: d90 = 90d0
388 REAL(8),
PARAMETER :: zero = 0.0d0
389 REAL(8),
PARAMETER :: half = 0.5d0
390 REAL(8),
PARAMETER :: one = 1.0d0
391 REAL(8),
PARAMETER :: two = 2.0d0
392 REAL(8),
PARAMETER :: four = 4.0d0
394 REAL(8),
PARAMETER :: rearth = 6371229.d0
396 REAL(8),
PARAMETER :: rearth = 4.d7/pi2
398 REAL(8),
PARAMETER :: d2m = rearth*d2r
399 REAL(8),
PARAMETER :: m2d = 1d0/d2m
402 REAL(8),
PARAMETER :: eps_default = 1.0d-6
404 REAL(8),
PARAMETER :: near_pole = 10.0d0
406 INTEGER,
PARAMETER :: ncb_default = 10
408 INTEGER,
PARAMETER :: nnp_default = 2
411 INTEGER,
PARAMETER :: max_fncl_level = 3
413 INTEGER,
PARAMETER :: nfd_default = 4
418 MODULE PROCEDURE w3gsuc_ptr_r4
419 MODULE PROCEDURE w3gsuc_ptr_r8
420 MODULE PROCEDURE w3gsuc_tgt_r4
421 MODULE PROCEDURE w3gsuc_tgt_r8
424 MODULE PROCEDURE w3bbox_gsu
425 MODULE PROCEDURE w3bbox_grd_ptr_r4
426 MODULE PROCEDURE w3bbox_grd_ptr_r8
427 MODULE PROCEDURE w3bbox_grd_tgt_r4
428 MODULE PROCEDURE w3bbox_grd_tgt_r8
431 MODULE PROCEDURE w3gfcl_r4
432 MODULE PROCEDURE w3gfcl_r8
435 MODULE PROCEDURE w3gfcd_r4
436 MODULE PROCEDURE w3gfcd_r8
439 MODULE PROCEDURE w3gfpt_r4
440 MODULE PROCEDURE w3gfpt_r8
443 MODULE PROCEDURE w3gfij_r4
444 MODULE PROCEDURE w3gfij_r8
447 MODULE PROCEDURE w3grmp_r4
448 MODULE PROCEDURE w3grmp_r8
451 MODULE PROCEDURE w3grmc_r4
452 MODULE PROCEDURE w3grmc_r8
455 MODULE PROCEDURE w3cgdm_r4
456 MODULE PROCEDURE w3cgdm_r8
459 MODULE PROCEDURE w3grd0_r4
460 MODULE PROCEDURE w3grd0_r8
463 MODULE PROCEDURE w3div1_r4
464 MODULE PROCEDURE w3div1_r8
467 MODULE PROCEDURE w3div2_r4
468 MODULE PROCEDURE w3div2_r8
471 MODULE PROCEDURE w3dist_r4
472 MODULE PROCEDURE w3dist_r8
475 MODULE PROCEDURE w3splx_0d_r4
476 MODULE PROCEDURE w3splx_0d_r8
477 MODULE PROCEDURE w3splx_1d_r4
478 MODULE PROCEDURE w3splx_1d_r8
479 MODULE PROCEDURE w3splx_2d_r4
480 MODULE PROCEDURE w3splx_2d_r8
483 MODULE PROCEDURE w3spxl_0d_r4
484 MODULE PROCEDURE w3spxl_0d_r8
485 MODULE PROCEDURE w3spxl_1d_r4
486 MODULE PROCEDURE w3spxl_1d_r8
487 MODULE PROCEDURE w3spxl_2d_r4
488 MODULE PROCEDURE w3spxl_2d_r8
491 MODULE PROCEDURE w3trll_0d_r4
492 MODULE PROCEDURE w3trll_0d_r8
493 MODULE PROCEDURE w3trll_1d_r4
494 MODULE PROCEDURE w3trll_1d_r8
495 MODULE PROCEDURE w3trll_2d_r4
496 MODULE PROCEDURE w3trll_2d_r8
499 MODULE PROCEDURE w3llaz_r4
500 MODULE PROCEDURE w3llaz_r8
503 MODULE PROCEDURE w3fdwt_r4
504 MODULE PROCEDURE w3fdwt_r8
507 MODULE PROCEDURE w3ckcl_r4
508 MODULE PROCEDURE w3ckcl_r8
511 MODULE PROCEDURE w3sort_r4
512 MODULE PROCEDURE w3sort_r8
515 MODULE PROCEDURE w3isrt_r4
516 MODULE PROCEDURE w3isrt_r8
519 MODULE PROCEDURE w3inan_r4
520 MODULE PROCEDURE w3inan_r8
625 FUNCTION w3gsuc_ptr_r4( IJG, LLG, ICLO, XG, YG, &
626 NCB, NNP, DEBUG )
RESULT(GSU)
629 LOGICAL,
INTENT(IN) :: ijg
630 LOGICAL,
INTENT(IN) :: llg
631 INTEGER,
INTENT(IN) :: iclo
632 REAL(4),
POINTER :: xg(:,:)
633 REAL(4),
POINTER :: yg(:,:)
634 INTEGER,
INTENT(IN),
OPTIONAL :: ncb
635 INTEGER,
INTENT(IN),
OPTIONAL :: nnp
636 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
639 INTEGER :: lb(2), ub(2)
641 INTEGER,
SAVE :: ient = 0
642 CALL strace (ient,
'W3GSUC_PTR_R4')
645 lb(1) = lbound(xg,1); lb(2) = lbound(xg,2)
646 ub(1) = ubound(xg,1); ub(2) = ubound(xg,2)
647 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg4=xg, yg4=yg, &
648 ncb=ncb, nnp=nnp, debug=debug)
650 END FUNCTION w3gsuc_ptr_r4
654 FUNCTION w3gsuc_ptr_r8( IJG, LLG, ICLO, XG, YG, &
655 NCB, NNP, DEBUG )
RESULT(GSU)
658 LOGICAL,
INTENT(IN) :: ijg
659 LOGICAL,
INTENT(IN) :: llg
660 INTEGER,
INTENT(IN) :: iclo
661 REAL(8),
POINTER :: xg(:,:)
662 REAL(8),
POINTER :: yg(:,:)
663 INTEGER,
INTENT(IN),
OPTIONAL :: ncb
664 INTEGER,
INTENT(IN),
OPTIONAL :: nnp
665 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
668 INTEGER :: lb(2), ub(2)
670 INTEGER,
SAVE :: ient = 0
671 CALL strace (ient,
'W3GSUC_PTR_R4')
674 lb(1) = lbound(xg,1); lb(2) = lbound(xg,2)
675 ub(1) = ubound(xg,1); ub(2) = ubound(xg,2)
676 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg8=xg, yg8=yg, &
677 ncb=ncb, nnp=nnp, debug=debug)
679 END FUNCTION w3gsuc_ptr_r8
683 FUNCTION w3gsuc_tgt_r4( IJG, LLG, ICLO, LB, UB, XG, YG, &
684 NCB, NNP, DEBUG )
RESULT(GSU)
687 LOGICAL,
INTENT(IN) :: ijg
688 LOGICAL,
INTENT(IN) :: llg
689 INTEGER,
INTENT(IN) :: iclo
690 INTEGER,
INTENT(IN) :: lb(2)
691 INTEGER,
INTENT(IN) :: ub(2)
692 REAL(4),
TARGET :: xg(lb(1):ub(1),lb(2):ub(2))
693 REAL(4),
TARGET :: yg(lb(1):ub(1),lb(2):ub(2))
694 INTEGER,
INTENT(IN),
OPTIONAL :: ncb
695 INTEGER,
INTENT(IN),
OPTIONAL :: nnp
696 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
700 INTEGER,
SAVE :: ient = 0
701 CALL strace (ient,
'W3GSUC_TGT_R4')
704 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg4=xg, yg4=yg, &
705 ncb=ncb, nnp=nnp, debug=debug)
707 END FUNCTION w3gsuc_tgt_r4
711 FUNCTION w3gsuc_tgt_r8( IJG, LLG, ICLO, LB, UB, XG, YG, &
712 NCB, NNP, DEBUG )
RESULT(GSU)
715 LOGICAL,
INTENT(IN) :: ijg
716 LOGICAL,
INTENT(IN) :: llg
717 INTEGER,
INTENT(IN) :: iclo
718 INTEGER,
INTENT(IN) :: lb(2)
719 INTEGER,
INTENT(IN) :: ub(2)
720 REAL(8),
TARGET :: xg(lb(1):ub(1),lb(2):ub(2))
721 REAL(8),
TARGET :: yg(lb(1):ub(1),lb(2):ub(2))
722 INTEGER,
INTENT(IN),
OPTIONAL :: ncb
723 INTEGER,
INTENT(IN),
OPTIONAL :: nnp
724 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
728 INTEGER,
SAVE :: ient = 0
729 CALL strace (ient,
'W3GSUC_TGT_R8')
732 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg8=xg, yg8=yg, &
733 ncb=ncb, nnp=nnp, debug=debug)
735 END FUNCTION w3gsuc_tgt_r8
795 INTEGER,
SAVE :: ient = 0
796 CALL strace (ient,
'W3GSUD')
799 IF (
ASSOCIATED(gsu%PTR) )
THEN
803 IF (
ASSOCIATED(gsu%PTR%B) )
THEN
806 IF ( gsu%PTR%B(jb,ib)%N .GT. 0 )
THEN
807 DEALLOCATE(gsu%PTR%B(jb,ib)%I)
808 NULLIFY(gsu%PTR%B(jb,ib)%I)
809 DEALLOCATE(gsu%PTR%B(jb,ib)%J)
810 NULLIFY(gsu%PTR%B(jb,ib)%J)
814 DEALLOCATE(gsu%PTR%B)
884 SUBROUTINE w3gsup( GSU, IUNIT, LFULL )
886 INTEGER,
OPTIONAL,
INTENT(IN) :: iunit
887 LOGICAL,
OPTIONAL,
INTENT(IN) :: lfull
890 INTEGER,
PARAMETER :: nbyte_ptr=4
891 INTEGER,
PARAMETER :: nbyte_int=4
892 TYPE(class_gsu),
POINTER :: ptr
893 INTEGER :: ndst, k, ib, jb, nbyte
895 INTEGER,
SAVE :: ient = 0
896 CALL strace (ient,
'W3GSUP')
902 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
903 WRITE(0,
'(/1A,1A/)')
'W3GSUP ERROR -- ', &
904 'grid search utility object not created'
908 IF (
PRESENT(iunit) )
THEN
919 nbyte = (nbyte_int+nbyte_ptr*2)*
SIZE(ptr%B)
922 nbyte = nbyte + nbyte_int*2*ptr%B(jb,ib)%N
929 WRITE(ndst,
'(//80A)') (
'-',k=1,80)
930 WRITE(ndst,
'(A)')
'Report on grid search utility object'
931 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
932 WRITE(ndst,
'(A,1L2)')
'Grid ijg:',ptr%IJG
933 WRITE(ndst,
'(A,1L2)')
'Grid llg:',ptr%LLG
934 WRITE(ndst,
'(A,1I2)')
'Grid iclo:',ptr%ICLO
935 WRITE(ndst,
'(A,1L2)')
'Grid lclo:',ptr%LCLO
936 WRITE(ndst,
'(A,1I2)')
'Grid precision:',ptr%GKIND
937 WRITE(ndst,
'(A,2I6)')
'Grid lbx,lby:',ptr%LBX,ptr%LBY
938 WRITE(ndst,
'(A,2I6)')
'Grid ubx,uby:',ptr%UBX,ptr%UBY
939 WRITE(ndst,
'(A,2I6)')
'Grid nx, ny:',ptr%NX,ptr%NY
940 IF (
PRESENT(lfull) )
THEN
942 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
943 WRITE(ndst,
'(A)')
'Nearest-neighbor point search indices'
944 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
948 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
949 WRITE(ndst,
'(A)')
'Bucket-search object'
950 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
951 WRITE(ndst,
'(A,4E24.16)')
'Spatial grid search domain: ', &
952 ptr%XMIN,ptr%YMIN,ptr%XMAX,ptr%YMAX
953 WRITE(ndst,
'(A,2I6)')
'nbx,nby:',ptr%NBX,ptr%NBY
954 WRITE(ndst,
'(A,2E24.16)')
'dxb,dyb:',ptr%DXB,ptr%DYB
955 WRITE(ndst,
'(A,1F10.1)')
'Approximate memory usage (MB):', &
957 IF (
PRESENT(lfull) )
THEN
959 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
960 WRITE(ndst,
'(A)')
'Search bucket bounds:'
961 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
962 WRITE(ndst,
'(2A4,4A24)')
'IB',
'JB',
'X1',
'Y1',
'X2',
'Y2'
965 WRITE(ndst,
'(2I4,4E24.16)') ib,jb, &
966 ptr%XMIN+(ib-1)*ptr%DXB,ptr%YMIN+(jb-1)*ptr%DYB, &
967 ptr%XMIN+(ib )*ptr%DXB,ptr%YMIN+(jb )*ptr%DYB
970 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
971 WRITE(ndst,
'(A)')
'Number of cells in each search bucket:'
972 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
974 WRITE(ndst,
'(500I4)') (ptr%B(jb,ib)%N,ib=1,ptr%NBX)
976 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
977 WRITE(ndst,
'(A)')
'Search bucket cell lists:'
978 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
979 WRITE(ndst,
'(3A4,A)')
'IB',
'JB',
'NC',
': ( IC, JC), ...'
982 WRITE(ndst,
'(3I4,A,500(A,I3,A,I3,A))') ib,jb, &
983 ptr%B(jb,ib)%N,
': ', &
984 (
'(',ptr%B(jb,ib)%I(k),
',',ptr%B(jb,ib)%J(k),
') ', &
988 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
989 WRITE(ndst,
'(A)')
'Nearest-neighbor bucket search indices'
990 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
994 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
995 WRITE(ndst,
'( 80A)') (
'-',k=1,80)
1082 SUBROUTINE w3bbox_gsu( GSU, XMIN, YMIN, XMAX, YMAX )
1083 TYPE(
t_gsu),
INTENT(IN) :: gsu
1084 REAL(8),
INTENT(OUT) :: xmin, ymin, xmax, ymax
1088 INTEGER,
SAVE :: ient = 0
1089 CALL strace (ient,
'W3BBOX_GSU')
1095 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
1096 WRITE(0,
'(/1A,1A/)')
'W3BBOX_GSU ERROR -- ', &
1097 'grid search utility object not created'
1109 END SUBROUTINE w3bbox_gsu
1113 SUBROUTINE w3bbox_grd_ptr_r4( IJG, LLG, ICLO, XG, YG, &
1114 XMIN, YMIN, XMAX, YMAX )
1115 LOGICAL,
INTENT(IN) :: ijg
1116 LOGICAL,
INTENT(IN) :: llg
1117 INTEGER,
INTENT(IN) :: iclo
1118 REAL(4),
POINTER :: xg(:,:)
1119 REAL(4),
POINTER :: yg(:,:)
1120 REAL(8),
INTENT(OUT) :: xmin, ymin, xmax, ymax
1124 INTEGER :: lb(2), ub(2)
1126 INTEGER,
SAVE :: ient = 0
1127 CALL strace (ient,
'W3BBOX_GRD_PTR_R4')
1133 lb(1) = lbound(xg,1); lb(2) = lbound(xg,2)
1134 ub(1) = ubound(xg,1); ub(2) = ubound(xg,2)
1135 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg4=xg, yg4=yg, bbox_only=.true. )
1142 END SUBROUTINE w3bbox_grd_ptr_r4
1146 SUBROUTINE w3bbox_grd_ptr_r8( IJG, LLG, ICLO, XG, YG, &
1147 XMIN, YMIN, XMAX, YMAX )
1148 LOGICAL,
INTENT(IN) :: ijg
1149 LOGICAL,
INTENT(IN) :: llg
1150 INTEGER,
INTENT(IN) :: iclo
1151 REAL(8),
POINTER :: xg(:,:)
1152 REAL(8),
POINTER :: yg(:,:)
1153 REAL(8),
INTENT(OUT) :: xmin, ymin, xmax, ymax
1157 INTEGER :: lb(2), ub(2)
1159 INTEGER,
SAVE :: ient = 0
1160 CALL strace (ient,
'W3BBOX_GRD_PTR_R8')
1166 lb(1) = lbound(xg,1); lb(2) = lbound(xg,2)
1167 ub(1) = ubound(xg,1); ub(2) = ubound(xg,2)
1168 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg8=xg, yg8=yg, bbox_only=.true. )
1175 END SUBROUTINE w3bbox_grd_ptr_r8
1179 SUBROUTINE w3bbox_grd_tgt_r4( IJG, LLG, ICLO, LB, UB, XG, YG, &
1180 XMIN, YMIN, XMAX, YMAX )
1181 LOGICAL,
INTENT(IN) :: ijg
1182 LOGICAL,
INTENT(IN) :: llg
1183 INTEGER,
INTENT(IN) :: iclo
1184 INTEGER,
INTENT(IN) :: lb(2), ub(2)
1185 REAL(4),
TARGET :: xg(lb(1):ub(1),lb(2):ub(2))
1186 REAL(4),
TARGET :: yg(lb(1):ub(1),lb(2):ub(2))
1187 REAL(8),
INTENT(OUT) :: xmin, ymin, xmax, ymax
1192 INTEGER,
SAVE :: ient = 0
1193 CALL strace (ient,
'W3BBOX_GRD_TGT_R4')
1199 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg4=xg, yg4=yg, bbox_only=.true. )
1206 END SUBROUTINE w3bbox_grd_tgt_r4
1210 SUBROUTINE w3bbox_grd_tgt_r8( IJG, LLG, ICLO, LB, UB, XG, YG, &
1211 XMIN, YMIN, XMAX, YMAX )
1212 LOGICAL,
INTENT(IN) :: ijg
1213 LOGICAL,
INTENT(IN) :: llg
1214 INTEGER,
INTENT(IN) :: iclo
1215 INTEGER,
INTENT(IN) :: lb(2), ub(2)
1216 REAL(8),
TARGET :: xg(lb(1):ub(1),lb(2):ub(2))
1217 REAL(8),
TARGET :: yg(lb(1):ub(1),lb(2):ub(2))
1218 REAL(8),
INTENT(OUT) :: xmin, ymin, xmax, ymax
1223 INTEGER,
SAVE :: ient = 0
1224 CALL strace (ient,
'W3BBOX_GRD_TGT_R8')
1230 gsu = gsu_create( ijg, llg, iclo, lb, ub, xg8=xg, yg8=yg, bbox_only=.true. )
1237 END SUBROUTINE w3bbox_grd_tgt_r8
1328 FUNCTION w3gfcl_r4( GSU, XT, YT, IS, JS, XS, YS, &
1329 POLE, EPS, FNCL, DEBUG )
RESULT(INGRID)
1332 TYPE(
t_gsu),
INTENT(IN) :: gsu
1333 REAL(4),
INTENT(INOUT) :: xt
1334 REAL(4),
INTENT(INOUT) :: yt
1335 INTEGER,
INTENT(INOUT) :: is(4), js(4)
1336 REAL(4),
INTENT(INOUT) :: xs(4), ys(4)
1337 LOGICAL,
INTENT(OUT),
OPTIONAL :: pole
1338 REAL(4),
INTENT(IN),
OPTIONAL :: eps
1339 LOGICAL,
INTENT(IN),
OPTIONAL :: fncl
1340 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
1343 REAL(8) :: xt8, yt8, xs8(4), ys8(4), eps8
1345 INTEGER,
SAVE :: ient = 0
1346 CALL strace (ient,
'W3GFCL_R4')
1351 IF (
PRESENT(eps) )
THEN
1358 ingrid = w3gfcl( gsu, xt8, yt8, is, js, xs8, ys8, pole=pole, &
1359 eps=eps8, fncl=fncl, debug=debug )
1365 END FUNCTION w3gfcl_r4
1369 FUNCTION w3gfcl_r8( GSU, XT, YT, IS, JS, XS, YS, &
1370 POLE, EPS, FNCL, DEBUG )
RESULT(INGRID)
1373 TYPE(
t_gsu),
INTENT(IN) :: gsu
1374 REAL(8),
INTENT(INOUT) :: xt
1375 REAL(8),
INTENT(INOUT) :: yt
1376 INTEGER,
INTENT(INOUT) :: is(4), js(4)
1377 REAL(8),
INTENT(INOUT) :: xs(4), ys(4)
1378 LOGICAL,
INTENT(OUT),
OPTIONAL :: pole
1379 REAL(8),
INTENT(IN),
OPTIONAL :: eps
1380 LOGICAL,
INTENT(IN),
OPTIONAL :: fncl
1381 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
1385 LOGICAL :: ldbg, lplc, lfncl, incell
1386 INTEGER :: i, j, k, l, n, ib, jb
1387 LOGICAL :: ijg, llg, lclo, l360
1388 INTEGER :: iclo, gkind
1389 INTEGER :: lbx, lby, ubx, uby, nx, ny
1390 REAL(4),
POINTER :: xg4(:,:), yg4(:,:)
1391 REAL(8),
POINTER :: xg8(:,:), yg8(:,:)
1393 REAL(8) :: dxb, dyb, xmin, xmax, ymin, ymax
1394 TYPE(t_bkt),
POINTER :: b(:,:)
1395 TYPE(
t_nns),
POINTER :: nnb
1397 INTEGER :: nlevel, lvl, lvl1, n1, ib0, jb0, ib1, jb1, k1
1398 INTEGER :: is1(4), js1(4)
1399 REAL(8) :: xs1(4), ys1(4), xsm, ysm, dd, dd1
1401 INTEGER,
SAVE :: ient = 0
1402 CALL strace (ient,
'W3GFCL_R8')
1408 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
1409 WRITE(0,
'(/2A/)')
'W3GFCL_R8 ERROR -- ', &
1410 'grid search utility object not created'
1413 IF (
PRESENT(eps) )
THEN
1414 IF ( eps .LT. zero )
THEN
1415 WRITE(0,
'(/2A/)')
'W3GFCL_R8 ERROR -- ', &
1416 'EPS parameter must be >= 0'
1427 IF (
PRESENT(fncl) )
THEN
1432 IF (
PRESENT(debug) )
THEN
1444 gkind = gsu%PTR%GKIND
1445 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
1446 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
1447 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
1448 IF ( gkind.EQ.4 )
THEN
1449 xg4 => gsu%PTR%XG4; yg4 => gsu%PTR%YG4;
1451 xg8 => gsu%PTR%XG8; yg8 => gsu%PTR%YG8;
1453 nbx = gsu%PTR%NBX; nby = gsu%PTR%NBY;
1454 dxb = gsu%PTR%DXB; dyb = gsu%PTR%DYB;
1455 xmin = gsu%PTR%XMIN; ymin = gsu%PTR%YMIN;
1456 xmax = gsu%PTR%XMAX; ymax = gsu%PTR%YMAX;
1464 xt = mod(xt,real(d360,8))
1465 IF ( lclo .OR. l360 )
THEN
1466 IF ( xt.LT.zero ) xt = xt + d360
1468 IF ( xt.GT.d180 ) xt = xt - d360
1471 IF ( ldbg )
WRITE(*,
'(/A,2E24.16)')
'W3GFCL_R8 - TARGET POINT:',xt,yt
1474 IF ( .NOT.lfncl )
THEN
1475 IF ( xt.LT.xmin-leps .OR. xt.GT.xmax+leps .OR. &
1476 yt.LT.ymin-leps .OR. yt.GT.ymax+leps )
THEN
1477 IF ( ldbg )
WRITE(*,
'(A)') &
1478 'W3GFCL_R8 - TARGET POINT OUTSIDE SEARCH DOMAIN'
1484 ib = max(int((xt-xmin)/dxb)+1,1);
IF ( .NOT.lclo ) ib = min(ib,nbx);
1485 jb = max(int((yt-ymin)/dyb)+1,1); jb = min(jb,nby);
1491 WRITE(*,
'(A,3I6,4E24.16)') &
1492 'W3GFCL_R8 - BUCKET SEARCH:',ib,jb,b(jb,ib)%N, &
1493 xmin+(ib-1)*dxb,ymin+(jb-1)*dyb,xmin+ib*dxb,ymin+jb*dyb
1494 cell_loop:
DO k=1,b(jb,ib)%N
1496 is(1) = b(jb,ib)%I(k) ; js(1) = b(jb,ib)%J(k) ;
1497 is(2) = b(jb,ib)%I(k)+1; js(2) = b(jb,ib)%J(k) ;
1498 is(3) = b(jb,ib)%I(k)+1; js(3) = b(jb,ib)%J(k)+1;
1499 is(4) = b(jb,ib)%I(k) ; js(4) = b(jb,ib)%J(k)+1;
1503 IF ( mod(iclo,2).EQ.0 ) &
1504 is(l) = lbx + mod(nx - 1 + mod(is(l) - lbx + 1, nx), nx)
1505 IF ( mod(iclo,3).EQ.0 ) &
1506 js(l) = lby + mod(ny - 1 + mod(js(l) - lby + 1, ny), ny)
1507 IF ( iclo.EQ.
iclo_trpl .AND. js(l).GT.uby )
THEN
1508 is(l) = ubx + lbx - is(l)
1509 js(l) = 2*uby - js(l) + 1
1513 IF ( gkind.EQ.4 )
THEN
1514 xs(l) = xg4(is(l),js(l)); ys(l) = yg4(is(l),js(l));
1516 xs(l) = xg8(is(l),js(l)); ys(l) = yg8(is(l),js(l));
1519 IF ( gkind.EQ.4 )
THEN
1520 xs(l) = xg4(js(l),is(l)); ys(l) = yg4(js(l),is(l));
1522 xs(l) = xg8(js(l),is(l)); ys(l) = yg8(js(l),is(l));
1527 xs(l) = mod(xs(l),real(d360,8))
1528 IF ( lclo .OR. l360 )
THEN
1529 IF ( xs(l).LT.zero ) xs(l) = xs(l) + d360
1531 IF ( xs(l).GT.d180 ) xs(l) = xs(l) - d360
1536 WRITE(*,
'(A,3I6,4(/A,1I1,A,2I6,2E24.16))') &
1537 'W3GFCL_R8 - CHECK CELL:',ib,jb,k, &
1538 (
' CORNER(',l,
'):',is(l),js(l),xs(l),ys(l),l=1,4)
1540 incell = w3ckcl(llg,xt,yt,4,xs,ys,lplc,leps,ldbg)
1541 IF ( ldbg )
WRITE(*,
'(A,1L2)')
'W3GFCL_R8 - INCELL:',incell
1545 WRITE(*,
'(A,3I6,4(2I6))') &
1546 'W3GFCL_R8 - ENCLOSING CELL:',ib,jb,k,(is(l),js(l),l=1,4)
1547 IF (
PRESENT(pole) ) pole = lplc
1552 IF ( ingrid )
RETURN
1553 IF ( .NOT.lfncl )
RETURN
1563 nnb =
w3nnsc(nint(half*max(nbx,nby)))
1564 IF ( ldbg )
WRITE(*,
'(A,3I6)') &
1565 'W3GFCL_R8 - CLOSEST CELL SEARCH:',ib0,jb0,nnb%NLVL
1566 level_loop:
DO lvl=0,nnb%NLVL
1568 nnbr_loop:
DO n=nnb%N1(lvl),nnb%N2(lvl)
1569 ib = ib0 + nnb%DI(n); jb = jb0 + nnb%DJ(n);
1570 IF ( ib.LT.1 .OR. ib.GT.nbx ) cycle nnbr_loop
1571 IF ( jb.LT.1 .OR. jb.GT.nby ) cycle nnbr_loop
1572 IF ( ldbg )
WRITE(*,
'(A,5I6)') &
1573 'W3GFCL_R8 - CHECK BUCKET:',lvl,n,ib,jb,b(jb,ib)%N
1574 cell_loop2:
DO k=1,b(jb,ib)%N
1576 is(1) = b(jb,ib)%I(k) ; js(1) = b(jb,ib)%J(k) ;
1577 is(2) = b(jb,ib)%I(k)+1; js(2) = b(jb,ib)%J(k) ;
1578 is(3) = b(jb,ib)%I(k)+1; js(3) = b(jb,ib)%J(k)+1;
1579 is(4) = b(jb,ib)%I(k) ; js(4) = b(jb,ib)%J(k)+1;
1583 IF ( mod(iclo,2).EQ.0 ) &
1584 is(l) = lbx + mod(nx - 1 + mod(is(l) - lbx + 1, nx), nx)
1585 IF ( mod(iclo,3).EQ.0 ) &
1586 js(l) = lby + mod(ny - 1 + mod(js(l) - lby + 1, ny), ny)
1587 IF ( iclo.EQ.
iclo_trpl .AND. js(l).GT.uby )
THEN
1588 is(l) = ubx + lbx - is(l)
1589 js(l) = 2*uby - js(l) + 1
1593 IF ( gkind.EQ.4 )
THEN
1594 xs(l) = xg4(is(l),js(l)); ys(l) = yg4(is(l),js(l));
1596 xs(l) = xg8(is(l),js(l)); ys(l) = yg8(is(l),js(l));
1599 IF ( gkind.EQ.4 )
THEN
1600 xs(l) = xg4(js(l),is(l)); ys(l) = yg4(js(l),is(l));
1602 xs(l) = xg8(js(l),is(l)); ys(l) = yg8(js(l),is(l));
1607 xs(l) = mod(xs(l),real(d360,8))
1608 IF ( lclo .OR. l360 )
THEN
1609 IF ( xs(l).LT.zero ) xs(l) = xs(l) + d360
1611 IF ( xs(l).GT.d180 ) xs(l) = xs(l) - d360
1616 xsm = sum(xs)/four; ysm = sum(ys)/four;
1617 dd = w3dist(llg,xt,yt,xsm,ysm)
1619 WRITE(*,
'(A,5I6,3E24.16,4(/A,1I1,A,2I6,2E24.16))') &
1620 'W3GFCL_R8 - CHECK CELL:',lvl,n,ib,jb,k,xsm,ysm,dd, &
1621 (
' CORNER(',l,
'):',is(l),js(l),xs(l),ys(l),l=1,4)
1637 IF ( found ) nlevel = nlevel + 1
1638 IF ( nlevel .GE. max_fncl_level )
EXIT level_loop
1647 WRITE(*,
'(A,5I6,1E24.16,4(/A,1I1,A,2I6,2E24.16))') &
1648 'W3GFCL_R8 - CLOSEST CELL:',lvl1,n1,ib1,jb1,k1,dd1, &
1649 (
' CORNER(',l,
'):',is(l),js(l),xs(l),ys(l),l=1,4)
1657 IF ( abs(xs(j)-xs(i)) .GT. d180 ) n = n + 1
1661 lplc = n.EQ.1 .OR. count(abs(ys).EQ.d90).EQ.1
1662 IF ( lplc .AND. ldbg ) &
1663 WRITE(*,
'(A)')
'W3GFCL_R8 - CELL INCLUDES A POLE'
1667 IF (
PRESENT(pole) ) pole = lplc
1669 END FUNCTION w3gfcl_r8
1755 FUNCTION w3gfcd_r4( GSU, XT, YT, IS, JS, XS, YS, &
1756 POLE, EPS, DEBUG )
RESULT(INGRID)
1759 TYPE(
t_gsu),
INTENT(IN) :: gsu
1760 REAL(4),
INTENT(INOUT) :: xt
1761 REAL(4),
INTENT(INOUT) :: yt
1762 INTEGER,
INTENT(INOUT) :: is(4), js(4)
1763 REAL(4),
INTENT(INOUT) :: xs(4), ys(4)
1764 LOGICAL,
INTENT(OUT),
OPTIONAL :: pole
1765 REAL(4),
INTENT(IN),
OPTIONAL :: eps
1766 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
1769 REAL(8) :: xt8, yt8, xs8(4), ys8(4), eps8
1771 INTEGER,
SAVE :: ient = 0
1772 CALL strace (ient,
'W3GFCD_R4')
1777 IF (
PRESENT(eps) )
THEN
1784 ingrid = w3gfcd( gsu, xt8, yt8, is, js, xs8, ys8, pole=pole, &
1785 eps=eps8, debug=debug )
1791 END FUNCTION w3gfcd_r4
1795 FUNCTION w3gfcd_r8( GSU, XT, YT, IS, JS, XS, YS, &
1796 POLE, EPS, DEBUG )
RESULT(INGRID)
1799 TYPE(
t_gsu),
INTENT(IN) :: gsu
1800 REAL(8),
INTENT(INOUT) :: xt
1801 REAL(8),
INTENT(INOUT) :: yt
1802 INTEGER,
INTENT(INOUT) :: is(4), js(4)
1803 REAL(8),
INTENT(INOUT) :: xs(4), ys(4)
1804 LOGICAL,
INTENT(OUT),
OPTIONAL :: pole
1805 REAL(8),
INTENT(IN),
OPTIONAL :: eps
1806 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
1810 LOGICAL :: ldbg, lplc
1812 LOGICAL :: ijg, llg, lclo, l360
1813 INTEGER :: iclo, gkind
1814 INTEGER :: lbx, lby, ubx, uby, nx, ny
1815 INTEGER :: lxc, lyc, uxc, uyc
1816 REAL(4),
POINTER :: xg4(:,:), yg4(:,:)
1817 REAL(8),
POINTER :: xg8(:,:), yg8(:,:)
1819 INTEGER,
SAVE :: ient = 0
1820 CALL strace (ient,
'W3GFCD_R8')
1826 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
1827 WRITE(0,
'(/2A/)')
'W3GFCD_R8 ERROR -- ', &
1828 'grid search utility object not created'
1831 IF (
PRESENT(eps) )
THEN
1832 IF ( eps .LT. zero )
THEN
1833 WRITE(0,
'(/2A/)')
'W3GFCD_R8 ERROR -- ', &
1834 'EPS parameter must be >= 0'
1845 IF (
PRESENT(debug) )
THEN
1857 gkind = gsu%PTR%GKIND
1858 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
1859 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
1860 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
1861 IF ( gkind.EQ.4 )
THEN
1862 xg4 => gsu%PTR%XG4; yg4 => gsu%PTR%YG4;
1864 xg8 => gsu%PTR%XG8; yg8 => gsu%PTR%YG8;
1871 xt = mod(xt,real(d360,8))
1872 IF ( lclo .OR. l360 )
THEN
1873 IF ( xt.LT.zero ) xt = xt + d360
1875 IF ( xt.GT.d180 ) xt = xt - d360
1879 WRITE(*,
'(/A,2E24.16)')
'W3GFCD_R8 - TARGET POINT:',xt,yt
1882 lxc = lbx; lyc = lby;
1883 SELECT CASE ( iclo )
1885 uxc = ubx-1; uyc = uby-1;
1887 uxc = ubx; uyc = uby-1;
1889 uxc = ubx-1; uyc = uby;
1891 uxc = ubx; uyc = uby;
1893 uxc = ubx; uyc = uby;
1899 cell_loop:
DO i=lxc,uxc
1902 is(1) = i ; js(1) = j ;
1903 is(2) = i+1; js(2) = j ;
1904 is(3) = i+1; js(3) = j+1;
1905 is(4) = i ; js(4) = j+1;
1909 IF ( mod(iclo,2).EQ.0 ) &
1910 is(l) = lbx + mod(nx - 1 + mod(is(l) - lbx + 1, nx), nx)
1911 IF ( mod(iclo,3).EQ.0 ) &
1912 js(l) = lby + mod(ny - 1 + mod(js(l) - lby + 1, ny), ny)
1913 IF ( iclo.EQ.
iclo_trpl .AND. js(l).GT.uby )
THEN
1914 is(l) = ubx + lbx - is(l)
1915 js(l) = 2*uby - js(l) + 1
1919 IF ( gkind.EQ.4 )
THEN
1920 xs(l) = xg4(is(l),js(l)); ys(l) = yg4(is(l),js(l));
1922 xs(l) = xg8(is(l),js(l)); ys(l) = yg8(is(l),js(l));
1925 IF ( gkind.EQ.4 )
THEN
1926 xs(l) = xg4(js(l),is(l)); ys(l) = yg4(js(l),is(l));
1928 xs(l) = xg8(js(l),is(l)); ys(l) = yg8(js(l),is(l));
1933 xs(l) = mod(xs(l),real(d360,8))
1934 IF ( lclo .OR. l360 )
THEN
1935 IF ( xs(l).LT.zero ) xs(l) = xs(l) + d360
1937 IF ( xs(l).GT.d180 ) xs(l) = xs(l) - d360
1942 WRITE(*,
'(A,4(/A,1I1,A,2I6,2E24.16))') &
1943 'W3GFCD_R8 - CHECK CELL:', &
1944 (
' CORNER(',l,
'):',is(l),js(l),xs(l),ys(l),l=1,4)
1946 ingrid = w3ckcl(llg,xt,yt,4,xs,ys,lplc,leps,ldbg)
1947 IF ( ldbg )
WRITE(*,
'(A,1L2)')
'W3GFCD_R8 - INGRID:',ingrid
1951 WRITE(*,
'(A,4(2I6))') &
1952 'W3GFCD_R8 - ENCLOSING CELL:',(is(l),js(l),l=1,4)
1953 IF (
PRESENT(pole) ) pole = lplc
1959 END FUNCTION w3gfcd_r8
2038 FUNCTION w3gfpt_r4( GSU, XTIN, YTIN, IX, IY, EPS, DCIN, DEBUG ) &
2042 TYPE(
t_gsu),
INTENT(IN) :: gsu
2043 REAL(4),
INTENT(IN) :: xtin
2044 REAL(4),
INTENT(IN) :: ytin
2045 INTEGER,
INTENT(OUT) :: ix, iy
2046 REAL(4),
INTENT(IN),
OPTIONAL :: eps
2047 REAL(4),
INTENT(IN),
OPTIONAL :: dcin
2048 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
2051 REAL(8) :: xt8, yt8, eps8, dcin8
2053 INTEGER,
SAVE :: ient = 0
2054 CALL strace (ient,
'W3GFPT_R4')
2058 xt8 = xtin; yt8 = ytin;
2059 IF (
PRESENT(eps) )
THEN
2064 IF (
PRESENT(dcin) )
THEN
2071 ingrid = w3gfpt( gsu, xt8, yt8, ix, iy, eps=eps8, dcin=dcin8, &
2074 END FUNCTION w3gfpt_r4
2078 FUNCTION w3gfpt_r8( GSU, XTIN, YTIN, IX, IY, EPS, DCIN, DEBUG ) &
2082 TYPE(
t_gsu),
INTENT(IN) :: gsu
2083 REAL(8),
INTENT(IN) :: xtin
2084 REAL(8),
INTENT(IN) :: ytin
2085 INTEGER,
INTENT(OUT) :: ix, iy
2086 REAL(8),
INTENT(IN),
OPTIONAL :: eps
2087 REAL(8),
INTENT(IN),
OPTIONAL :: dcin
2088 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
2091 REAL(8),
PARAMETER :: big = 1d16
2092 REAL(8) :: leps, ldcin
2093 LOGICAL :: ldbg, fncl
2095 INTEGER :: is(4), js(4)
2096 REAL(8) :: xt, yt, xs(4), ys(4)
2097 REAL(8) :: xtc, ytc, xsc(4), ysc(4)
2098 REAL(8) :: ixr, jxr, dd, lon0, lat0, dmin
2101 INTEGER,
SAVE :: ient = 0
2102 CALL strace (ient,
'W3GFPT_R8')
2108 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
2109 WRITE(0,
'(/2A/)')
'W3GFPT_R8 ERROR -- ', &
2110 'grid search utility object not created'
2113 IF (
PRESENT(eps) )
THEN
2114 IF ( eps .LT. zero )
THEN
2115 WRITE(0,
'(/2A/)')
'W3GFPT_R8 ERROR -- ', &
2116 'EPS parameter must be >= 0'
2123 IF (
PRESENT(dcin) )
THEN
2124 IF ( dcin .LT. zero )
THEN
2125 WRITE(0,
'(/2A/)')
'W3GFPT_R8 ERROR -- ', &
2126 'DCIN parameter must be >= 0'
2137 IF (
PRESENT(debug) )
THEN
2151 xt = xtin; yt = ytin;
2153 WRITE(*,
'(/A,2E24.16)')
'W3GFPT_R8 - TARGET POINT:',xt,yt
2158 fncl = ldcin .GT. zero
2159 ingrid = w3gfcl( gsu, xt, yt, is, js, xs, ys, eps=leps, fncl=fncl, debug=ldbg )
2160 IF ( .NOT.ingrid .AND. .NOT.fncl )
RETURN
2163 IF ( .NOT.ingrid .AND. fncl )
THEN
2165 lon0 = sum(xs)/four; lat0 = sum(ys)/four;
2166 IF ( d90-abs(lat0).GT.near_pole )
THEN
2168 CALL getpqr(xt,yt,xs,ys,ixr,jxr,eps=leps,debug=ldbg)
2171 CALL w3splx(lon0,lat0,zero,xt,yt,xtc,ytc)
2173 CALL w3splx(lon0,lat0,zero,xs(i),ys(i),xsc(i),ysc(i))
2175 CALL getpqr(xtc,ytc,xsc,ysc,ixr,jxr,eps=leps,debug=ldbg)
2178 ingrid = abs(ixr-half).LE.dd .AND. abs(jxr-half).LE.dd
2185 dd = w3dist( llg, xt, yt, xs(l), ys(l) )
2186 IF ( dd .LT. dmin )
THEN
2187 dmin = dd; ix = is(l); iy = js(l);
2192 END FUNCTION w3gfpt_r8
2276 FUNCTION w3gfij_r4( GSU, XTIN, YTIN, IX, JX, EPS, DCIN, DEBUG ) &
2280 TYPE(
t_gsu),
INTENT(IN) :: gsu
2281 REAL(4),
INTENT(IN) :: xtin
2282 REAL(4),
INTENT(IN) :: ytin
2283 REAL(4),
INTENT(OUT) :: ix
2284 REAL(4),
INTENT(OUT) :: jx
2285 REAL(4),
INTENT(IN),
OPTIONAL :: eps
2286 REAL(4),
INTENT(IN),
OPTIONAL :: dcin
2287 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
2290 REAL(8) :: xt8, yt8, ix8, jx8, eps8, dcin8
2292 INTEGER,
SAVE :: ient = 0
2293 CALL strace (ient,
'W3GFIJ_R4')
2297 xt8 = xtin; yt8 = ytin;
2298 IF (
PRESENT(eps) )
THEN
2303 IF (
PRESENT(dcin) )
THEN
2310 ingrid = w3gfij( gsu, xt8, yt8, ix8, jx8, eps=eps8, dcin=dcin8, &
2316 END FUNCTION w3gfij_r4
2320 FUNCTION w3gfij_r8( GSU, XTIN, YTIN, IX, JX, EPS, DCIN, DEBUG ) &
2324 TYPE(
t_gsu),
INTENT(IN) :: gsu
2325 REAL(8),
INTENT(IN) :: xtin
2326 REAL(8),
INTENT(IN) :: ytin
2327 REAL(8),
INTENT(OUT) :: ix
2328 REAL(8),
INTENT(OUT) :: jx
2329 REAL(8),
INTENT(IN),
OPTIONAL :: eps
2330 REAL(8),
INTENT(IN),
OPTIONAL :: dcin
2331 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
2334 REAL(8) :: leps, ldcin
2336 LOGICAL :: ldbg, fncl, pole
2337 INTEGER :: is(4), js(4)
2338 REAL(8) :: xt, yt, xs(4), ys(4)
2339 REAL(8) :: xtc, ytc, xsc(4), ysc(4)
2340 REAL(8) :: ixr, jxr, dd, lon0, lat0
2342 INTEGER,
SAVE :: ient = 0
2343 CALL strace (ient,
'W3GFIJ_R8')
2349 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
2350 WRITE(0,
'(/2A/)')
'W3GFIJ_R8 ERROR -- ', &
2351 'grid search utility object not created'
2354 IF (
PRESENT(eps) )
THEN
2355 IF ( eps .LT. zero )
THEN
2356 WRITE(0,
'(/2A/)')
'W3GFIJ_R8 ERROR -- ', &
2357 'EPS parameter must be >= 0'
2364 IF (
PRESENT(dcin) )
THEN
2365 IF ( dcin .LT. zero )
THEN
2366 WRITE(0,
'(/2A/)')
'W3GFIJ_R8 ERROR -- ', &
2367 'DCIN parameter must be >= 0'
2378 IF (
PRESENT(debug) )
THEN
2384 xt = xtin; yt = ytin;
2385 IF ( ldbg )
WRITE(*,
'(/A,2E24.16)')
'W3GFIJ_R8 - TARGET POINT:',xt,yt
2390 fncl = ldcin .GT. zero
2391 ingrid = w3gfcl(gsu,xt,yt,is,js,xs,ys,pole=pole,eps=leps,fncl=fncl,debug=ldbg)
2392 IF ( .NOT.ingrid .AND. .NOT.fncl )
RETURN
2395 lon0 = sum(xs)/four; lat0 = sum(ys)/four;
2396 IF ( d90-abs(lat0).GT.near_pole )
THEN
2398 CALL getpqr(xt,yt,xs,ys,ixr,jxr,eps=leps,debug=ldbg)
2401 CALL w3splx(lon0,lat0,zero,xt,yt,xtc,ytc)
2403 CALL w3splx(lon0,lat0,zero,xs(i),ys(i),xsc(i),ysc(i))
2405 CALL getpqr(xtc,ytc,xsc,ysc,ixr,jxr,eps=leps,debug=ldbg)
2408 WRITE(*,
'(A,2L2,2E24.16)')
'W3GFIJ_R8 - RELATIVE:',ingrid,fncl,ixr,jxr
2411 IF ( .NOT.ingrid .AND. fncl )
THEN
2413 ingrid = abs(ixr-half).LE.dd .AND. abs(jxr-half).LE.dd
2417 ix = is(1)+ixr; jx = js(1)+jxr;
2419 WRITE(*,
'(A,2L2,2E24.16)')
'W3GFIJ_R8 - ABSOLUTE:',ingrid,fncl,ix,jx
2421 END FUNCTION w3gfij_r8
2521 FUNCTION w3grmp_r4( GSU, XTIN, YTIN, IS, JS, RW, EPS, &
2522 DCIN, MASK, MSKC, NNBR, DEBUG )
RESULT(INGRID)
2525 TYPE(
t_gsu),
INTENT(IN) :: gsu
2526 REAL(4),
INTENT(IN) :: xtin
2527 REAL(4),
INTENT(IN) :: ytin
2528 INTEGER,
INTENT(OUT) :: is(4)
2529 INTEGER,
INTENT(OUT) :: js(4)
2530 REAL(4),
INTENT(OUT) :: rw(4)
2531 REAL(4),
INTENT(IN) ,
OPTIONAL :: eps
2532 REAL(4),
INTENT(IN) ,
OPTIONAL :: dcin
2533 LOGICAL,
INTENT(IN) ,
OPTIONAL :: mask(:,:)
2534 INTEGER,
INTENT(OUT) ,
OPTIONAL :: mskc
2535 INTEGER,
INTENT(INOUT),
OPTIONAL :: nnbr
2536 LOGICAL,
INTENT(IN) ,
OPTIONAL :: debug
2539 REAL(8) :: xt8, yt8, rw8(4), eps8, dcin8
2541 INTEGER,
SAVE :: ient = 0
2542 CALL strace (ient,
'W3GRMP_R4')
2546 xt8 = xtin; yt8 = ytin;
2547 IF (
PRESENT(eps) )
THEN
2552 IF (
PRESENT(dcin) )
THEN
2559 ingrid = w3grmp( gsu, xt8, yt8, is, js, rw8, &
2560 eps=eps8, dcin=dcin8, &
2561 mask=mask, mskc=mskc, nnbr=nnbr, debug=debug )
2566 END FUNCTION w3grmp_r4
2570 FUNCTION w3grmp_r8( GSU, XTIN, YTIN, IS, JS, RW, EPS, &
2571 DCIN, MASK, MSKC, NNBR, DEBUG )
RESULT(INGRID)
2574 TYPE(
t_gsu),
INTENT(IN) :: gsu
2575 REAL(8),
INTENT(IN) :: xtin
2576 REAL(8),
INTENT(IN) :: ytin
2577 INTEGER,
INTENT(OUT) :: is(4)
2578 INTEGER,
INTENT(OUT) :: js(4)
2579 REAL(8),
INTENT(OUT) :: rw(4)
2580 REAL(8),
INTENT(IN) ,
OPTIONAL :: eps
2581 REAL(8),
INTENT(IN) ,
OPTIONAL :: dcin
2582 LOGICAL,
INTENT(IN) ,
OPTIONAL :: mask(:,:)
2583 INTEGER,
INTENT(OUT) ,
OPTIONAL :: mskc
2584 INTEGER,
INTENT(INOUT),
OPTIONAL :: nnbr
2585 LOGICAL,
INTENT(IN) ,
OPTIONAL :: debug
2588 REAL(8),
PARAMETER :: big = 1d16
2589 REAL(8),
PARAMETER :: small = 1d-6
2591 LOGICAL :: ldbg, fncl, pole
2593 LOGICAL :: m, msk(4)
2594 INTEGER :: lvl, n, ns, icc, jcc
2595 REAL(8) :: xt, yt, xs(4), ys(4), dw(4)
2596 REAL(8) :: xtc, ytc, xsc(4), ysc(4)
2597 REAL(8) :: ldcin, ixr, jxr, x, y, d(4), dd, dmin, dsum, lon0, lat0
2598 LOGICAL :: ijg, llg, lclo
2599 INTEGER :: iclo, gkind
2600 INTEGER :: lbx, lby, ubx, uby, nx, ny
2601 REAL(4),
POINTER :: xg4(:,:), yg4(:,:)
2602 REAL(8),
POINTER :: xg8(:,:), yg8(:,:)
2603 TYPE(
t_nns),
POINTER :: nnp
2605 INTEGER,
SAVE :: ient = 0
2606 CALL strace (ient,
'W3GRMP_R8')
2612 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
2613 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2614 'grid search utility object not created'
2618 IF (
PRESENT(eps) )
THEN
2619 IF ( eps .LT. zero )
THEN
2620 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2621 'EPS parameter must be >= 0'
2629 IF (
PRESENT(dcin) )
THEN
2630 IF ( dcin .LT. zero )
THEN
2631 WRITE(0,
'(/2A/)')
'W3GRMP_R4 ERROR -- ', &
2632 'DCIN parameter must be >= 0'
2640 IF (
PRESENT(mask) )
THEN
2641 IF ( .NOT.
PRESENT(mskc) )
THEN
2642 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2643 'MSKC must be specified with MASK'
2646 IF (
PRESENT(nnbr) )
THEN
2647 IF ( .NOT.
ASSOCIATED(gsu%PTR%NNP) )
THEN
2648 WRITE(0,
'(/3A/)')
'W3GRMP_R8 ERROR -- ', &
2649 'MASK and NNBR input specified, ', &
2650 'but grid point-search object not created'
2653 IF ( nnbr .LE. 0 .OR. nnbr .GT. 4 )
THEN
2654 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2655 'NNBR must be >= 1 AND <= 4'
2664 IF (
PRESENT(debug) )
THEN
2675 gkind = gsu%PTR%GKIND
2676 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
2677 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
2678 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
2679 IF ( gkind.EQ.4 )
THEN
2680 xg4 => gsu%PTR%XG4; yg4 => gsu%PTR%YG4;
2682 xg8 => gsu%PTR%XG8; yg8 => gsu%PTR%YG8;
2686 IF (
PRESENT(mask) )
THEN
2688 IF ( .NOT.(ubound(mask,1).EQ.nx.AND. &
2689 ubound(mask,2).EQ.ny) )
THEN
2690 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2691 'MASK array size does not agree with GSU index bounds'
2695 IF ( .NOT.(ubound(mask,2).EQ.nx.AND. &
2696 ubound(mask,1).EQ.ny) )
THEN
2697 WRITE(0,
'(/2A/)')
'W3GRMP_R8 ERROR -- ', &
2698 'MASK array size does not agree with GSU index bounds'
2706 xt = xtin; yt = ytin;
2707 IF ( ldbg )
WRITE(*,
'(/A,2E24.16)')
'W3GRMP_R8 - TARGET POINT:',xt,yt
2712 fncl = ldcin .GT. zero
2713 ingrid = w3gfcl(gsu,xt,yt,is,js,xs,ys,pole=pole,eps=leps,fncl=fncl,debug=ldbg)
2714 IF ( .NOT.ingrid .AND. .NOT.fncl )
RETURN
2717 lon0 = sum(xs)/four; lat0 = sum(ys)/four;
2718 IF ( d90-abs(lat0).GT.near_pole )
THEN
2720 CALL getpqr(xt,yt,xs,ys,ixr,jxr,eps=leps,debug=ldbg)
2723 CALL w3splx(lon0,lat0,zero,xt,yt,xtc,ytc)
2725 CALL w3splx(lon0,lat0,zero,xs(i),ys(i),xsc(i),ysc(i))
2727 CALL getpqr(xtc,ytc,xsc,ysc,ixr,jxr,eps=leps,debug=ldbg)
2729 dw(1) = (one-ixr)*(one-jxr)
2730 dw(2) = ixr*(one-jxr)
2732 dw(4) = (one-ixr)*jxr
2735 WRITE(*,
'(A,2E24.16)')
'W3GRMP_R8 - REMAP (TGT):',xt,yt
2737 WRITE(*,
'(A,3I6,E24.16)')
'W3GRMP_R8 - REMAP (SRC):', &
2743 IF ( .NOT.ingrid .AND. fncl )
THEN
2745 ingrid = abs(ixr-half).LE.dd .AND. abs(jxr-half).LE.dd
2747 IF ( .NOT.ingrid )
RETURN
2749 IF ( .NOT.
PRESENT(mask) )
RETURN
2757 msk(l) = mask(is(l)-lbx+1,js(l)-lby+1)
2761 msk(l) = mask(js(l)-lby+1,is(l)-lbx+1)
2775 IF ( ns .EQ. 4 )
THEN
2779 IF ( ns .GT. 0 .AND. dsum .GT. small )
THEN
2783 WRITE(*,
'(A,2E24.16,4(2I6,E24.16))') &
2784 'W3GRMP_R8 - PARTIAL MASKED CELL:', &
2785 xt,yt,(is(l),js(l),dw(l),l=1,4)
2790 IF ( .NOT.
PRESENT(nnbr) )
RETURN
2799 dd = w3dist(llg,xt,yt,xs(l),ys(l))
2800 IF ( dd .LT. dmin )
THEN
2801 dmin = dd; icc = is(l); jcc = js(l);
2808 WRITE(*,
'(A,2I6)') &
2809 'W3GRMP_R8 - BEGIN POINT NNBR SEARCH:',icc,jcc
2811 level_loop:
DO lvl=0,nnp%NLVL
2812 nnbr_loop:
DO n=nnp%N1(lvl),nnp%N2(lvl)
2813 i = icc + nnp%DI(n); j = jcc + nnp%DJ(n);
2815 IF ( i.LT.lbx .OR. i.GT.ubx ) cycle nnbr_loop
2816 IF ( j.LT.lby .OR. j.GT.uby ) cycle nnbr_loop
2819 IF ( mod(iclo,2).EQ.0 ) &
2820 i = lbx + mod(nx - 1 + mod(i - lbx + 1, nx), nx)
2821 IF ( mod(iclo,3).EQ.0 ) &
2822 j = lby + mod(ny - 1 + mod(j - lby + 1, ny), ny)
2823 IF ( iclo.EQ.
iclo_trpl .AND. j.GT.uby )
THEN
2829 m = mask(i-lbx+1,j-lby+1)
2831 m = mask(j-lby+1,i-lbx+1)
2834 WRITE(*,
'(A,4I6,1L6)') &
2835 'W3GRMP_R8 - POINT NNBR SEARCH:',lvl,n,i,j,m
2837 IF ( m ) cycle nnbr_loop
2840 IF ( gkind.EQ.4 )
THEN
2841 x = xg4(i,j); y = yg4(i,j);
2843 x = xg8(i,j); y = yg8(i,j);
2846 IF ( gkind.EQ.4 )
THEN
2847 x = xg4(j,i); y = yg4(j,i);
2849 x = xg8(j,i); y = yg8(j,i);
2852 dd = w3dist(llg,xt,yt,x,y)
2854 IF ( ns .LT. nnbr )
THEN
2857 is(ns) = i; js(ns) = j; d(ns) = dd;
2859 IF ( ns .EQ. nnbr )
CALL w3sort(ns,is,js,d)
2863 CALL w3isrt(i,j,dd,ns,is,js,d)
2866 WRITE(*,
'(A,I2,I3,I6,4(2I6,E24.16))') &
2867 'W3GRMP_R8 - POINT NNBR LIST:', &
2868 lvl,n,ns,(is(l),js(l),d(l),l=1,ns)
2871 IF ( ns .EQ. nnbr )
EXIT level_loop
2876 IF ( nnbr .EQ. 0 )
RETURN
2881 dsum = dsum + one/(d(l)+small)
2883 dw(1:nnbr) = one/(d(1:nnbr)+small)/dsum
2886 WRITE(*,
'(A,2E24.16,I6)') &
2887 'W3GRMP_R8 - FULLY MASKED CELL (TGT):',xt,yt,nnbr
2889 WRITE(*,
'(A,3I6,E24.16)') &
2890 'W3GRMP_R8 - FULLY MASKED CELL (SRC):', &
2895 END FUNCTION w3grmp_r8
2993 FUNCTION w3grmc_r4( GSU, XTIN, YTIN, RTYP, NS, IS, JS, CS, EPS, &
2994 DCIN, WDTH, MASK, NMSK, DEBUG )
RESULT(INGRID)
2997 TYPE(
t_gsu),
INTENT(IN) :: gsu
2998 REAL(4),
INTENT(IN) :: xtin
2999 REAL(4),
INTENT(IN) :: ytin
3000 CHARACTER(6),
INTENT(IN):: rtyp
3001 INTEGER,
INTENT(OUT) :: ns
3002 INTEGER,
INTENT(INOUT),
POINTER :: is(:)
3003 INTEGER,
INTENT(INOUT),
POINTER :: js(:)
3004 REAL(4),
INTENT(INOUT),
POINTER :: cs(:)
3005 REAL(4),
INTENT(IN) ,
OPTIONAL :: eps
3006 REAL(4),
INTENT(IN) ,
OPTIONAL :: dcin
3007 REAL(4),
INTENT(IN) ,
OPTIONAL :: wdth
3008 LOGICAL,
INTENT(IN) ,
OPTIONAL :: mask(:,:)
3009 INTEGER,
INTENT(IN) ,
OPTIONAL :: nmsk
3010 LOGICAL,
INTENT(IN) ,
OPTIONAL :: debug
3013 REAL(8) :: leps, ldcin, lwdth=zero
3015 REAL(8),
POINTER :: cs8(:) => null()
3017 INTEGER,
SAVE :: ient = 0
3018 CALL strace (ient,
'W3GRMC_R4')
3024 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
3025 WRITE(0,
'(/2A/)')
'W3GRMC_R4 ERROR -- ', &
3026 'grid search utility object not created'
3035 IF ( .NOT.
PRESENT(wdth) )
THEN
3036 WRITE(0,
'(/2A/)')
'W3GRMC_R4 ERROR -- ', &
3037 'WDTH parameter is required with RTYP = filter'
3043 WRITE(0,
'(/2A/)')
'W3GRMC_R4 ERROR -- ', &
3044 'RTYP = '//rtyp//
' not supported'
3048 IF (
PRESENT(eps) )
THEN
3049 IF ( eps .LT. zero )
THEN
3050 WRITE(0,
'(/2A/)')
'W3GRMC_R4 ERROR -- ', &
3051 'EPS parameter must be >= 0'
3059 IF (
PRESENT(dcin) )
THEN
3060 IF ( dcin .LT. zero )
THEN
3061 WRITE(0,
'(/2A/)')
'W3GRMC_R4 ERROR -- ', &
3062 'DCIN parameter must be >= 0'
3073 xt = xtin; yt = ytin;
3074 ingrid = w3grmc( gsu, xt, yt, rtyp, ns, is, js, cs8, &
3075 eps=leps, dcin=ldcin, wdth=lwdth, &
3076 mask=mask, nmsk=nmsk, debug=debug )
3083 END FUNCTION w3grmc_r4
3087 FUNCTION w3grmc_r8( GSU, XTIN, YTIN, RTYP, NS, IS, JS, CS, EPS, &
3088 DCIN, WDTH, MASK, NMSK, DEBUG )
RESULT(INGRID)
3091 TYPE(
t_gsu),
INTENT(IN) :: gsu
3092 REAL(8),
INTENT(IN) :: xtin
3093 REAL(8),
INTENT(IN) :: ytin
3094 CHARACTER(6),
INTENT(IN):: rtyp
3095 INTEGER,
INTENT(OUT) :: ns
3096 INTEGER,
INTENT(INOUT),
POINTER :: is(:)
3097 INTEGER,
INTENT(INOUT),
POINTER :: js(:)
3098 REAL(8),
INTENT(INOUT),
POINTER :: cs(:)
3099 REAL(8),
INTENT(IN) ,
OPTIONAL :: eps
3100 REAL(8),
INTENT(IN) ,
OPTIONAL :: dcin
3101 REAL(8),
INTENT(IN) ,
OPTIONAL :: wdth
3102 LOGICAL,
INTENT(IN) ,
OPTIONAL :: mask(:,:)
3103 INTEGER,
INTENT(IN) ,
OPTIONAL :: nmsk
3104 LOGICAL,
INTENT(IN) ,
OPTIONAL :: debug
3107 LOGICAL,
PARAMETER :: lcmp = .true.
3108 INTEGER,
PARAMETER :: nmsk_default = 2
3109 REAL(8),
PARAMETER :: big = 1d16
3110 REAL(8) :: leps, lwdth=zero
3111 LOGICAL :: ldbg, fncl, pole, doblc, lmsk
3112 INTEGER :: i, ii, jj, k, kk, mcs, mcsmax
3113 INTEGER :: ic(4), jc(4)
3114 REAL(8) :: xt, yt, xc(4), yc(4)
3115 REAL(8) :: xtc, ytc, xsc(4), ysc(4)
3116 REAL(8) :: ldcin, ixr, jxr, dd, lon0, lat0, dmin
3117 REAL(8) :: ix, jx, czs
3119 LOGICAL,
POINTER :: lz(:)=>null()
3120 INTEGER,
POINTER :: iz(:)=>null(), jz(:)=>null()
3121 REAL(8),
POINTER :: cz(:)=>null()
3122 LOGICAL :: ijg, llg, lclo
3123 INTEGER :: iclo, gkind
3124 INTEGER :: lbx, lby, ubx, uby, nx, ny
3126 INTEGER,
SAVE :: ient = 0
3127 CALL strace (ient,
'W3GRMC_R8')
3133 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
3134 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3135 'grid search utility object not created'
3144 IF ( .NOT.
PRESENT(wdth) )
THEN
3145 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3146 'WDTH parameter is required with RTYP = filter'
3152 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3153 'RTYP = '//rtyp//
' not supported'
3157 IF (
PRESENT(eps) )
THEN
3158 IF ( eps .LT. zero )
THEN
3159 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3160 'EPS parameter must be >= 0'
3168 IF (
PRESENT(dcin) )
THEN
3169 IF ( dcin .LT. zero )
THEN
3170 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3171 'DCIN parameter must be >= 0'
3179 IF (
PRESENT(nmsk) )
THEN
3180 IF ( nmsk .LT. zero .OR. nmsk .GE. 4 )
THEN
3181 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3182 'NMSK parameter must be >= 0 and < 4'
3187 mcsmax = nmsk_default
3193 IF (
PRESENT(debug) )
THEN
3204 gkind = gsu%PTR%GKIND
3205 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
3206 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
3207 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
3209 IF (
PRESENT(mask) )
THEN
3211 IF ( .NOT.(ubound(mask,1).EQ.nx.AND. &
3212 ubound(mask,2).EQ.ny) )
THEN
3213 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3214 'MASK array size does not agree with GSU index bounds'
3218 IF ( .NOT.(ubound(mask,2).EQ.nx.AND. &
3219 ubound(mask,1).EQ.ny) )
THEN
3220 WRITE(0,
'(/2A/)')
'W3GRMC_R8 ERROR -- ', &
3221 'MASK array size does not agree with GSU index bounds'
3228 IF (
ASSOCIATED(is) )
THEN
3232 IF (
ASSOCIATED(js) )
THEN
3236 IF (
ASSOCIATED(cs) )
THEN
3241 xt = xtin; yt = ytin;
3242 IF ( ldbg )
WRITE(*,
'(/A,2E24.16)')
'W3GRMC_R8 - TARGET POINT:',xt,yt
3247 fncl = ldcin .GT. zero
3248 ingrid = w3gfcl(gsu,xt,yt,ic,jc,xc,yc,pole=pole,eps=leps,fncl=fncl,debug=ldbg)
3249 IF ( .NOT.ingrid .AND. .NOT.fncl )
RETURN
3252 lon0 = sum(xc)/four; lat0 = sum(yc)/four;
3253 IF ( d90-abs(lat0).GT.near_pole )
THEN
3255 CALL getpqr(xt,yt,xc,yc,ixr,jxr,eps=leps,debug=ldbg)
3258 CALL w3splx(lon0,lat0,zero,xt,yt,xtc,ytc)
3260 CALL w3splx(lon0,lat0,zero,xc(i),yc(i),xsc(i),ysc(i))
3262 CALL getpqr(xtc,ytc,xsc,ysc,ixr,jxr,eps=leps,debug=ldbg)
3265 WRITE(*,
'(A,2L2,2E24.16)')
'W3GRMC_R8 - RELATIVE:',ingrid,fncl,ixr,jxr
3268 IF ( .NOT.ingrid .AND. fncl )
THEN
3270 ingrid = abs(ixr-half).LE.dd .AND. abs(jxr-half).LE.dd
3272 IF ( .NOT.ingrid )
RETURN
3275 ix = ic(1) + ixr; jx = jc(1) + jxr;
3280 IF ( abs(ic(kk)-ix).LE.leps .AND. &
3281 abs(jc(kk)-jx).LE.leps )
THEN
3282 IF (
PRESENT(mask) )
THEN
3284 IF ( .NOT.mask(ic(kk)-lbx+1,jc(kk)-lby+1) )
EXIT kk_loop
3286 IF ( .NOT.mask(jc(kk)-lby+1,ic(kk)-lbx+1) )
EXIT kk_loop
3296 IF (
PRESENT(mask) )
THEN
3299 IF ( mask(ic(k)-lbx+1,jc(k)-lby+1) ) mcs = mcs+1
3301 IF ( mask(jc(k)-lby+1,ic(k)-lbx+1) ) mcs = mcs+1
3314 dd = (ix - ic(k))**2 + (jx - jc(k))**2
3315 IF ( dd .LT. dmin )
THEN
3316 dmin = dd; ii = ic(k); jj = jc(k);
3320 IF (
PRESENT(mask) )
THEN
3322 IF ( mask(ii-lbx+1,jj-lby+1) ) nz = 0
3324 IF ( mask(jj-lby+1,ii-lbx+1) ) nz = 0
3330 ALLOCATE( lz(nz), iz(nz), jz(nz), cz(nz) )
3347 ALLOCATE( lz(nz), iz(nz), jz(nz), cz(nz) )
3354 IF ( mcs.LE.mcsmax )
THEN
3357 CALL getblc( gsu, ic(1), jc(1), ixr, jxr, &
3358 lcmp, nz, lz, iz, jz, cz )
3372 ALLOCATE( lz(nz), iz(nz), jz(nz), cz(nz) )
3379 IF ( mcs.EQ.0 )
THEN
3382 CALL getbcc( gsu, ic(1), jc(1), ixr, jxr, &
3383 lcmp, nz, lz, iz, jz, cz )
3386 IF (
PRESENT(mask) )
THEN
3390 lmsk = mask(iz(k)-lbx+1,jz(k)-lby+1)
3392 lmsk = mask(jz(k)-lby+1,iz(k)-lbx+1)
3404 CALL getblc( gsu, ic(1), jc(1), ixr, jxr, &
3405 lcmp, nz, lz, iz, jz, cz )
3407 ELSE IF ( mcs.LE.mcsmax )
THEN
3410 CALL getblc( gsu, ic(1), jc(1), ixr, jxr, &
3411 lcmp, nz, lz, iz, jz, cz )
3421 IF ( mcs.LE.mcsmax )
THEN
3424 CALL getgfc( gsu, ic(1), jc(1), ixr, jxr, &
3425 lwdth, lcmp, nz, lz, iz, jz, cz )
3437 IF ( nz .GT. 1 )
THEN
3441 IF (
PRESENT(mask) )
THEN
3443 lmsk = mask(iz(k)-lbx+1,jz(k)-lby+1)
3445 lmsk = mask(jz(k)-lby+1,iz(k)-lbx+1)
3458 IF ( czs .GT. zero )
THEN
3460 IF ( lz(k) ) cz(k) = cz(k)/czs
3470 IF ( lz(k) ) ns = ns + 1
3473 ALLOCATE( is(ns), js(ns), cs(ns) )
3485 DEALLOCATE( lz, iz, jz, cz )
3487 END FUNCTION w3grmc_r8
3576 FUNCTION w3ckcl_r4( LLG, XT, YT, NS, XS, YS, POLE, EPS, DEBUG ) &
3580 LOGICAL,
INTENT(IN) :: llg
3581 REAL(4),
INTENT(INOUT) :: xt, yt
3582 INTEGER,
INTENT(IN) :: ns
3583 REAL(4),
INTENT(INOUT) :: xs(ns), ys(ns)
3584 LOGICAL,
INTENT(OUT) :: pole
3585 REAL(4),
INTENT(IN),
OPTIONAL :: eps
3586 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
3589 REAL(8) :: xt8, yt8, xs8(ns), ys8(ns), eps8
3591 INTEGER,
SAVE :: ient = 0
3592 CALL strace (ient,
'W3CKCL_R4')
3598 IF (
PRESENT(eps) )
THEN
3605 incell = w3ckcl( llg, xt8, yt8, ns, xs8, ys8, pole, &
3606 eps=eps8, debug=debug )
3611 END FUNCTION w3ckcl_r4
3615 FUNCTION w3ckcl_r8( LLG, XT, YT, NS, XS, YS, POLE, EPS, DEBUG ) &
3619 LOGICAL,
INTENT(IN) :: llg
3620 REAL(8),
INTENT(INOUT) :: xt, yt
3621 INTEGER,
INTENT(IN) :: ns
3622 REAL(8),
INTENT(INOUT) :: xs(ns), ys(ns)
3623 LOGICAL,
INTENT(OUT) :: pole
3624 REAL(8),
INTENT(IN),
OPTIONAL :: eps
3625 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
3629 LOGICAL :: ldbg, lsbc, bcut
3630 INTEGER :: i, j, k, n
3631 REAL(8) :: xxt, yyt, xxs(ns), yys(ns)
3632 REAL(8) :: xct, yct, xcs(ns), ycs(ns)
3633 REAL(8) :: v1x, v1y, v2x, v2y, s90
3637 INTEGER,
SAVE :: ient = 0
3638 CALL strace (ient,
'W3CKCL_R8')
3644 IF ( ns .LT. 3 )
THEN
3649 IF (
PRESENT(eps) )
THEN
3650 IF ( eps .LT. zero )
THEN
3651 WRITE(0,
'(/2A/)')
'W3CKCL_R8 ERROR -- ', &
3652 'EPS parameter must be >= 0'
3659 IF (
PRESENT(debug) )
THEN
3675 IF ( abs(xxs(j)-xxs(i)) .GT. d180 ) n = n + 1
3681 pole = n.EQ.1 .OR. count(abs(d90-abs(yys)).LE.leps).EQ.1
3689 IF ( minval(xxs) .GE. zero )
THEN
3690 WHERE ( xxs .GT. d180 ) xxs = xxs - d360
3691 IF ( xxt .GT. d180 ) xxt = xxt - d360
3693 WHERE ( xxs .LT. zero ) xxs = xxs + d360
3694 IF ( xxt .LT. zero ) xxt = xxt + d360
3697 WRITE(*,
'(A)')
'W3CKCL_R8 - CELL INCLUDES A BRANCH CUT'
3698 WRITE(*,
'(A,2E24.16,4(/A,1I1,A,2E24.16))') &
3699 'W3CKCL_R8 - SHIFT BRANCH CUT:',xxt,yyt, &
3700 (
' CORNER(',k,
'):',xxs(k),yys(k),k=1,4)
3708 IF ( abs(xxt-xxs(i)).LE.leps .AND. abs(yyt-yys(i)).LE.leps )
THEN
3710 WRITE(*,
'(A,I1,A,2E24.16)') &
3711 'W3CKCL_R8 - COINCIDENT WITH CORNER(',i,
'): ', &
3712 abs(xxt-xxs(i)),abs(yyt-yys(i))
3726 WRITE(*,
'(A)')
'W3CKCL_R8 - CELL INCLUDES A POLE'
3727 s90 = d90;
IF ( maxval(ys).LT.zero ) s90 = -d90;
3728 subcell_loop:
DO i=1,ns
3736 v1x = xxs(j) - xxs(i)
3737 v1y = yys(j) - yys(i)
3743 v1x = xxs(j) - xxs(j)
3750 v1x = xxs(i) - xxs(j)
3757 v1x = xxs(i) - xxs(i)
3764 IF ( abs(v1x) .GT. d180 )
THEN
3765 v1x = v1x - sign(d360,v1x)
3767 IF ( abs(v2x) .GT. d180 )
THEN
3768 v2x = v2x - sign(d360,v2x)
3771 cross = v1x*v2y - v1y*v2x
3773 IF ( abs(cross) .LT. leps ) cross = zero
3775 WRITE(*,
'(A,3(I1,A),5E24.16)')
'W3CKCL_R8 - CROSS(', &
3776 i,
',',j,
',',k,
'):',v1x,v1y,v2x,v2y,cross
3779 IF ( abs(sign1) .LE. leps )
THEN
3780 IF (abs(cross) .GT. leps) sign1 = sign(one,cross)
3784 IF ( abs(cross) .GT. leps )
THEN
3785 IF ( sign(one,cross) .NE. sign1 )
THEN
3798 xct = xxt; yct = yyt;
3799 xcs = xxs; ycs = yys;
3807 v1x = xcs(j) - xcs(i)
3808 v1y = ycs(j) - ycs(i)
3813 cross = v1x*v2y - v1y*v2x
3815 IF ( abs(cross) .LT. leps ) cross = zero
3817 WRITE(*,
'(A,2(I1,A),5E24.16)')
'W3CKCL_R8 - CROSS(', &
3818 i,
',',j,
'):',v1x,v1y,v2x,v2y,cross
3821 IF ( abs(sign1) .LE. leps )
THEN
3822 IF (abs(cross) .GT. leps) sign1 = sign(one,cross)
3826 IF ( abs(cross) .GT. leps )
THEN
3827 IF ( sign(one,cross) .NE. sign1 )
THEN
3840 END FUNCTION w3ckcl_r8
4014 SUBROUTINE w3cgdm_r4( IJG, LLG, ICLO, PTILED, QTILED, &
4015 PRANGE, QRANGE, LBI, UBI, LBO, UBO, X, Y, &
4016 MASK, NFD, SPHERE, RADIUS, DX, DY, &
4017 GPPC, GQQC, GPQC, GSQR, &
4018 HPFC, HQFC, APPC, AQQC, APQC, &
4019 DXDP, DYDP, DXDQ, DYDQ, &
4020 DPDX, DPDY, DQDX, DQDY, &
4021 COSA, COSC, SINC, ANGL, RC )
4023 LOGICAL,
INTENT(IN) :: ijg
4024 LOGICAL,
INTENT(IN) :: llg
4025 INTEGER,
INTENT(IN) :: iclo
4026 LOGICAL,
INTENT(IN) :: ptiled, qtiled
4027 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
4028 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
4029 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
4030 REAL(4),
INTENT(IN) :: x(lbi(1):ubi(1),lbi(2):ubi(2))
4031 REAL(4),
INTENT(IN) :: y(lbi(1):ubi(1),lbi(2):ubi(2))
4032 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
4033 INTEGER,
INTENT(IN),
OPTIONAL :: nfd
4034 LOGICAL,
INTENT(IN),
OPTIONAL :: sphere
4035 REAL(4),
INTENT(IN),
OPTIONAL :: radius
4036 REAL(4),
INTENT(IN),
OPTIONAL :: dx, dy
4037 REAL(4),
INTENT(OUT),
OPTIONAL :: gppc(lbo(1):ubo(1),lbo(2):ubo(2))
4038 REAL(4),
INTENT(OUT),
OPTIONAL :: gqqc(lbo(1):ubo(1),lbo(2):ubo(2))
4039 REAL(4),
INTENT(OUT),
OPTIONAL :: gpqc(lbo(1):ubo(1),lbo(2):ubo(2))
4040 REAL(4),
INTENT(OUT),
OPTIONAL :: gsqr(lbo(1):ubo(1),lbo(2):ubo(2))
4041 REAL(4),
INTENT(OUT),
OPTIONAL :: hpfc(lbo(1):ubo(1),lbo(2):ubo(2))
4042 REAL(4),
INTENT(OUT),
OPTIONAL :: hqfc(lbo(1):ubo(1),lbo(2):ubo(2))
4043 REAL(4),
INTENT(OUT),
OPTIONAL :: appc(lbo(1):ubo(1),lbo(2):ubo(2))
4044 REAL(4),
INTENT(OUT),
OPTIONAL :: aqqc(lbo(1):ubo(1),lbo(2):ubo(2))
4045 REAL(4),
INTENT(OUT),
OPTIONAL :: apqc(lbo(1):ubo(1),lbo(2):ubo(2))
4046 REAL(4),
INTENT(OUT),
OPTIONAL :: dxdp(lbo(1):ubo(1),lbo(2):ubo(2))
4047 REAL(4),
INTENT(OUT),
OPTIONAL :: dydp(lbo(1):ubo(1),lbo(2):ubo(2))
4048 REAL(4),
INTENT(OUT),
OPTIONAL :: dxdq(lbo(1):ubo(1),lbo(2):ubo(2))
4049 REAL(4),
INTENT(OUT),
OPTIONAL :: dydq(lbo(1):ubo(1),lbo(2):ubo(2))
4050 REAL(4),
INTENT(OUT),
OPTIONAL :: dpdx(lbo(1):ubo(1),lbo(2):ubo(2))
4051 REAL(4),
INTENT(OUT),
OPTIONAL :: dpdy(lbo(1):ubo(1),lbo(2):ubo(2))
4052 REAL(4),
INTENT(OUT),
OPTIONAL :: dqdx(lbo(1):ubo(1),lbo(2):ubo(2))
4053 REAL(4),
INTENT(OUT),
OPTIONAL :: dqdy(lbo(1):ubo(1),lbo(2):ubo(2))
4054 REAL(4),
INTENT(OUT),
OPTIONAL :: cosa(lbo(1):ubo(1),lbo(2):ubo(2))
4055 REAL(4),
INTENT(OUT),
OPTIONAL :: cosc(lbo(1):ubo(1),lbo(2):ubo(2))
4056 REAL(4),
INTENT(OUT),
OPTIONAL :: sinc(lbo(1):ubo(1),lbo(2):ubo(2))
4057 REAL(4),
INTENT(OUT),
OPTIONAL :: angl(lbo(1):ubo(1),lbo(2):ubo(2))
4058 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
4061 INTEGER,
PARAMETER :: m = 1
4062 REAL(8),
PARAMETER :: small = 1d-15
4063 INTEGER :: istat=0, n, np, nq, i1, i2, p, q
4065 REAL(8) :: r, facx, facy
4066 INTEGER,
ALLOCATABLE :: k(:,:,:), k2(:,:,:)
4067 REAL(8),
ALLOCATABLE :: c(:,:,:), c2(:,:,:)
4068 REAL(8) :: gppcl, gqqcl, gpqcl
4069 REAL(8) :: gsqrl, hpfcl, hqfcl
4070 REAL(8) :: appcl, aqqcl, apqcl
4071 REAL(8) :: dxdpl, dydpl, dxdql, dydql
4072 REAL(8) :: dpdxl, dpdyl, dqdxl, dqdyl
4073 REAL(8) :: cosal, sinal, costp, sintp, coscl, sincl, angll
4075 INTEGER,
SAVE :: ient = 0
4076 CALL strace (ient,
'W3CGDM_R4')
4081 IF (
PRESENT(rc) ) rc = 0
4083 IF (
PRESENT(nfd) )
THEN
4088 IF ( n.LE.0 .OR. mod(n,2).NE.0 )
THEN
4089 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4090 'NFD must be even and greater than zero'
4092 IF (
PRESENT(rc) )
THEN
4100 np = prange(2) - prange(1) + 1
4101 nq = qrange(2) - qrange(1) + 1
4103 SELECT CASE ( iclo )
4107 WRITE(0,
'(/1A,1A,1I2/)')
'W3CGDM ERROR -- ', &
4108 'unsupported ICLO: ',iclo
4110 IF (
PRESENT(rc) )
THEN
4118 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
4119 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4120 'tripole grid closure requires NP even'
4122 IF (
PRESENT(rc) )
THEN
4130 IF (
PRESENT(sphere) )
THEN
4136 IF (
PRESENT(radius) )
THEN
4143 IF (
PRESENT(dx) )
THEN
4144 IF ( dx.LE.zero )
THEN
4145 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ',
'DX must be > 0'
4147 IF (
PRESENT(rc) )
THEN
4156 IF (
PRESENT(dy) )
THEN
4157 IF ( dy.LE.zero )
THEN
4158 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ',
'DY must be > 0'
4160 IF (
PRESENT(rc) )
THEN
4172 ALLOCATE ( k(0:n,0:n,1:n), c(0:n,0:n,1:n), stat=istat )
4173 IF ( istat .NE. 0 )
THEN
4174 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4175 'finite difference coeff allocation failed'
4176 IF (
PRESENT(rc) )
THEN
4183 CALL get_fdw3 ( n, m, k, c )
4185 ALLOCATE ( k2(0:2,0:2,1:2), c2(0:2,0:2,1:2), stat=istat )
4186 IF ( istat .NE. 0 )
THEN
4187 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4188 'finite difference coeff allocation for N=2 failed'
4189 IF (
PRESENT(rc) )
THEN
4196 CALL get_fdw3 ( 2, m, k2, c2 )
4201 DO i2 = lbo(2), ubo(2)
4202 DO i1 = lbo(1), ubo(1)
4210 IF (
PRESENT(dx) )
THEN
4214 CALL dxydp( n, k, c, ijg, llg, iclo, ptiled, qtiled, &
4215 prange, qrange, lbi, ubi, p, q, dxdpl, dydpl, &
4216 mask=mask, x4=x, y4=y, rc=istat )
4217 IF ( istat .NE. 0 )
THEN
4218 IF (
PRESENT(rc) )
THEN
4226 IF (
PRESENT(dy) )
THEN
4230 CALL dxydq( n, k, c, ijg, llg, iclo, ptiled, qtiled, &
4231 prange, qrange, lbi, ubi, p, q, dxdql, dydql, &
4232 mask=mask, x4=x, y4=y, rc=istat )
4233 IF ( istat .NE. 0 )
THEN
4234 IF (
PRESENT(rc) )
THEN
4242 IF ( llg .AND. sphr )
THEN
4243 facx = facy*cos(real(y(i1,i2),8)*d2r)
4249 gsqrl = dxdpl*dydql - dxdql*dydpl
4250 IF ( gsqrl .LT. zero .AND. n .GT. 2 )
THEN
4255 IF (
PRESENT(dx) )
THEN
4259 CALL dxydp( 2, k2, c2, ijg, llg, iclo, ptiled, qtiled, &
4260 prange, qrange, lbi, ubi, p, q, dxdpl, dydpl, &
4261 mask=mask, x4=x, y4=y, rc=istat )
4262 IF ( istat .NE. 0 )
THEN
4263 IF (
PRESENT(rc) )
THEN
4271 IF (
PRESENT(dy) )
THEN
4275 CALL dxydq( 2, k2, c2, ijg, llg, iclo, ptiled, qtiled, &
4276 prange, qrange, lbi, ubi, p, q, dxdql, dydql, &
4277 mask=mask, x4=x, y4=y, rc=istat )
4278 IF ( istat .NE. 0 )
THEN
4279 IF (
PRESENT(rc) )
THEN
4287 IF ( llg .AND. sphr )
THEN
4288 facx = facy*cos(real(y(i1,i2),8)*d2r)
4294 gsqrl = dxdpl*dydql - dxdql*dydpl
4296 IF ( gsqrl .LT. zero )
THEN
4298 WRITE(0,
'(/1A,1A)')
'W3CGDM ERROR -- ', &
4299 'input coordinates do not define a '// &
4300 'right-handed coordinate system'
4301 WRITE(0,
'(1A,2A6,5A16)')
'W3CGDM ERROR --', &
4302 'P',
'Q',
'GSQRL',
'DXDPL',
'DYDQL',
'DXDQL',
'DYDPL'
4303 WRITE(0,
'(1A,2I6,5E16.8/)')
'W3CGDM ERROR --', &
4304 p,q,gsqrl,dxdpl,dydql,dxdql,dydpl
4305 IF (
PRESENT(rc) )
THEN
4312 gppcl = dxdpl*dxdpl + dydpl*dydpl
4313 gqqcl = dxdql*dxdql + dydql*dydql
4314 gpqcl = dxdpl*dxdql + dydpl*dydql
4315 gsqrl = max(gsqrl,small)
4316 gppcl = max(gppcl,small)
4317 gqqcl = max(gqqcl,small)
4322 appcl = dpdxl*dpdxl + dpdyl*dpdyl
4323 aqqcl = dqdxl*dqdxl + dqdyl*dqdyl
4324 apqcl = dpdxl*dqdxl + dpdyl*dqdyl
4327 cosal = gpqcl/(hpfcl*hqfcl)
4328 sinal = gsqrl**2/(gppcl*gqqcl)
4331 coscl = sinal*costp + cosal*sintp
4332 sincl = sinal*sintp - cosal*costp
4333 angll = atan2(sincl,coscl)*r2d
4334 IF (
PRESENT(gppc)) gppc(i1,i2) = gppcl
4335 IF (
PRESENT(gqqc)) gqqc(i1,i2) = gqqcl
4336 IF (
PRESENT(gpqc)) gpqc(i1,i2) = gpqcl
4337 IF (
PRESENT(appc)) appc(i1,i2) = appcl
4338 IF (
PRESENT(aqqc)) aqqc(i1,i2) = aqqcl
4339 IF (
PRESENT(apqc)) apqc(i1,i2) = apqcl
4340 IF (
PRESENT(gsqr)) gsqr(i1,i2) = gsqrl
4341 IF (
PRESENT(hpfc)) hpfc(i1,i2) = hpfcl
4342 IF (
PRESENT(hqfc)) hqfc(i1,i2) = hqfcl
4343 IF (
PRESENT(dxdp)) dxdp(i1,i2) = dxdpl
4344 IF (
PRESENT(dydp)) dydp(i1,i2) = dydpl
4345 IF (
PRESENT(dxdq)) dxdq(i1,i2) = dxdql
4346 IF (
PRESENT(dydq)) dydq(i1,i2) = dydql
4347 IF (
PRESENT(dpdx)) dpdx(i1,i2) = dpdxl
4348 IF (
PRESENT(dpdy)) dpdy(i1,i2) = dpdyl
4349 IF (
PRESENT(dqdx)) dqdx(i1,i2) = dqdxl
4350 IF (
PRESENT(dqdy)) dqdy(i1,i2) = dqdyl
4351 IF (
PRESENT(cosa)) cosa(i1,i2) = cosal
4352 IF (
PRESENT(cosc)) cosc(i1,i2) = coscl
4353 IF (
PRESENT(sinc)) sinc(i1,i2) = sincl
4354 IF (
PRESENT(angl)) angl(i1,i2) = angll
4361 DEALLOCATE ( k, c, k2, c2 )
4363 END SUBROUTINE w3cgdm_r4
4367 SUBROUTINE w3cgdm_r8( IJG, LLG, ICLO, PTILED, QTILED, &
4368 PRANGE, QRANGE, LBI, UBI, LBO, UBO, X, Y, &
4369 MASK, NFD, SPHERE, RADIUS, DX, DY, &
4370 GPPC, GQQC, GPQC, GSQR, &
4371 HPFC, HQFC, APPC, AQQC, APQC, &
4372 DXDP, DYDP, DXDQ, DYDQ, &
4373 DPDX, DPDY, DQDX, DQDY, &
4374 COSA, COSC, SINC, ANGL, RC )
4376 LOGICAL,
INTENT(IN) :: ijg
4377 LOGICAL,
INTENT(IN) :: llg
4378 INTEGER,
INTENT(IN) :: iclo
4379 LOGICAL,
INTENT(IN) :: ptiled, qtiled
4380 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
4381 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
4382 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
4383 REAL(8),
INTENT(IN) :: x(lbi(1):ubi(1),lbi(2):ubi(2))
4384 REAL(8),
INTENT(IN) :: y(lbi(1):ubi(1),lbi(2):ubi(2))
4385 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
4386 INTEGER,
INTENT(IN),
OPTIONAL :: nfd
4387 LOGICAL,
INTENT(IN),
OPTIONAL :: sphere
4388 REAL(8),
INTENT(IN),
OPTIONAL :: radius
4389 REAL(8),
INTENT(IN),
OPTIONAL :: dx, dy
4390 REAL(8),
INTENT(OUT),
OPTIONAL :: gppc(lbo(1):ubo(1),lbo(2):ubo(2))
4391 REAL(8),
INTENT(OUT),
OPTIONAL :: gqqc(lbo(1):ubo(1),lbo(2):ubo(2))
4392 REAL(8),
INTENT(OUT),
OPTIONAL :: gpqc(lbo(1):ubo(1),lbo(2):ubo(2))
4393 REAL(8),
INTENT(OUT),
OPTIONAL :: gsqr(lbo(1):ubo(1),lbo(2):ubo(2))
4394 REAL(8),
INTENT(OUT),
OPTIONAL :: hpfc(lbo(1):ubo(1),lbo(2):ubo(2))
4395 REAL(8),
INTENT(OUT),
OPTIONAL :: hqfc(lbo(1):ubo(1),lbo(2):ubo(2))
4396 REAL(8),
INTENT(OUT),
OPTIONAL :: appc(lbo(1):ubo(1),lbo(2):ubo(2))
4397 REAL(8),
INTENT(OUT),
OPTIONAL :: aqqc(lbo(1):ubo(1),lbo(2):ubo(2))
4398 REAL(8),
INTENT(OUT),
OPTIONAL :: apqc(lbo(1):ubo(1),lbo(2):ubo(2))
4399 REAL(8),
INTENT(OUT),
OPTIONAL :: dxdp(lbo(1):ubo(1),lbo(2):ubo(2))
4400 REAL(8),
INTENT(OUT),
OPTIONAL :: dydp(lbo(1):ubo(1),lbo(2):ubo(2))
4401 REAL(8),
INTENT(OUT),
OPTIONAL :: dxdq(lbo(1):ubo(1),lbo(2):ubo(2))
4402 REAL(8),
INTENT(OUT),
OPTIONAL :: dydq(lbo(1):ubo(1),lbo(2):ubo(2))
4403 REAL(8),
INTENT(OUT),
OPTIONAL :: dpdx(lbo(1):ubo(1),lbo(2):ubo(2))
4404 REAL(8),
INTENT(OUT),
OPTIONAL :: dpdy(lbo(1):ubo(1),lbo(2):ubo(2))
4405 REAL(8),
INTENT(OUT),
OPTIONAL :: dqdx(lbo(1):ubo(1),lbo(2):ubo(2))
4406 REAL(8),
INTENT(OUT),
OPTIONAL :: dqdy(lbo(1):ubo(1),lbo(2):ubo(2))
4407 REAL(8),
INTENT(OUT),
OPTIONAL :: cosa(lbo(1):ubo(1),lbo(2):ubo(2))
4408 REAL(8),
INTENT(OUT),
OPTIONAL :: cosc(lbo(1):ubo(1),lbo(2):ubo(2))
4409 REAL(8),
INTENT(OUT),
OPTIONAL :: sinc(lbo(1):ubo(1),lbo(2):ubo(2))
4410 REAL(8),
INTENT(OUT),
OPTIONAL :: angl(lbo(1):ubo(1),lbo(2):ubo(2))
4411 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
4414 INTEGER,
PARAMETER :: m = 1
4415 REAL(8),
PARAMETER :: small = 1d-15
4416 INTEGER :: istat=0, n, np, nq, i1, i2, p, q
4418 REAL(8) :: r, facx, facy
4419 INTEGER,
ALLOCATABLE :: k(:,:,:), k2(:,:,:)
4420 REAL(8),
ALLOCATABLE :: c(:,:,:), c2(:,:,:)
4421 REAL(8) :: gppcl, gqqcl, gpqcl
4422 REAL(8) :: gsqrl, hpfcl, hqfcl
4423 REAL(8) :: appcl, aqqcl, apqcl
4424 REAL(8) :: dxdpl, dydpl, dxdql, dydql
4425 REAL(8) :: dpdxl, dpdyl, dqdxl, dqdyl
4426 REAL(8) :: cosal, sinal, costp, sintp, coscl, sincl, angll
4428 INTEGER,
SAVE :: ient = 0
4429 CALL strace (ient,
'W3CGDM_R8')
4434 IF (
PRESENT(rc) ) rc = 0
4436 IF (
PRESENT(nfd) )
THEN
4441 IF ( n.LE.0 .OR. mod(n,2).NE.0 )
THEN
4442 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4443 'NFD must be even and greater than zero'
4445 IF (
PRESENT(rc) )
THEN
4453 np = prange(2) - prange(1) + 1
4454 nq = qrange(2) - qrange(1) + 1
4456 SELECT CASE ( iclo )
4460 WRITE(0,
'(/1A,1A,1I2/)')
'W3CGDM ERROR -- ', &
4461 'unsupported ICLO: ',iclo
4463 IF (
PRESENT(rc) )
THEN
4471 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
4472 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4473 'tripole grid closure requires NP even'
4475 IF (
PRESENT(rc) )
THEN
4483 IF (
PRESENT(sphere) )
THEN
4489 IF (
PRESENT(radius) )
THEN
4496 IF (
PRESENT(dx) )
THEN
4497 IF ( dx.LE.zero )
THEN
4498 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ',
'DX must be > 0'
4500 IF (
PRESENT(rc) )
THEN
4509 IF (
PRESENT(dy) )
THEN
4510 IF ( dy.LE.zero )
THEN
4511 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ',
'DY must be > 0'
4513 IF (
PRESENT(rc) )
THEN
4525 ALLOCATE ( k(0:n,0:n,1:n), c(0:n,0:n,1:n), stat=istat )
4526 IF ( istat .NE. 0 )
THEN
4527 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4528 'finite difference coeff allocation failed'
4529 IF (
PRESENT(rc) )
THEN
4536 CALL get_fdw3 ( n, m, k, c )
4538 ALLOCATE ( k2(0:2,0:2,1:2), c2(0:2,0:2,1:2), stat=istat )
4539 IF ( istat .NE. 0 )
THEN
4540 WRITE(0,
'(/1A,1A/)')
'W3CGDM ERROR -- ', &
4541 'finite difference coeff allocation for N=2 failed'
4542 IF (
PRESENT(rc) )
THEN
4549 CALL get_fdw3 ( 2, m, k2, c2 )
4554 DO i2 = lbo(2), ubo(2)
4555 DO i1 = lbo(1), ubo(1)
4563 IF (
PRESENT(dx) )
THEN
4567 CALL dxydp( n, k, c, ijg, llg, iclo, ptiled, qtiled, &
4568 prange, qrange, lbi, ubi, p, q, dxdpl, dydpl, &
4569 mask=mask, x8=x, y8=y, rc=istat )
4570 IF ( istat .NE. 0 )
THEN
4571 IF (
PRESENT(rc) )
THEN
4579 IF (
PRESENT(dy) )
THEN
4583 CALL dxydq( n, k, c, ijg, llg, iclo, ptiled, qtiled, &
4584 prange, qrange, lbi, ubi, p, q, dxdql, dydql, &
4585 mask=mask, x8=x, y8=y, rc=istat )
4586 IF ( istat .NE. 0 )
THEN
4587 IF (
PRESENT(rc) )
THEN
4595 IF ( llg .AND. sphr )
THEN
4596 facx = facy*cos(real(y(i1,i2),8)*d2r)
4602 gsqrl = dxdpl*dydql - dxdql*dydpl
4603 IF ( gsqrl .LT. zero .AND. n .GT. 2 )
THEN
4608 IF (
PRESENT(dx) )
THEN
4612 CALL dxydp( 2, k2, c2, ijg, llg, iclo, ptiled, qtiled, &
4613 prange, qrange, lbi, ubi, p, q, dxdpl, dydpl, &
4614 mask=mask, x8=x, y8=y, rc=istat )
4615 IF ( istat .NE. 0 )
THEN
4616 IF (
PRESENT(rc) )
THEN
4624 IF (
PRESENT(dy) )
THEN
4628 CALL dxydq( 2, k2, c2, ijg, llg, iclo, ptiled, qtiled, &
4629 prange, qrange, lbi, ubi, p, q, dxdql, dydql, &
4630 mask=mask, x8=x, y8=y, rc=istat )
4631 IF ( istat .NE. 0 )
THEN
4632 IF (
PRESENT(rc) )
THEN
4640 IF ( llg .AND. sphr )
THEN
4641 facx = facy*cos(real(y(i1,i2),8)*d2r)
4647 gsqrl = dxdpl*dydql - dxdql*dydpl
4649 IF ( gsqrl .LT. zero )
THEN
4651 WRITE(0,
'(/1A,1A)')
'W3CGDM ERROR -- ', &
4652 'input coordinates do not define a '// &
4653 'right-handed coordinate system'
4654 WRITE(0,
'(1A,2A6,5A16)')
'W3CGDM ERROR --', &
4655 'P',
'Q',
'GSQRL',
'DXDPL',
'DYDQL',
'DXDQL',
'DYDPL'
4656 WRITE(0,
'(1A,2I6,5E16.8/)')
'W3CGDM ERROR --', &
4657 p,q,gsqrl,dxdpl,dydql,dxdql,dydpl
4658 IF (
PRESENT(rc) )
THEN
4665 gppcl = dxdpl*dxdpl + dydpl*dydpl
4666 gqqcl = dxdql*dxdql + dydql*dydql
4667 gpqcl = dxdpl*dxdql + dydpl*dydql
4668 gsqrl = max(gsqrl,small)
4669 gppcl = max(gppcl,small)
4670 gqqcl = max(gqqcl,small)
4675 appcl = dpdxl*dpdxl + dpdyl*dpdyl
4676 aqqcl = dqdxl*dqdxl + dqdyl*dqdyl
4677 apqcl = dpdxl*dqdxl + dpdyl*dqdyl
4680 cosal = gpqcl/(hpfcl*hqfcl)
4681 sinal = gsqrl**2/(gppcl*gqqcl)
4684 coscl = sinal*costp + cosal*sintp
4685 sincl = sinal*sintp - cosal*costp
4686 angll = atan2(sincl,coscl)*r2d
4687 IF (
PRESENT(gppc)) gppc(i1,i2) = gppcl
4688 IF (
PRESENT(gqqc)) gqqc(i1,i2) = gqqcl
4689 IF (
PRESENT(gpqc)) gpqc(i1,i2) = gpqcl
4690 IF (
PRESENT(appc)) appc(i1,i2) = appcl
4691 IF (
PRESENT(aqqc)) aqqc(i1,i2) = aqqcl
4692 IF (
PRESENT(apqc)) apqc(i1,i2) = apqcl
4693 IF (
PRESENT(gsqr)) gsqr(i1,i2) = gsqrl
4694 IF (
PRESENT(hpfc)) hpfc(i1,i2) = hpfcl
4695 IF (
PRESENT(hqfc)) hqfc(i1,i2) = hqfcl
4696 IF (
PRESENT(dxdp)) dxdp(i1,i2) = dxdpl
4697 IF (
PRESENT(dydp)) dydp(i1,i2) = dydpl
4698 IF (
PRESENT(dxdq)) dxdq(i1,i2) = dxdql
4699 IF (
PRESENT(dydq)) dydq(i1,i2) = dydql
4700 IF (
PRESENT(dpdx)) dpdx(i1,i2) = dpdxl
4701 IF (
PRESENT(dpdy)) dpdy(i1,i2) = dpdyl
4702 IF (
PRESENT(dqdx)) dqdx(i1,i2) = dqdxl
4703 IF (
PRESENT(dqdy)) dqdy(i1,i2) = dqdyl
4704 IF (
PRESENT(cosa)) cosa(i1,i2) = cosal
4705 IF (
PRESENT(cosc)) cosc(i1,i2) = coscl
4706 IF (
PRESENT(sinc)) sinc(i1,i2) = sincl
4707 IF (
PRESENT(angl)) angl(i1,i2) = atan2(sincl,coscl)*r2d
4708 IF (
PRESENT(angl)) angl(i1,i2) = angll
4715 DEALLOCATE ( k, c, k2, c2 )
4717 END SUBROUTINE w3cgdm_r8
4808 SUBROUTINE w3grd0_r4( NFD, IJG, ICLO, PTILED, QTILED, &
4809 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
4810 DPDX, DPDY, DQDX, DQDY, &
4811 F, DFDX, DFDY, MASK, RC )
4813 INTEGER,
INTENT(IN) :: nfd
4814 LOGICAL,
INTENT(IN) :: ijg
4815 INTEGER,
INTENT(IN) :: iclo
4816 LOGICAL,
INTENT(IN) :: ptiled, qtiled
4817 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
4818 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
4819 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
4820 REAL(4),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
4821 REAL(4),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
4822 REAL(4),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
4823 REAL(4),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
4824 REAL(4),
INTENT(IN) :: f(lbi(1):ubi(1),lbi(2):ubi(2))
4825 REAL(4),
INTENT(OUT) :: dfdx(lbo(1):ubo(1),lbo(2):ubo(2))
4826 REAL(4),
INTENT(OUT) :: dfdy(lbo(1):ubo(1),lbo(2):ubo(2))
4827 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
4828 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
4831 INTEGER,
PARAMETER :: m = 1
4832 INTEGER :: np, nq, i1, i2, p, q
4834 INTEGER :: k(0:nfd,0:nfd,1:nfd)
4835 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
4836 REAL(8) :: dfdp, dfdq
4838 INTEGER,
SAVE :: ient = 0
4839 CALL strace (ient,
'W3GRD0_R4')
4844 IF (
PRESENT(rc) ) rc = 0
4846 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
4847 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
4848 'NFD must be even and greater than zero'
4850 IF (
PRESENT(rc) )
THEN
4858 np = prange(2) - prange(1) + 1
4859 nq = qrange(2) - qrange(1) + 1
4861 SELECT CASE ( iclo )
4865 WRITE(0,
'(/1A,1A,1I2/)')
'W3GRD0 ERROR -- ', &
4866 'unsupported ICLO: ',iclo
4868 IF (
PRESENT(rc) )
THEN
4876 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
4877 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
4878 'tripole grid closure requires NP even'
4880 IF (
PRESENT(rc) )
THEN
4891 CALL get_fdw3 ( nfd, m, k, c )
4896 DO i2 = lbo(2), ubo(2)
4897 DO i1 = lbo(1), ubo(1)
4898 IF (
PRESENT(mask) )
THEN
4899 IF ( mask(i1,i2) ) cycle
4908 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
4909 prange, qrange, lbi, ubi, p, q, &
4910 f4=f, dfdp=dfdp, dfdq=dfdq, &
4911 mask=mask, rc=istat )
4912 IF ( istat .NE. 0 )
THEN
4913 IF (
PRESENT(rc) )
THEN
4920 dfdx(i1,i2) = dfdp*dpdx(i1,i2) + dfdq*dqdx(i1,i2)
4921 dfdy(i1,i2) = dfdp*dpdy(i1,i2) + dfdq*dqdy(i1,i2)
4925 END SUBROUTINE w3grd0_r4
4929 SUBROUTINE w3grd0_r8( NFD, IJG, ICLO, PTILED, QTILED, &
4930 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
4931 DPDX, DPDY, DQDX, DQDY, &
4932 F, DFDX, DFDY, MASK, RC )
4934 INTEGER,
INTENT(IN) :: nfd
4935 LOGICAL,
INTENT(IN) :: ijg
4936 INTEGER,
INTENT(IN) :: iclo
4937 LOGICAL,
INTENT(IN) :: ptiled, qtiled
4938 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
4939 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
4940 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
4941 REAL(8),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
4942 REAL(8),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
4943 REAL(8),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
4944 REAL(8),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
4945 REAL(8),
INTENT(IN) :: f(lbi(1):ubi(1),lbi(2):ubi(2))
4946 REAL(8),
INTENT(OUT) :: dfdx(lbo(1):ubo(1),lbo(2):ubo(2))
4947 REAL(8),
INTENT(OUT) :: dfdy(lbo(1):ubo(1),lbo(2):ubo(2))
4948 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
4949 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
4952 INTEGER,
PARAMETER :: m = 1
4953 INTEGER :: np, nq, i1, i2, p, q
4955 INTEGER :: k(0:nfd,0:nfd,1:nfd)
4956 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
4957 REAL(8) :: dfdp, dfdq
4959 INTEGER,
SAVE :: ient = 0
4960 CALL strace (ient,
'W3GRD0_R8')
4965 IF (
PRESENT(rc) ) rc = 0
4967 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
4968 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
4969 'NFD must be even and greater than zero'
4971 IF (
PRESENT(rc) )
THEN
4979 np = prange(2) - prange(1) + 1
4980 nq = qrange(2) - qrange(1) + 1
4982 SELECT CASE ( iclo )
4986 WRITE(0,
'(/1A,1A,1I2/)')
'W3GRD0 ERROR -- ', &
4987 'unsupported ICLO: ',iclo
4989 IF (
PRESENT(rc) )
THEN
4997 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
4998 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
4999 'tripole grid closure requires NP even'
5001 IF (
PRESENT(rc) )
THEN
5012 CALL get_fdw3 ( nfd, m, k, c )
5017 DO i2 = lbo(2), ubo(2)
5018 DO i1 = lbo(1), ubo(1)
5019 IF (
PRESENT(mask) )
THEN
5020 IF ( mask(i1,i2) ) cycle
5029 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
5030 prange, qrange, lbi, ubi, p, q, &
5031 f8=f, dfdp=dfdp, dfdq=dfdq, &
5032 mask=mask, rc=istat )
5033 IF ( istat .NE. 0 )
THEN
5034 IF (
PRESENT(rc) )
THEN
5041 dfdx(i1,i2) = dfdp*dpdx(i1,i2) + dfdq*dqdx(i1,i2)
5042 dfdy(i1,i2) = dfdp*dpdy(i1,i2) + dfdq*dqdy(i1,i2)
5046 END SUBROUTINE w3grd0_r8
5139 SUBROUTINE w3div1_r4( NFD, IJG, ICLO, PTILED, QTILED, &
5140 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
5141 DPDX, DPDY, DQDX, DQDY, &
5142 VX, VY, DIVV, MASK, RC )
5144 INTEGER,
INTENT(IN) :: nfd
5145 LOGICAL,
INTENT(IN) :: ijg
5146 INTEGER,
INTENT(IN) :: iclo
5147 LOGICAL,
INTENT(IN) :: ptiled, qtiled
5148 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
5149 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
5150 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
5151 REAL(4),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
5152 REAL(4),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
5153 REAL(4),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
5154 REAL(4),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
5155 REAL(4),
INTENT(IN) :: vx(lbi(1):ubi(1),lbi(2):ubi(2))
5156 REAL(4),
INTENT(IN) :: vy(lbi(1):ubi(1),lbi(2):ubi(2))
5157 REAL(4),
INTENT(OUT) :: divv(lbo(1):ubo(1),lbo(2):ubo(2))
5158 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
5159 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
5162 INTEGER,
PARAMETER :: m = 1
5163 INTEGER :: np, nq, i1, i2, p, q
5165 INTEGER :: k(0:nfd,0:nfd,1:nfd)
5166 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
5167 REAL(8) :: dvxdp, dvxdq, dvydp, dvydq
5168 REAL(8) :: dvxdx, dvydy
5170 INTEGER,
SAVE :: ient = 0
5171 CALL strace (ient,
'W3DIV1_R4')
5176 IF (
PRESENT(rc) ) rc = 0
5178 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
5179 WRITE(0,
'(/1A,1A/)')
'W3DIV1 ERROR -- ', &
5180 'NFD must be even and greater than zero'
5182 IF (
PRESENT(rc) )
THEN
5190 np = prange(2) - prange(1) + 1
5191 nq = qrange(2) - qrange(1) + 1
5193 SELECT CASE ( iclo )
5197 WRITE(0,
'(/1A,1A,1I2/)')
'W3DIV1 ERROR -- ', &
5198 'unsupported ICLO: ',iclo
5200 IF (
PRESENT(rc) )
THEN
5208 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
5209 WRITE(0,
'(/1A,1A/)')
'W3DIV1 ERROR -- ', &
5210 'tripole grid closure requires NP even'
5212 IF (
PRESENT(rc) )
THEN
5223 CALL get_fdw3 ( nfd, m, k, c )
5228 DO i2 = lbo(2), ubo(2)
5229 DO i1 = lbo(1), ubo(1)
5230 IF (
PRESENT(mask) )
THEN
5231 IF ( mask(i1,i2) ) cycle
5240 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
5241 prange, qrange, lbi, ubi, p, q, &
5242 f4=vx, dfdp=dvxdp, dfdq=dvxdq, &
5243 g4=vy, dgdp=dvydp, dgdq=dvydq, &
5244 mask=mask, rc=istat )
5245 IF ( istat .NE. 0 )
THEN
5246 IF (
PRESENT(rc) )
THEN
5253 dvxdx = dvxdp*dpdx(i1,i2) + dvxdq*dqdx(i1,i2)
5254 dvydy = dvydp*dpdy(i1,i2) + dvydq*dqdy(i1,i2)
5255 divv(i1,i2) = dvxdx + dvydy
5259 END SUBROUTINE w3div1_r4
5263 SUBROUTINE w3div1_r8( NFD, IJG, ICLO, PTILED, QTILED, &
5264 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
5265 DPDX, DPDY, DQDX, DQDY, &
5266 VX, VY, DIVV, MASK, RC )
5268 INTEGER,
INTENT(IN) :: nfd
5269 LOGICAL,
INTENT(IN) :: ijg
5270 INTEGER,
INTENT(IN) :: iclo
5271 LOGICAL,
INTENT(IN) :: ptiled, qtiled
5272 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
5273 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
5274 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
5275 REAL(8),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
5276 REAL(8),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
5277 REAL(8),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
5278 REAL(8),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
5279 REAL(8),
INTENT(IN) :: vx(lbi(1):ubi(1),lbi(2):ubi(2))
5280 REAL(8),
INTENT(IN) :: vy(lbi(1):ubi(1),lbi(2):ubi(2))
5281 REAL(8),
INTENT(OUT) :: divv(lbo(1):ubo(1),lbo(2):ubo(2))
5282 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
5283 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
5286 INTEGER,
PARAMETER :: m = 1
5287 INTEGER :: np, nq, i1, i2, p, q
5289 INTEGER :: k(0:nfd,0:nfd,1:nfd)
5290 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
5291 REAL(8) :: dvxdp, dvxdq, dvydp, dvydq
5292 REAL(8) :: dvxdx, dvydy
5294 INTEGER,
SAVE :: ient = 0
5295 CALL strace (ient,
'W3DIV1_R8')
5300 IF (
PRESENT(rc) ) rc = 0
5302 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
5303 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
5304 'NFD must be even and greater than zero'
5306 IF (
PRESENT(rc) )
THEN
5314 np = prange(2) - prange(1) + 1
5315 nq = qrange(2) - qrange(1) + 1
5317 SELECT CASE ( iclo )
5321 WRITE(0,
'(/1A,1A,1I2/)')
'W3GRD0 ERROR -- ', &
5322 'unsupported ICLO: ',iclo
5324 IF (
PRESENT(rc) )
THEN
5332 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
5333 WRITE(0,
'(/1A,1A/)')
'W3GRD0 ERROR -- ', &
5334 'tripole grid closure requires NP even'
5336 IF (
PRESENT(rc) )
THEN
5347 CALL get_fdw3 ( nfd, m, k, c )
5352 DO i2 = lbo(2), ubo(2)
5353 DO i1 = lbo(1), ubo(1)
5354 IF (
PRESENT(mask) )
THEN
5355 IF ( mask(i1,i2) ) cycle
5364 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
5365 prange, qrange, lbi, ubi, p, q, &
5366 f8=vx, dfdp=dvxdp, dfdq=dvxdq, &
5367 g8=vy, dgdp=dvydp, dgdq=dvydq, &
5368 mask=mask, rc=istat )
5369 IF ( istat .NE. 0 )
THEN
5370 IF (
PRESENT(rc) )
THEN
5377 dvxdx = dvxdp*dpdx(i1,i2) + dvxdq*dqdx(i1,i2)
5378 dvydy = dvydp*dpdy(i1,i2) + dvydq*dqdy(i1,i2)
5379 divv(i1,i2) = dvxdx + dvydy
5383 END SUBROUTINE w3div1_r8
5479 SUBROUTINE w3div2_r4( NFD, IJG, ICLO, PTILED, QTILED, &
5480 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
5481 DPDX, DPDY, DQDX, DQDY, &
5482 SXX, SYY, SXY, DSX, DSY, MASK, RC )
5484 INTEGER,
INTENT(IN) :: nfd
5485 LOGICAL,
INTENT(IN) :: ijg
5486 INTEGER,
INTENT(IN) :: iclo
5487 LOGICAL,
INTENT(IN) :: ptiled, qtiled
5488 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
5489 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
5490 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
5491 REAL(4),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
5492 REAL(4),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
5493 REAL(4),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
5494 REAL(4),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
5495 REAL(4),
INTENT(IN) :: sxx(lbi(1):ubi(1),lbi(2):ubi(2))
5496 REAL(4),
INTENT(IN) :: syy(lbi(1):ubi(1),lbi(2):ubi(2))
5497 REAL(4),
INTENT(IN) :: sxy(lbi(1):ubi(1),lbi(2):ubi(2))
5498 REAL(4),
INTENT(OUT) :: dsx(lbo(1):ubo(1),lbo(2):ubo(2))
5499 REAL(4),
INTENT(OUT) :: dsy(lbo(1):ubo(1),lbo(2):ubo(2))
5500 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
5501 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
5504 INTEGER,
PARAMETER :: m = 1
5505 INTEGER :: np, nq, i1, i2, p, q
5507 INTEGER :: k(0:nfd,0:nfd,1:nfd)
5508 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
5509 REAL(8) :: dxxdp, dxxdq, dyydp, dyydq, dxydp, dxydq
5510 REAL(8) :: dxxdx, dyydy, dxydx, dxydy
5512 INTEGER,
SAVE :: ient = 0
5513 CALL strace (ient,
'W3DIV2_R4')
5518 IF (
PRESENT(rc) ) rc = 0
5520 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
5521 WRITE(0,
'(/1A,1A/)')
'W3DIV2 ERROR -- ', &
5522 'NFD must be even and greater than zero'
5524 IF (
PRESENT(rc) )
THEN
5532 np = prange(2) - prange(1) + 1
5533 nq = qrange(2) - qrange(1) + 1
5535 SELECT CASE ( iclo )
5539 WRITE(0,
'(/1A,1A,1I2/)')
'W3DIV2 ERROR -- ', &
5540 'unsupported ICLO: ',iclo
5542 IF (
PRESENT(rc) )
THEN
5550 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
5551 WRITE(0,
'(/1A,1A/)')
'W3DIV2 ERROR -- ', &
5552 'tripole grid closure requires NP even'
5554 IF (
PRESENT(rc) )
THEN
5565 CALL get_fdw3 ( nfd, m, k, c )
5570 DO i2 = lbo(2), ubo(2)
5571 DO i1 = lbo(1), ubo(1)
5572 IF (
PRESENT(mask) )
THEN
5573 IF ( mask(i1,i2) ) cycle
5582 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
5583 prange, qrange, lbi, ubi, p, q, &
5584 f4=sxx, dfdp=dxxdp, dfdq=dxxdq, &
5585 g4=syy, dgdp=dyydp, dgdq=dyydq, &
5586 h4=sxy, dhdp=dxydp, dhdq=dxydq, &
5587 mask=mask, rc=istat )
5588 IF ( istat .NE. 0 )
THEN
5589 IF (
PRESENT(rc) )
THEN
5596 dxxdx = dxxdp*dpdx(i1,i2) + dxxdq*dqdx(i1,i2)
5597 dyydy = dyydp*dpdy(i1,i2) + dyydq*dqdy(i1,i2)
5598 dxydx = dxydp*dpdx(i1,i2) + dxydq*dqdx(i1,i2)
5599 dxydy = dxydp*dpdy(i1,i2) + dxydq*dqdy(i1,i2)
5600 dsx(i1,i2) = dxxdx + dxydy
5601 dsy(i1,i2) = dxydx + dyydy
5605 END SUBROUTINE w3div2_r4
5609 SUBROUTINE w3div2_r8( NFD, IJG, ICLO, PTILED, QTILED, &
5610 PRANGE, QRANGE, LBI, UBI, LBO, UBO, &
5611 DPDX, DPDY, DQDX, DQDY, &
5612 SXX, SYY, SXY, DSX, DSY, MASK, RC )
5614 INTEGER,
INTENT(IN) :: nfd
5615 LOGICAL,
INTENT(IN) :: ijg
5616 INTEGER,
INTENT(IN) :: iclo
5617 LOGICAL,
INTENT(IN) :: ptiled, qtiled
5618 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
5619 INTEGER,
INTENT(IN) :: lbi(2), ubi(2)
5620 INTEGER,
INTENT(IN) :: lbo(2), ubo(2)
5621 REAL(8),
INTENT(IN) :: dpdx(lbi(1):ubi(1),lbi(2):ubi(2))
5622 REAL(8),
INTENT(IN) :: dpdy(lbi(1):ubi(1),lbi(2):ubi(2))
5623 REAL(8),
INTENT(IN) :: dqdx(lbi(1):ubi(1),lbi(2):ubi(2))
5624 REAL(8),
INTENT(IN) :: dqdy(lbi(1):ubi(1),lbi(2):ubi(2))
5625 REAL(8),
INTENT(IN) :: sxx(lbi(1):ubi(1),lbi(2):ubi(2))
5626 REAL(8),
INTENT(IN) :: syy(lbi(1):ubi(1),lbi(2):ubi(2))
5627 REAL(8),
INTENT(IN) :: sxy(lbi(1):ubi(1),lbi(2):ubi(2))
5628 REAL(8),
INTENT(OUT) :: dsx(lbo(1):ubo(1),lbo(2):ubo(2))
5629 REAL(8),
INTENT(OUT) :: dsy(lbo(1):ubo(1),lbo(2):ubo(2))
5630 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lbi(1):ubi(1),lbi(2):ubi(2))
5631 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
5634 INTEGER,
PARAMETER :: m = 1
5635 INTEGER :: np, nq, i1, i2, p, q
5637 INTEGER :: k(0:nfd,0:nfd,1:nfd)
5638 REAL(8) :: c(0:nfd,0:nfd,1:nfd)
5639 REAL(8) :: dxxdp, dxxdq, dyydp, dyydq, dxydp, dxydq
5640 REAL(8) :: dxxdx, dyydy, dxydx, dxydy
5642 INTEGER,
SAVE :: ient = 0
5643 CALL strace (ient,
'W3DIV2_R8')
5648 IF (
PRESENT(rc) ) rc = 0
5650 IF ( nfd.LE.0 .OR. mod(nfd,2).NE.0 )
THEN
5651 WRITE(0,
'(/1A,1A/)')
'W3DIV2 ERROR -- ', &
5652 'NFD must be even and greater than zero'
5654 IF (
PRESENT(rc) )
THEN
5662 np = prange(2) - prange(1) + 1
5663 nq = qrange(2) - qrange(1) + 1
5665 SELECT CASE ( iclo )
5669 WRITE(0,
'(/1A,1A,1I2/)')
'W3DIV2 ERROR -- ', &
5670 'unsupported ICLO: ',iclo
5672 IF (
PRESENT(rc) )
THEN
5680 IF ( iclo.EQ.
iclo_trpl .AND. mod(np,2).NE.0 )
THEN
5681 WRITE(0,
'(/1A,1A/)')
'W3DIV2 ERROR -- ', &
5682 'tripole grid closure requires NP even'
5684 IF (
PRESENT(rc) )
THEN
5695 CALL get_fdw3 ( nfd, m, k, c )
5700 DO i2 = lbo(2), ubo(2)
5701 DO i1 = lbo(1), ubo(1)
5702 IF (
PRESENT(mask) )
THEN
5703 IF ( mask(i1,i2) ) cycle
5712 CALL dfdpq ( nfd, k, c, ijg, iclo, ptiled, qtiled, &
5713 prange, qrange, lbi, ubi, p, q, &
5714 f8=sxx, dfdp=dxxdp, dfdq=dxxdq, &
5715 g8=syy, dgdp=dyydp, dgdq=dyydq, &
5716 h8=sxy, dhdp=dxydp, dhdq=dxydq, &
5717 mask=mask, rc=istat )
5718 IF ( istat .NE. 0 )
THEN
5719 IF (
PRESENT(rc) )
THEN
5726 dxxdx = dxxdp*dpdx(i1,i2) + dxxdq*dqdx(i1,i2)
5727 dyydy = dyydp*dpdy(i1,i2) + dyydq*dqdy(i1,i2)
5728 dxydx = dxydp*dpdx(i1,i2) + dxydq*dqdx(i1,i2)
5729 dxydy = dxydp*dpdy(i1,i2) + dxydq*dqdy(i1,i2)
5730 dsx(i1,i2) = dxxdx + dxydy
5731 dsy(i1,i2) = dxydx + dyydy
5735 END SUBROUTINE w3div2_r8
5802 FUNCTION w3dist_r4( LLG, XT, YT, XS, YS )
RESULT(DIST)
5805 LOGICAL,
INTENT(IN) :: llg
5806 REAL(4),
INTENT(IN) :: xt, yt
5807 REAL(4),
INTENT(IN) :: xs, ys
5810 REAL(8) :: xt8, yt8, xs8, ys8
5812 INTEGER,
SAVE :: ient = 0
5813 CALL strace (ient,
'W3DIST_R4')
5821 dist = w3dist( llg, xt8, yt8, xs8, ys8 )
5823 END FUNCTION w3dist_r4
5827 #define DIST_WITH_SINE
5828 #define DIST_CHECK_NAN____disabled
5829 FUNCTION w3dist_r8( LLG, XT, YT, XS, YS )
RESULT(DIST)
5832 LOGICAL,
INTENT(IN) :: llg
5833 REAL(8),
INTENT(IN) :: xt, yt
5834 REAL(8),
INTENT(IN) :: xs, ys
5837 REAL(8) :: dx, dy, slam, sphi, argd
5839 INTEGER,
SAVE :: ient = 0
5840 CALL strace (ient,
'W3DIST_R8')
5849 IF ( abs(dx) .GT. d270 )
THEN
5850 dx = dx - sign(d360,dx)
5852 #ifdef DIST_WITH_SINE
5855 slam = sin(half*dx*d2r)
5856 sphi = sin(half*dy*d2r)
5857 argd = sqrt( cos(yt*d2r)*cos(ys*d2r)*slam*slam + sphi*sphi )
5858 dist = r2d*two*asin( argd )
5862 argd = min( one, cos(yt*d2r)*cos(ys*d2r)*cos(dx*d2r) &
5863 + sin(yt*d2r)*sin(ys*d2r) )
5864 dist = r2d*acos( argd )
5868 dist = sqrt( dx**2 + dy**2 )
5870 #ifdef DIST_CHECK_NAN
5871 IF ( w3inan(dist) )
THEN
5872 WRITE(0,
'(/1A/)')
'W3DIST_R8 ERROR -- result is NaN'
5877 END FUNCTION w3dist_r8
5943 SUBROUTINE w3splx_0d_r4( LAM0, PHI0, C0, LAM, PHI, X, Y )
5945 REAL(4),
INTENT(IN) :: lam0, phi0, c0
5946 REAL(4),
INTENT(IN) :: lam, phi
5947 REAL(4),
INTENT(OUT):: x, y
5950 REAL(8) :: k, k0, clam, slam, cphi0, cphi, sphi0, sphi
5952 INTEGER,
SAVE :: ient = 0
5953 CALL strace (ient,
'W3SPLX_0D_R4')
5956 clam = cos((lam-lam0)*d2r)
5957 slam = sin((lam-lam0)*d2r)
5958 cphi0 = cos(phi0*d2r)
5960 sphi0 = sin(phi0*d2r)
5962 k0 = cos(half*c0*d2r)**2
5963 k = two*k0*rearth/(one+sphi0*sphi+cphi0*cphi*clam)
5965 y = k*(cphi0*sphi-sphi0*cphi*clam)
5967 END SUBROUTINE w3splx_0d_r4
5971 SUBROUTINE w3splx_0d_r8( LAM0, PHI0, C0, LAM, PHI, X, Y )
5973 REAL(8),
INTENT(IN) :: lam0, phi0, c0
5974 REAL(8),
INTENT(IN) :: lam, phi
5975 REAL(8),
INTENT(OUT):: x, y
5978 REAL(8) :: k, k0, clam, slam, cphi0, cphi, sphi0, sphi
5980 INTEGER,
SAVE :: ient = 0
5981 CALL strace (ient,
'W3SPLX_0D_R8')
5984 clam = cos((lam-lam0)*d2r)
5985 slam = sin((lam-lam0)*d2r)
5986 cphi0 = cos(phi0*d2r)
5988 sphi0 = sin(phi0*d2r)
5990 k0 = cos(half*c0*d2r)**2
5991 k = two*k0*rearth/(one+sphi0*sphi+cphi0*cphi*clam)
5993 y = k*(cphi0*sphi-sphi0*cphi*clam)
5995 END SUBROUTINE w3splx_0d_r8
5999 SUBROUTINE w3splx_1d_r4( LAM0, PHI0, C0, LAM, PHI, X, Y )
6001 REAL(4),
INTENT(IN) :: lam0, phi0, c0
6002 REAL(4),
INTENT(IN) :: lam(:), phi(:)
6003 REAL(4),
INTENT(OUT):: x(:), y(:)
6008 INTEGER,
SAVE :: ient = 0
6009 CALL strace (ient,
'W3SPLX_1D_R4')
6012 DO i = lbound(lam,1),ubound(lam,1)
6013 CALL w3splx( lam0, phi0, c0, lam(i), phi(i), x(i), y(i) )
6016 END SUBROUTINE w3splx_1d_r4
6020 SUBROUTINE w3splx_1d_r8( LAM0, PHI0, C0, LAM, PHI, X, Y )
6022 REAL(8),
INTENT(IN) :: lam0, phi0, c0
6023 REAL(8),
INTENT(IN) :: lam(:), phi(:)
6024 REAL(8),
INTENT(OUT):: x(:), y(:)
6029 INTEGER,
SAVE :: ient = 0
6030 CALL strace (ient,
'W3SPLX_1D_R8')
6033 DO i = lbound(lam,1),ubound(lam,1)
6034 CALL w3splx( lam0, phi0, c0, lam(i), phi(i), x(i), y(i) )
6037 END SUBROUTINE w3splx_1d_r8
6041 SUBROUTINE w3splx_2d_r4( LAM0, PHI0, C0, LAM, PHI, X, Y )
6043 REAL(4),
INTENT(IN) :: lam0, phi0, c0
6044 REAL(4),
INTENT(IN) :: lam(:,:), phi(:,:)
6045 REAL(4),
INTENT(OUT):: x(:,:), y(:,:)
6050 INTEGER,
SAVE :: ient = 0
6051 CALL strace (ient,
'W3SPLX_2D_R4')
6054 DO j = lbound(lam,2),ubound(lam,2)
6055 DO i = lbound(lam,1),ubound(lam,1)
6056 CALL w3splx( lam0, phi0, c0, lam(i,j), phi(i,j), x(i,j), y(i,j) )
6060 END SUBROUTINE w3splx_2d_r4
6064 SUBROUTINE w3splx_2d_r8( LAM0, PHI0, C0, LAM, PHI, X, Y )
6066 REAL(8),
INTENT(IN) :: lam0, phi0, c0
6067 REAL(8),
INTENT(IN) :: lam(:,:), phi(:,:)
6068 REAL(8),
INTENT(OUT):: x(:,:), y(:,:)
6073 INTEGER,
SAVE :: ient = 0
6074 CALL strace (ient,
'W3SPLX_2D_R8')
6077 DO j = lbound(lam,2),ubound(lam,2)
6078 DO i = lbound(lam,1),ubound(lam,1)
6079 CALL w3splx( lam0, phi0, c0, lam(i,j), phi(i,j), x(i,j), y(i,j) )
6083 END SUBROUTINE w3splx_2d_r8
6149 SUBROUTINE w3spxl_0d_r4( LAM0, PHI0, C0, X, Y, LAM, PHI )
6151 REAL(4),
INTENT(IN) :: lam0, phi0, c0
6152 REAL(4),
INTENT(IN) :: x, y
6153 REAL(4),
INTENT(OUT):: lam, phi
6156 REAL(8) :: k0, rho, c, cosc, sinc, cphi0, sphi0
6158 INTEGER,
SAVE :: ient = 0
6159 CALL strace (ient,
'W3SPXL_0D_R4')
6162 k0 = cos(half*c0*d2r)**2
6164 c = two*atan2(rho,two*rearth*k0)
6167 cphi0 = cos(phi0*d2r)
6168 sphi0 = sin(phi0*d2r)
6169 phi = asin(cosc*sphi0+y*sinc*cphi0/rho)*r2d
6170 lam = lam0 + atan2(x*sinc,rho*cphi0*cosc-y*sphi0*sinc)*r2d
6172 END SUBROUTINE w3spxl_0d_r4
6176 SUBROUTINE w3spxl_0d_r8( LAM0, PHI0, C0, X, Y, LAM, PHI )
6178 REAL(8),
INTENT(IN) :: lam0, phi0, c0
6179 REAL(8),
INTENT(IN) :: x, y
6180 REAL(8),
INTENT(OUT):: lam, phi
6183 REAL(8) :: k0, rho, c, cosc, sinc, cphi0, sphi0
6185 INTEGER,
SAVE :: ient = 0
6186 CALL strace (ient,
'W3SPXL_0D_R8')
6189 k0 = cos(half*c0*d2r)**2
6191 c = two*atan2(rho,two*rearth*k0)
6194 cphi0 = cos(phi0*d2r)
6195 sphi0 = sin(phi0*d2r)
6196 phi = asin(cosc*sphi0+y*sinc*cphi0/rho)*r2d
6197 lam = lam0 + atan2(x*sinc,rho*cphi0*cosc-y*sphi0*sinc)*r2d
6199 END SUBROUTINE w3spxl_0d_r8
6203 SUBROUTINE w3spxl_1d_r4( LAM0, PHI0, C0, X, Y, LAM, PHI )
6205 REAL(4),
INTENT(IN) :: lam0, phi0, c0
6206 REAL(4),
INTENT(IN) :: x(:), y(:)
6207 REAL(4),
INTENT(OUT):: lam(:), phi(:)
6212 INTEGER,
SAVE :: ient = 0
6213 CALL strace (ient,
'W3SPXL_1D_R4')
6216 DO i = lbound(x,1),ubound(x,1)
6217 CALL w3spxl( lam0, phi0, c0, x(i), y(i), lam(i), phi(i) )
6220 END SUBROUTINE w3spxl_1d_r4
6224 SUBROUTINE w3spxl_1d_r8( LAM0, PHI0, C0, X, Y, LAM, PHI )
6226 REAL(8),
INTENT(IN) :: lam0, phi0, c0
6227 REAL(8),
INTENT(IN) :: x(:), y(:)
6228 REAL(8),
INTENT(OUT):: lam(:), phi(:)
6233 INTEGER,
SAVE :: ient = 0
6234 CALL strace (ient,
'W3SPXL_1D_R8')
6237 DO i = lbound(x,1),ubound(x,1)
6238 CALL w3spxl( lam0, phi0, c0, x(i), y(i), lam(i), phi(i) )
6241 END SUBROUTINE w3spxl_1d_r8
6245 SUBROUTINE w3spxl_2d_r4( LAM0, PHI0, C0, X, Y, LAM, PHI )
6247 REAL(4),
INTENT(IN) :: lam0, phi0, c0
6248 REAL(4),
INTENT(IN) :: x(:,:), y(:,:)
6249 REAL(4),
INTENT(OUT):: lam(:,:), phi(:,:)
6254 INTEGER,
SAVE :: ient = 0
6255 CALL strace (ient,
'W3SPXL_2D_R4')
6258 DO j = lbound(x,2),ubound(x,2)
6259 DO i = lbound(x,1),ubound(x,1)
6260 CALL w3spxl( lam0, phi0, c0, x(i,j), y(i,j), lam(i,j), phi(i,j) )
6264 END SUBROUTINE w3spxl_2d_r4
6268 SUBROUTINE w3spxl_2d_r8( LAM0, PHI0, C0, X, Y, LAM, PHI )
6270 REAL(8),
INTENT(IN) :: lam0, phi0, c0
6271 REAL(8),
INTENT(IN) :: x(:,:), y(:,:)
6272 REAL(8),
INTENT(OUT):: lam(:,:), phi(:,:)
6277 INTEGER,
SAVE :: ient = 0
6278 CALL strace (ient,
'W3SPXL_2D_R8')
6281 DO j = lbound(x,2),ubound(x,2)
6282 DO i = lbound(x,1),ubound(x,1)
6283 CALL w3spxl( lam0, phi0, c0, x(i,j), y(i,j), lam(i,j), phi(i,j) )
6287 END SUBROUTINE w3spxl_2d_r8
6350 SUBROUTINE w3trll_0d_r4( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6352 REAL(4),
INTENT(IN) :: lam0, phi0
6353 REAL(4),
INTENT(IN) :: lam1, phi1
6354 REAL(4),
INTENT(OUT):: lam, phi
6357 REAL(8) :: clam, slam, calp, salp, cphi, sphi
6359 INTEGER,
SAVE :: ient = 0
6360 CALL strace (ient,
'W3TRLL_0D_R4')
6363 clam = cos((lam1-lam0)*d2r)
6364 slam = sin((lam1-lam0)*d2r)
6365 calp = cos(phi0*d2r)
6366 salp = sin(phi0*d2r)
6367 cphi = cos(phi1*d2r)
6368 sphi = sin(phi1*d2r)
6369 lam = lam0 + atan2(cphi*slam,salp*cphi*clam+calp*sphi)*r2d
6370 phi = asin(salp*sphi-calp*cphi*clam)*r2d
6372 END SUBROUTINE w3trll_0d_r4
6376 SUBROUTINE w3trll_0d_r8( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6378 REAL(8),
INTENT(IN) :: lam0, phi0
6379 REAL(8),
INTENT(IN) :: lam1, phi1
6380 REAL(8),
INTENT(OUT):: lam, phi
6383 REAL(8) :: clam, slam, calp, salp, cphi, sphi
6385 INTEGER,
SAVE :: ient = 0
6386 CALL strace (ient,
'W3TRLL_0D_R8')
6389 clam = cos((lam1-lam0)*d2r)
6390 slam = sin((lam1-lam0)*d2r)
6391 calp = cos(phi0*d2r)
6392 salp = sin(phi0*d2r)
6393 cphi = cos(phi1*d2r)
6394 sphi = sin(phi1*d2r)
6395 lam = lam0 + atan2(cphi*slam,salp*cphi*clam+calp*sphi)*r2d
6396 phi = asin(salp*sphi-calp*cphi*clam)*r2d
6398 END SUBROUTINE w3trll_0d_r8
6402 SUBROUTINE w3trll_1d_r4( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6404 REAL(4),
INTENT(IN) :: lam0, phi0
6405 REAL(4),
INTENT(IN) :: lam1(:), phi1(:)
6406 REAL(4),
INTENT(OUT):: lam(:), phi(:)
6411 INTEGER,
SAVE :: ient = 0
6412 CALL strace (ient,
'W3TRLL_1D_R4')
6415 DO i = lbound(lam1,1),ubound(lam1,1)
6416 CALL w3trll( lam0, phi0, lam1(i), phi1(i), lam(i), phi(i) )
6419 END SUBROUTINE w3trll_1d_r4
6423 SUBROUTINE w3trll_1d_r8( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6425 REAL(8),
INTENT(IN) :: lam0, phi0
6426 REAL(8),
INTENT(IN) :: lam1(:), phi1(:)
6427 REAL(8),
INTENT(OUT):: lam(:), phi(:)
6432 INTEGER,
SAVE :: ient = 0
6433 CALL strace (ient,
'W3TRLL_1D_R8')
6436 DO i = lbound(lam1,1),ubound(lam1,1)
6437 CALL w3trll( lam0, phi0, lam1(i), phi1(i), lam(i), phi(i) )
6440 END SUBROUTINE w3trll_1d_r8
6444 SUBROUTINE w3trll_2d_r4( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6446 REAL(4),
INTENT(IN) :: lam0, phi0
6447 REAL(4),
INTENT(IN) :: lam1(:,:), phi1(:,:)
6448 REAL(4),
INTENT(OUT):: lam(:,:), phi(:,:)
6453 INTEGER,
SAVE :: ient = 0
6454 CALL strace (ient,
'W3TRLL_2D_R4')
6457 DO j = lbound(lam1,2),ubound(lam1,2)
6458 DO i = lbound(lam1,1),ubound(lam1,1)
6459 CALL w3trll( lam0, phi0, lam1(i,j), phi1(i,j), lam(i,j), phi(i,j) )
6463 END SUBROUTINE w3trll_2d_r4
6467 SUBROUTINE w3trll_2d_r8( LAM0, PHI0, LAM1, PHI1, LAM, PHI )
6469 REAL(8),
INTENT(IN) :: lam0, phi0
6470 REAL(8),
INTENT(IN) :: lam1(:,:), phi1(:,:)
6471 REAL(8),
INTENT(OUT):: lam(:,:), phi(:,:)
6476 INTEGER,
SAVE :: ient = 0
6477 CALL strace (ient,
'W3TRLL_2D_R8')
6480 DO j = lbound(lam1,2),ubound(lam1,2)
6481 DO i = lbound(lam1,1),ubound(lam1,1)
6482 CALL w3trll( lam0, phi0, lam1(i,j), phi1(i,j), lam(i,j), phi(i,j) )
6486 END SUBROUTINE w3trll_2d_r8
6551 FUNCTION w3llaz_r4( LAM1, PHI1, LAM2, PHI2 )
RESULT(AZ)
6554 REAL(4),
INTENT(IN):: lam1, phi1
6555 REAL(4),
INTENT(IN):: lam2, phi2
6558 REAL(8) :: clam, slam, cph1, sph1, cph2, sph2
6560 INTEGER,
SAVE :: ient = 0
6561 CALL strace (ient,
'W3LLAZ_R4')
6564 clam = cos((lam2-lam1)*d2r)
6565 slam = sin((lam2-lam1)*d2r)
6566 cph1 = cos(phi1*d2r)
6567 sph1 = sin(phi1*d2r)
6568 cph2 = cos(phi2*d2r)
6569 sph2 = sin(phi2*d2r)
6570 az = atan2(cph2*slam,cph1*sph2-sph1*cph2*clam)*r2d
6572 END FUNCTION w3llaz_r4
6576 FUNCTION w3llaz_r8( LAM1, PHI1, LAM2, PHI2 )
RESULT(AZ)
6579 REAL(8),
INTENT(IN):: lam1, phi1
6580 REAL(8),
INTENT(IN):: lam2, phi2
6583 REAL(8) :: clam, slam, cph1, sph1, cph2, sph2
6585 INTEGER,
SAVE :: ient = 0
6586 CALL strace (ient,
'W3LLAZ_R8')
6589 clam = cos((lam2-lam1)*d2r)
6590 slam = sin((lam2-lam1)*d2r)
6591 cph1 = cos(phi1*d2r)
6592 sph1 = sin(phi1*d2r)
6593 cph2 = cos(phi2*d2r)
6594 sph2 = sin(phi2*d2r)
6595 az = atan2(cph2*slam,cph1*sph2-sph1*cph2*clam)*r2d
6597 END FUNCTION w3llaz_r8
6661 SUBROUTINE w3fdwt_r4 ( N, ND, M, Z, X, C )
6663 INTEGER,
INTENT(IN) :: n, nd, m
6664 REAL(4),
INTENT(IN) :: z
6665 REAL(4),
INTENT(IN) :: x(0:nd)
6666 REAL(4),
INTENT(OUT) :: c(0:nd,0:m)
6669 INTEGER :: i, j, k, mn
6670 REAL(8) :: c1, c2, c3, c4, c5
6672 INTEGER,
SAVE :: ient = 0
6673 CALL strace (ient,
'W3FDWT_R4')
6688 IF ( j.EQ.i-1 )
THEN
6690 c(i,k) = c1*(k*c(i-1,k-1)-c5*c(i-1,k))/c2
6692 c(i,0) = -c1*c5*c(i-1,0)/c2
6695 c(j,k) = (c4*c(j,k)-k*c(j,k-1))/c3
6697 c(j,0) = c4*c(j,0)/c3
6702 END SUBROUTINE w3fdwt_r4
6706 SUBROUTINE w3fdwt_r8 ( N, ND, M, Z, X, C )
6708 INTEGER,
INTENT(IN) :: n, nd, m
6709 REAL(8),
INTENT(IN) :: z
6710 REAL(8),
INTENT(IN) :: x(0:nd)
6711 REAL(8),
INTENT(OUT) :: c(0:nd,0:m)
6714 INTEGER :: i, j, k, mn
6715 REAL(8) :: c1, c2, c3, c4, c5
6717 INTEGER,
SAVE :: ient = 0
6718 CALL strace (ient,
'W3FDWT_R4')
6733 IF ( j.EQ.i-1 )
THEN
6735 c(i,k) = c1*(k*c(i-1,k-1)-c5*c(i-1,k))/c2
6737 c(i,0) = -c1*c5*c(i-1,0)/c2
6740 c(j,k) = (c4*c(j,k)-k*c(j,k-1))/c3
6742 c(j,0) = c4*c(j,0)/c3
6747 END SUBROUTINE w3fdwt_r8
6819 FUNCTION w3nnsc( NLVL )
RESULT(NNS)
6821 INTEGER,
INTENT(IN) :: nlvl
6824 INTEGER :: i, j, l, n
6826 INTEGER,
SAVE :: ient = 0
6827 CALL strace (ient,
'W3NNSC')
6835 nns%NNBR = (2*nlvl+1)**2
6838 ALLOCATE(nns%N1(0:nns%NLVL))
6839 ALLOCATE(nns%N2(0:nns%NLVL))
6840 ALLOCATE(nns%DI(0:nns%NNBR-1))
6841 ALLOCATE(nns%DJ(0:nns%NNBR-1))
6847 nns%N1(l) = 0; nns%N2(l) = (2*l+1)**2-1;
6848 nns%DI(n) = 0; nns%DJ(n) = 0;
6852 nns%N1(l) = (2*l-1)**2; nns%N2(l) = (2*l+1)**2-1;
6857 nns%DI(n) = i; nns%DJ(n) = j;
6863 nns%DI(n) = i; nns%DJ(n) = j;
6869 nns%DI(n) = i; nns%DJ(n) = j;
6875 nns%DI(n) = i; nns%DJ(n) = j;
6935 INTEGER,
SAVE :: ient = 0
6936 CALL strace (ient,
'W3NNSD')
6939 IF (
ASSOCIATED(nns) )
THEN
6942 IF (
ASSOCIATED(nns%N1) )
THEN
6943 DEALLOCATE(nns%N1);
NULLIFY(nns%N1);
6945 IF (
ASSOCIATED(nns%N2) )
THEN
6946 DEALLOCATE(nns%N2);
NULLIFY(nns%N2);
6948 IF (
ASSOCIATED(nns%DI) )
THEN
6949 DEALLOCATE(nns%DI);
NULLIFY(nns%DI);
6951 IF (
ASSOCIATED(nns%DJ) )
THEN
6952 DEALLOCATE(nns%DJ);
NULLIFY(nns%DJ);
7011 SUBROUTINE w3nnsp(NNS, IUNIT)
7013 INTEGER,
OPTIONAL,
INTENT(IN) :: iunit
7016 INTEGER :: ndst, l, n
7018 INTEGER,
SAVE :: ient = 0
7019 CALL strace (ient,
'W3NNSP')
7022 IF (
PRESENT(iunit) )
THEN
7028 WRITE(ndst,
'(A,2I6)')
'nlvl,nnbr:',nns%NLVL,nns%NNBR
7030 DO n=nns%N1(l),nns%N2(l)
7031 WRITE(ndst,
'(A,4I6)')
'l,n,di,dj:',l,n,nns%DI(n),nns%DJ(n)
7086 SUBROUTINE w3sort_r4( N, I, J, D )
7088 INTEGER,
INTENT(IN) :: n
7089 INTEGER,
INTENT(INOUT) :: i(n)
7090 INTEGER,
INTENT(INOUT) :: j(n)
7091 REAL(4),
INTENT(INOUT) :: d(n)
7094 INTEGER :: k, l, im, jm
7097 INTEGER,
SAVE :: ient = 0
7098 CALL strace (ient,
'W3SORT_R4')
7103 IF ( d(l) .LT. d(k) )
THEN
7104 im = i(k); jm = j(k); dm = d(k);
7105 i(k) = i(l); j(k) = j(l); d(k) = d(l);
7106 i(l) = im; j(l) = jm; d(l) = dm;
7111 END SUBROUTINE w3sort_r4
7115 SUBROUTINE w3sort_r8( N, I, J, D )
7117 INTEGER,
INTENT(IN) :: n
7118 INTEGER,
INTENT(INOUT) :: i(n)
7119 INTEGER,
INTENT(INOUT) :: j(n)
7120 REAL(8),
INTENT(INOUT) :: d(n)
7123 INTEGER :: k, l, im, jm
7126 INTEGER,
SAVE :: ient = 0
7127 CALL strace (ient,
'W3SORT_R8')
7132 IF ( d(l) .LT. d(k) )
THEN
7133 im = i(k); jm = j(k); dm = d(k);
7134 i(k) = i(l); j(k) = j(l); d(k) = d(l);
7135 i(l) = im; j(l) = jm; d(l) = dm;
7140 END SUBROUTINE w3sort_r8
7191 SUBROUTINE w3isrt_r4( II, JJ, DD, N, I, J, D )
7193 INTEGER,
INTENT(IN) :: ii
7194 INTEGER,
INTENT(IN) :: jj
7195 REAL(4),
INTENT(IN) :: dd
7196 INTEGER,
INTENT(IN) :: n
7197 INTEGER,
INTENT(INOUT) :: i(n)
7198 INTEGER,
INTENT(INOUT) :: j(n)
7199 REAL(4),
INTENT(INOUT) :: d(n)
7204 INTEGER,
SAVE :: ient = 0
7205 CALL strace (ient,
'W3ISRT_R4')
7209 IF ( dd .LT. d(k) )
THEN
7212 i(l) = i(l-1); j(l) = j(l-1); d(l) = d(l-1);
7215 i(k) = ii; j(k) = jj; d(k) = dd;
7220 END SUBROUTINE w3isrt_r4
7224 SUBROUTINE w3isrt_r8( II, JJ, DD, N, I, J, D )
7226 INTEGER,
INTENT(IN) :: ii
7227 INTEGER,
INTENT(IN) :: jj
7228 REAL(8),
INTENT(IN) :: dd
7229 INTEGER,
INTENT(IN) :: n
7230 INTEGER,
INTENT(INOUT) :: i(n)
7231 INTEGER,
INTENT(INOUT) :: j(n)
7232 REAL(8),
INTENT(INOUT) :: d(n)
7237 INTEGER,
SAVE :: ient = 0
7238 CALL strace (ient,
'W3ISRT_R8')
7242 IF ( dd .LT. d(k) )
THEN
7245 i(l) = i(l-1); j(l) = j(l-1); d(l) = d(l-1);
7248 i(k) = ii; j(k) = jj; d(k) = dd;
7253 END SUBROUTINE w3isrt_r8
7304 FUNCTION w3inan_r4( X )
RESULT(INAN)
7307 REAL(4),
INTENT(IN) :: x
7311 INTEGER,
SAVE :: ient = 0
7312 CALL strace (ient,
'W3INAN_R4')
7316 inan = .NOT. ( x .GE. -huge(x) .AND. x .LE. huge(x) )
7318 END FUNCTION w3inan_r4
7322 FUNCTION w3inan_r8( X )
RESULT(INAN)
7325 REAL(8),
INTENT(IN) :: x
7329 INTEGER,
SAVE :: ient = 0
7330 CALL strace (ient,
'W3INAN_R8')
7334 inan = .NOT. ( x .GE. -huge(x) .AND. x .LE. huge(x) )
7336 END FUNCTION w3inan_r8
7354 FUNCTION gsu_create( IJG, LLG, ICLO, LB, UB, XG4, YG4, XG8, YG8, &
7355 BBOX_ONLY, NCB, NNP, DEBUG )
RESULT(GSU)
7358 LOGICAL,
INTENT(IN) :: ijg
7359 LOGICAL,
INTENT(IN) :: llg
7360 INTEGER,
INTENT(IN) :: iclo
7361 INTEGER,
INTENT(IN) :: lb(2)
7362 INTEGER,
INTENT(IN) :: ub(2)
7363 REAL(4),
TARGET,
OPTIONAL :: xg4(lb(1):ub(1),lb(2):ub(2))
7364 REAL(4),
TARGET,
OPTIONAL :: yg4(lb(1):ub(1),lb(2):ub(2))
7365 REAL(8),
TARGET,
OPTIONAL :: xg8(lb(1):ub(1),lb(2):ub(2))
7366 REAL(8),
TARGET,
OPTIONAL :: yg8(lb(1):ub(1),lb(2):ub(2))
7367 LOGICAL,
INTENT(IN),
OPTIONAL :: bbox_only
7368 INTEGER,
INTENT(IN),
OPTIONAL :: ncb
7369 INTEGER,
INTENT(IN),
OPTIONAL :: nnp
7370 LOGICAL,
INTENT(IN),
OPTIONAL :: debug
7373 TYPE(class_gsu),
POINTER :: ptr
7374 LOGICAL :: type_r4, type_r8
7375 LOGICAL :: ldbg, lbbox, lbc, lpl, lnpl, lspl
7376 INTEGER :: lbx, lby, ubx, uby, nx, ny
7377 INTEGER :: lxc, lyc, uxc, uyc
7378 INTEGER :: i, j, k, l, n, ic(4), jc(4), ib, jb
7379 INTEGER :: ns, ib1(2), ib2(2), jb1(2), jb2(2), ibc(4), jbc(4)
7380 INTEGER :: istep, istat
7381 REAL(8) :: xc(4), yc(4)
7383 INTEGER,
SAVE :: ient = 0
7384 CALL strace (ient,
'W3GSUC')
7389 type_r4 =
PRESENT(xg4).AND.
PRESENT(yg4)
7390 type_r8 =
PRESENT(xg8).AND.
PRESENT(yg8)
7391 IF ( .NOT.type_r4.AND..NOT.type_r8 )
THEN
7392 WRITE(0,
'(/1A,1A,1I2/)')
'W3GSUC ERROR -- ', &
7393 'no input grid coordinates specified'
7411 SELECT CASE ( iclo )
7415 WRITE(0,
'(/1A,1A,1I2/)')
'W3GSUC ERROR -- ', &
7416 'unsupported ICLO: ',iclo
7420 IF ( iclo.EQ.
iclo_trpl .AND. mod(nx,2).NE.0 )
THEN
7421 WRITE(0,
'(/1A,1A/)')
'W3GSUC ERROR -- ', &
7422 'tripole grid closure requires NX=UBX-LBX+1 be even'
7426 IF (
PRESENT(bbox_only) )
THEN
7432 IF (
PRESENT(ncb) )
THEN
7433 IF ( ncb .LE. 0 )
THEN
7434 WRITE(0,
'(/1A,1A/)')
'W3GSUC ERROR -- ', &
7435 'NCB must be greater than zero'
7440 IF (
PRESENT(debug) )
THEN
7449 ALLOCATE(ptr, stat=istat)
7450 IF ( istat .NE. 0 )
THEN
7451 WRITE(0,
'(/1A,1A/)')
'W3GSUC ERROR -- ', &
7452 'gsu object allocation failed'
7480 IF ( .NOT.lbbox )
THEN
7481 IF (
PRESENT(nnp) )
THEN
7484 ptr%NNP =>
w3nnsc(nnp_default)
7492 lxc = lbx; lyc = lby;
7493 SELECT CASE ( iclo )
7495 uxc = ubx-1; uyc = uby-1;
7497 uxc = ubx; uyc = uby-1;
7499 uxc = ubx-1; uyc = uby;
7501 uxc = ubx; uyc = uby;
7503 uxc = ubx; uyc = uby;
7516 WRITE(*,
'(/A)')
'W3GSUC - check source grid'
7523 ic(1) = i ; jc(1) = j ;
7524 ic(2) = i+1; jc(2) = j ;
7525 ic(3) = i+1; jc(3) = j+1;
7526 ic(4) = i ; jc(4) = j+1;
7529 IF ( mod(iclo,2).EQ.0 ) &
7530 ic(l) = lbx + mod(nx - 1 + mod(ic(l) - lbx + 1, nx), nx)
7531 IF ( mod(iclo,3).EQ.0 ) &
7532 jc(l) = lby + mod(ny - 1 + mod(jc(l) - lby + 1, ny), ny)
7533 IF ( iclo.EQ.
iclo_trpl .AND. jc(l).GT.uby )
THEN
7534 ic(l) = ubx + lbx - ic(l)
7535 jc(l) = 2*uby - jc(l) + 1
7540 xc(l) = xg4(ic(l),jc(l))
7541 yc(l) = yg4(ic(l),jc(l))
7543 xc(l) = xg8(ic(l),jc(l))
7544 yc(l) = yg8(ic(l),jc(l))
7548 xc(l) = xg4(jc(l),ic(l))
7549 yc(l) = yg4(jc(l),ic(l))
7551 xc(l) = xg8(jc(l),ic(l))
7552 yc(l) = yg8(jc(l),ic(l))
7564 IF ( abs(xc(k)-xc(l)) .GT. d180 ) n = n + 1
7568 IF ( lbc .AND. ldbg ) &
7569 WRITE(*,
'(A,8I6)') &
7570 'W3GSUC -- cell includes branch cut:',ic(:),jc(:)
7573 lpl = n.EQ.1 .OR. count(abs(yc).EQ.d90).EQ.1
7574 IF ( lpl.AND.minval(yc).GT.zero )
THEN
7576 WRITE(*,
'(A,8I6)') &
7577 'W3GSUC -- cell includes N-pole:',ic(:),jc(:)
7580 IF ( lpl.AND.maxval(yc).LT.zero )
THEN
7582 WRITE(*,
'(A,8I6)') &
7583 'W3GSUC -- cell includes S-pole:',ic(:),jc(:)
7587 IF ( n.GT.0 ) ptr%LCLO = .true.
7597 ptr%XMIN = minval(xg4); ptr%XMAX = maxval(xg4);
7598 ptr%YMIN = minval(yg4); ptr%YMAX = maxval(yg4);
7600 ptr%XMIN = minval(xg8); ptr%XMAX = maxval(xg8);
7601 ptr%YMIN = minval(yg8); ptr%YMAX = maxval(yg8);
7603 IF ( ptr%LCLO )
THEN
7604 ptr%XMIN = zero; ptr%XMAX = d360;
7606 IF ( lspl ) ptr%YMIN = -d90
7607 IF ( lnpl ) ptr%YMAX = d90
7608 ptr%L360 = ptr%XMIN.GE.zero
7617 IF (
PRESENT(ncb) )
THEN
7618 ptr%NBX = max(1,nx/ncb)
7619 ptr%NBY = max(1,ny/ncb)
7621 ptr%NBX = max(1,nx/ncb_default)
7622 ptr%NBY = max(1,ny/ncb_default)
7624 ptr%DXB = (ptr%XMAX-ptr%XMIN)/real(ptr%NBX)
7625 ptr%DYB = (ptr%YMAX-ptr%YMIN)/real(ptr%NBY)
7629 WRITE(*,
'(/A,1I2,1L2,1I2)')
'W3GSUC - ICLO,LCLO,GKIND: ', &
7630 ptr%ICLO,ptr%LCLO,ptr%GKIND
7631 WRITE(*,
'(A,4E24.16)')
'W3GSUC - grid search domain:', &
7632 ptr%XMIN,ptr%YMIN,ptr%XMAX,ptr%YMAX
7633 WRITE(*,
'(A,2I6)')
'W3GSUC - number of search buckets:', &
7635 WRITE(*,
'(A,2E24.16)')
'W3GSUC - search bucket size:', &
7640 ALLOCATE(ptr%B(ptr%NBY,ptr%NBX),stat=istat)
7641 IF ( istat .NE. 0 )
THEN
7642 WRITE(0,
'(/1A,1A/)')
'W3GSUC ERROR -- ', &
7643 'search bucket array allocation failed'
7650 istep_loop:
DO istep=1,2
7653 IF ( istep .EQ. 2 )
THEN
7656 NULLIFY(ptr%B(jb,ib)%I)
7657 NULLIFY(ptr%B(jb,ib)%J)
7658 IF ( ptr%B(jb,ib)%N .GT. 0 )
THEN
7659 ALLOCATE(ptr%B(jb,ib)%I(ptr%B(jb,ib)%N),stat=istat)
7660 IF ( istat .NE. 0 )
THEN
7661 WRITE(0,
'(/1A,2A,3I6/)')
'W3GSUC ERROR -- ', &
7662 'search bucket cell-i list allocation failed -- ', &
7666 ALLOCATE(ptr%B(jb,ib)%J(ptr%B(jb,ib)%N),stat=istat)
7667 IF ( istat .NE. 0 )
THEN
7668 WRITE(0,
'(/1A,2A,3I6/)')
'W3GSUC ERROR -- ', &
7669 'search bucket cell-j list allocation failed -- ', &
7683 IF ( j.EQ.uyc .AND. i.GT.lbx+nx/2 ) cycle
7686 ic(1) = i ; jc(1) = j ;
7687 ic(2) = i+1; jc(2) = j ;
7688 ic(3) = i+1; jc(3) = j+1;
7689 ic(4) = i ; jc(4) = j+1;
7692 IF ( mod(iclo,2).EQ.0 ) &
7693 ic(l) = lbx + mod(nx - 1 + mod(ic(l) - lbx + 1, nx), nx)
7694 IF ( mod(iclo,3).EQ.0 ) &
7695 jc(l) = lby + mod(ny - 1 + mod(jc(l) - lby + 1, ny), ny)
7696 IF ( iclo.EQ.
iclo_trpl .AND. jc(l).GT.uby )
THEN
7697 ic(l) = ubx + lbx - ic(l)
7698 jc(l) = 2*uby - jc(l) + 1
7703 xc(l) = xg4(ic(l),jc(l))
7704 yc(l) = yg4(ic(l),jc(l))
7706 xc(l) = xg8(ic(l),jc(l))
7707 yc(l) = yg8(ic(l),jc(l))
7711 xc(l) = xg4(jc(l),ic(l))
7712 yc(l) = yg4(jc(l),ic(l))
7714 xc(l) = xg8(jc(l),ic(l))
7715 yc(l) = yg8(jc(l),ic(l))
7725 IF ( ptr%LCLO .OR. ptr%L360 )
THEN
7726 WHERE ( xc.LT.zero ) xc = xc + d360
7728 WHERE ( xc.GT.d180 ) xc = xc - d360
7734 IF ( abs(xc(k)-xc(l)) .GT. d180 ) n = n + 1
7740 lpl = n.EQ.1 .OR. count(abs(yc).EQ.d90).EQ.1
7744 ibc(l) = int((xc(l)-ptr%XMIN)/ptr%DXB)+1
7745 IF ( .NOT.ptr%LCLO ) ibc(l) = min(ibc(l),ptr%NBX)
7746 jbc(l) = min(int((yc(l)-ptr%YMIN)/ptr%DYB)+1,ptr%NBY)
7754 IF ( minval(yc).GT.zero )
THEN
7755 jb1(1) = max(1,minval(jbc))
7758 IF ( maxval(yc).LT.zero )
THEN
7760 jb2(1) = min(ptr%NBY,maxval(jbc))
7766 ELSE IF ( lbc )
THEN
7774 IF ( ibc(l) .GT. ptr%NBX/2 )
THEN
7775 ib1(1) = min(ib1(1),ibc(l))
7777 ib2(2) = max(ib2(2),ibc(l))
7780 jb1(:) = max(1,minval(jbc))
7781 jb2(:) = min(ptr%NBY,maxval(jbc))
7785 ib1(1) = max(1,minval(ibc))
7786 ib2(1) = min(ptr%NBX,maxval(ibc))
7787 jb1(1) = max(1,minval(jbc))
7788 jb2(1) = min(ptr%NBY,maxval(jbc))
7795 IF ( ldbg .AND. istep.EQ.1 )
THEN
7796 WRITE(*,
'(/A,2I6)')
'W3GSUC -- BUCKET SORT:',i,j
7797 WRITE(*,
'(A,2L6,1I6)')
'W3GSUC -- LBC,LPL:',lbc,lpl
7798 WRITE(*,
'(A,4I6)')
'W3GSUC -- IC:',ic(:)
7799 WRITE(*,
'(A,4I6)')
'W3GSUC -- JC:',jc(:)
7800 WRITE(*,
'(A,4E24.16)')
'W3GSUC -- XC:',xc(:)
7801 WRITE(*,
'(A,4E24.16)')
'W3GSUC -- YC:',yc(:)
7802 WRITE(*,
'(A,4I6)')
'W3GSUC -- IBC:',ibc(:)
7803 WRITE(*,
'(A,4I6)')
'W3GSUC -- JBC:',jbc(:)
7804 WRITE(*,
'(A,1I6)')
'W3GSUC -- NS:',ns
7805 WRITE(*,
'(A,4I6)')
'W3GSUC -- IB1:',ib1(:)
7806 WRITE(*,
'(A,4I6)')
'W3GSUC -- JB1:',jb1(:)
7807 WRITE(*,
'(A,4I6)')
'W3GSUC -- IB2:',ib2(:)
7808 WRITE(*,
'(A,4I6)')
'W3GSUC -- JB2:',jb2(:)
7814 ptr%B(jb,ib)%N = ptr%B(jb,ib)%N + 1
7815 IF ( istep .EQ. 2 )
THEN
7816 ptr%B(jb,ib)%I(ptr%B(jb,ib)%N) = ic(1)
7817 ptr%B(jb,ib)%J(ptr%B(jb,ib)%N) = jc(1)
7829 ptr%NNB =>
w3nnsc(nint(half*max(ptr%NBX,ptr%NBY)))
7833 WRITE(*,
'(/A,3I6,4E24.16)')
'W3GSUC - search bucket list:'
7834 WRITE(*,
'(3A6,4A14)')
'I',
'J',
'N',
'X1',
'Y1',
'X2',
'Y2'
7837 WRITE(*,
'(3I6,4E24.16)') ib,jb,ptr%B(jb,ib)%N, &
7838 ptr%XMIN+(ib-1)*ptr%DXB,ptr%YMIN+(jb-1)*ptr%DYB, &
7839 ptr%XMIN+(ib-0)*ptr%DXB,ptr%YMIN+(jb-0)*ptr%DYB
7849 END FUNCTION gsu_create
7853 SUBROUTINE getpqr( XT, YT, XS, YS, PR, QR, EPS, DEBUG )
7856 REAL(8),
INTENT(IN) :: xt
7857 REAL(8),
INTENT(IN) :: yt
7858 REAL(8),
INTENT(IN) :: xs(4)
7859 REAL(8),
INTENT(IN) :: ys(4)
7860 REAL(8),
INTENT(OUT) :: pr
7861 REAL(8),
INTENT(OUT) :: qr
7862 REAL(8),
INTENT(IN),
OPTIONAL :: eps
7863 LOGICAL,
INTENT(IN) ,
OPTIONAL :: debug
7866 INTEGER,
PARAMETER :: max_iter = 10
7867 REAL(8),
PARAMETER :: converge = 1d-6
7871 REAL(8) :: dxt, dx1, dx2, dx3, dxp, dyt, dy1, dy2, dy3, dyp
7872 REAL(8) :: mat1, mat2, mat3, mat4, delp, delq, det
7874 INTEGER,
SAVE :: ient = 0
7875 CALL strace (ient,
'GETPQR')
7878 IF (
PRESENT(eps) )
THEN
7879 IF ( eps .LT. zero )
THEN
7880 WRITE(0,
'(/2A/)')
'GETPQR ERROR -- ', &
7881 'EPS parameter must be >= 0'
7888 IF (
PRESENT(debug) )
THEN
7896 IF ( abs(xt-xs(k)).LE.leps .AND. abs(yt-ys(k)).LE.leps )
THEN
7899 pr = zero; qr = zero;
7901 pr = one; qr = zero;
7905 pr = zero; qr = one;
7908 WRITE(*,
'(A,I3,4E24.16)')
'GETPQR - COINCIDENT:', &
7909 k,abs(xt-xs(k)),abs(yt-ys(k)),pr,qr
7920 dy3 = ys(3) - ys(2) - dy2
7924 dx3 = xs(3) - xs(2) - dx2
7927 iter_loop:
DO iter=1,max_iter
7928 dyp = dyt - dy1*pr - dy2*qr - dy3*pr*qr
7929 dxp = dxt - dx1*pr - dx2*qr - dx3*pr*qr
7934 det = mat1*mat4 - mat2*mat3
7935 delp = (dyp*mat4 - mat2*dxp)/det
7936 delq = (mat1*dxp - dyp*mat3)/det
7938 WRITE(*,
'(A,I3,4E24.16)')
'GETPQR - ITER:', &
7939 iter,pr,qr,delp,delq
7942 IF ( abs(delp) < converge .AND. &
7943 abs(delq) < converge )
EXIT iter_loop
7947 IF ( iter .GT. max_iter )
THEN
7949 'GETPQR -- ERROR: exceeded max iteration count'
7950 WRITE(0,
'(A,2E24.16)')
'GETPQR - DEST POINT COORDS: ',xt,yt
7952 WRITE(0,
'(A,I1,A,2E24.16)') &
7953 'GETPQR - SRC POINT ',k,
': ',xs(k),ys(k)
7955 WRITE(0,
'(A,4E24.16)') &
7956 'GETPQR - CURRENT PR,QR,DELP,DELQ: ',pr,qr,delp,delq
7960 END SUBROUTINE getpqr
7964 SUBROUTINE getblc( GSU, I, J, PR, QR, LCMP, NS, LS, IS, JS, CS )
7970 TYPE(
t_gsu),
INTENT(IN) :: gsu
7971 INTEGER,
INTENT(IN) :: i, j
7972 REAL(8),
INTENT(IN) :: pr, qr
7973 LOGICAL,
INTENT(IN) :: lcmp
7974 INTEGER,
INTENT(OUT) :: ns
7975 LOGICAL,
POINTER,
INTENT(INOUT) :: ls(:)
7976 INTEGER,
POINTER,
INTENT(INOUT) :: is(:), js(:)
7977 REAL(8),
POINTER,
INTENT(INOUT) :: cs(:)
7980 LOGICAL :: ijg, llg, lclo
7981 INTEGER :: iclo, gkind
7982 INTEGER :: lbx, lby, ubx, uby, nx, ny
7986 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
7987 WRITE(0,
'(/2A/)')
'GETBLC ERROR -- ', &
7988 'grid search utility object not created'
7995 gkind = gsu%PTR%GKIND
7996 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
7997 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
7998 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
8001 IF (
ASSOCIATED(ls) )
THEN
8002 DEALLOCATE(ls);
NULLIFY(ls);
8004 IF (
ASSOCIATED(is) )
THEN
8005 DEALLOCATE(is);
NULLIFY(is);
8007 IF (
ASSOCIATED(js) )
THEN
8008 DEALLOCATE(js);
NULLIFY(js);
8010 IF (
ASSOCIATED(cs) )
THEN
8011 DEALLOCATE(cs);
NULLIFY(cs);
8016 ALLOCATE( ls(ns), is(ns), js(ns), cs(ns), stat=istat )
8017 IF ( istat .NE. 0 )
THEN
8018 WRITE(0,
'(/1A,1A/)')
'GETBLC ERROR -- ', &
8019 'array allocation failed'
8034 is(1) = i ; js(1) = j ;
8035 is(2) = i+1; js(2) = j ;
8036 is(3) = i+1; js(3) = j+1;
8037 is(4) = i ; js(4) = j+1;
8041 IF ( mod(iclo,2).EQ.0 ) &
8042 is(k) = lbx + mod(nx - 1 + mod(is(k) - lbx + 1, nx), nx)
8043 IF ( mod(iclo,3).EQ.0 ) &
8044 js(k) = lby + mod(ny - 1 + mod(js(k) - lby + 1, ny), ny)
8045 IF ( iclo.EQ.
iclo_trpl .AND. js(k).GT.uby )
THEN
8046 is(k) = ubx + lbx - is(k)
8047 js(k) = 2*uby - js(k) + 1
8053 cs(1) = (one-pr)*(one-qr)
8059 END SUBROUTINE getblc
8063 SUBROUTINE getbcc( GSU, I, J, PR, QR, LCMP, NS, LS, IS, JS, CS )
8068 TYPE(
t_gsu),
INTENT(IN) :: gsu
8069 INTEGER,
INTENT(IN) :: i, j
8070 REAL(8),
INTENT(IN) :: pr, qr
8071 LOGICAL,
INTENT(IN) :: lcmp
8072 INTEGER,
INTENT(OUT) :: ns
8073 LOGICAL,
POINTER,
INTENT(INOUT) :: ls(:)
8074 INTEGER,
POINTER,
INTENT(INOUT) :: is(:), js(:)
8075 REAL(8),
POINTER,
INTENT(INOUT) :: cs(:)
8078 REAL(8),
PARAMETER :: small = 1d-6
8079 LOGICAL :: ijg, llg, lclo
8080 INTEGER :: iclo, gkind
8081 INTEGER :: lbx, lby, ubx, uby, nx, ny
8082 INTEGER :: istat, p, q, ii, jj, k, l, n, m
8083 REAL(8) :: pv(0:3), qv(0:3), pw(0:3), qw(0:3)
8084 REAL(8) :: a(0:1,0:1,0:3)
8085 REAL(8) :: w(0:3,0:3) = reshape((/ 1, 0, -3, 2, &
8090 INTEGER,
PARAMETER :: nfd = 2
8091 INTEGER :: kfd(0:nfd,0:nfd) = reshape((/ 0, 1, 2, &
8095 REAL(8) :: cfd(0:nfd,0:nfd) = half* reshape((/ -3, 4, -1, &
8099 REAL(8) :: cs2d(-nfd/2:nfd,-nfd/2:nfd)
8102 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
8103 WRITE(0,
'(/2A/)')
'GETBCC ERROR -- ', &
8104 'grid search utility object not created'
8111 gkind = gsu%PTR%GKIND
8112 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
8113 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
8114 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
8117 IF (
ASSOCIATED(ls) )
THEN
8118 DEALLOCATE(ls);
NULLIFY(ls);
8120 IF (
ASSOCIATED(is) )
THEN
8121 DEALLOCATE(is);
NULLIFY(is);
8123 IF (
ASSOCIATED(js) )
THEN
8124 DEALLOCATE(js);
NULLIFY(js);
8126 IF (
ASSOCIATED(cs) )
THEN
8127 DEALLOCATE(cs);
NULLIFY(cs);
8166 a(ii,jj,0) = pw(ii) *qw(jj)
8167 a(ii,jj,1) = pw(ii+2)*qw(jj)
8168 a(ii,jj,2) = pw(ii) *qw(jj+2)
8169 a(ii,jj,3) = pw(ii+2)*qw(jj+2)
8212 IF ( mod(iclo,2).EQ.0 )
THEN
8215 IF (p-lbx.LT.nfd/2)
THEN
8217 ELSE IF (ubx-p.LT.nfd/2)
THEN
8223 IF ( mod(iclo,3).EQ.0 )
THEN
8226 IF (q-lby.LT.nfd/2)
THEN
8232 IF (q-lby.LT.nfd/2)
THEN
8234 ELSE IF (uby-q.LT.nfd/2)
THEN
8240 cs2d(ii,jj) = cs2d(ii,jj) + a(ii,jj,0)
8242 cs2d(ii+kfd(n,k),jj) = cs2d(ii+kfd(n,k),jj) &
8243 + a(ii,jj,1)*cfd(n,k)
8244 cs2d(ii,jj+kfd(n,l)) = cs2d(ii,jj+kfd(n,l)) &
8245 + a(ii,jj,2)*cfd(n,l)
8247 cs2d(ii+kfd(n,k),jj+kfd(m,l)) = &
8248 cs2d(ii+kfd(n,k),jj+kfd(m,l)) &
8249 + a(ii,jj,3)*cfd(n,k)*cfd(m,l)
8256 ns = count( abs(cs2d) .GT. small )
8257 ALLOCATE( ls(ns), is(ns), js(ns), cs(ns), stat=istat )
8258 IF ( istat .NE. 0 )
THEN
8259 WRITE(0,
'(/1A,1A/)')
'GETBCC ERROR -- ', &
8260 'array allocation failed'
8270 IF ( abs(cs2d(ii,jj)) .GT. small )
THEN
8274 cs(ns) = cs2d(ii,jj)
8275 IF ( mod(iclo,2).EQ.0 ) &
8276 is(ns) = lbx + mod(nx - 1 + mod(is(ns) - lbx + 1, nx), nx)
8277 IF ( mod(iclo,3).EQ.0 ) &
8278 js(ns) = lby + mod(ny - 1 + mod(js(ns) - lby + 1, ny), ny)
8279 IF ( iclo.EQ.
iclo_trpl .AND. js(ns).GT.uby )
THEN
8280 is(ns) = ubx + lbx - is(ns)
8281 js(ns) = 2*uby - js(ns) + 1
8287 END SUBROUTINE getbcc
8291 SUBROUTINE getgfc( GSU, I, J, PR, QR, WIDTH, LCMP, NS, LS, IS, JS, CS )
8297 TYPE(
t_gsu),
INTENT(IN) :: gsu
8298 INTEGER,
INTENT(IN) :: i, j
8299 REAL(8),
INTENT(IN) :: pr, qr
8300 REAL(8),
INTENT(IN) :: width
8301 LOGICAL,
INTENT(IN) :: lcmp
8302 INTEGER,
INTENT(OUT) :: ns
8303 LOGICAL,
POINTER,
INTENT(INOUT) :: ls(:)
8304 INTEGER,
POINTER,
INTENT(INOUT) :: is(:), js(:)
8305 REAL(8),
POINTER,
INTENT(INOUT) :: cs(:)
8310 REAL(8),
PARAMETER :: nsig = 6.0d0
8311 REAL(8),
PARAMETER :: width_min = 1.5d0
8312 LOGICAL :: ijg, llg, lclo
8313 INTEGER :: iclo, gkind
8314 INTEGER :: lbx, lby, ubx, uby, nx, ny
8316 INTEGER :: ii, jj, imin, jmin, imax, jmax
8317 REAL(8) :: wdth, sig2, rmax, r2mx, sfac, r2, gij, gsum
8320 IF ( .NOT.
ASSOCIATED(gsu%PTR) )
THEN
8321 WRITE(0,
'(/2A/)')
'GETBLC ERROR -- ', &
8322 'grid search utility object not created'
8329 gkind = gsu%PTR%GKIND
8330 lbx = gsu%PTR%LBX; lby = gsu%PTR%LBY;
8331 ubx = gsu%PTR%UBX; uby = gsu%PTR%UBY;
8332 nx = gsu%PTR%NX; ny = gsu%PTR%NY;
8333 wdth = max(width,width_min)
8334 sig2 = (wdth/nsig)**2
8338 imin = int(min(zero,pr)-rmax)
8339 jmin = int(min(zero,qr)-rmax)
8340 imax = ceiling(max(zero,pr)+rmax)
8341 jmax = ceiling(max(zero,qr)+rmax)
8344 IF (
ASSOCIATED(ls) )
THEN
8345 DEALLOCATE(ls);
NULLIFY(ls);
8347 IF (
ASSOCIATED(is) )
THEN
8348 DEALLOCATE(is);
NULLIFY(is);
8350 IF (
ASSOCIATED(js) )
THEN
8351 DEALLOCATE(js);
NULLIFY(js);
8353 IF (
ASSOCIATED(cs) )
THEN
8354 DEALLOCATE(cs);
NULLIFY(cs);
8358 ns = (imax-imin+1)*(jmax-jmin+1)
8359 ALLOCATE( ls(ns), is(ns), js(ns), cs(ns), stat=istat )
8360 IF ( istat .NE. 0 )
THEN
8361 WRITE(0,
'(/1A,1A/)')
'GETGFC ERROR -- ', &
8362 'array allocation failed'
8372 k = (imax-imin+1)*(jj-jmin) + ii - imin + 1
8377 IF ( mod(iclo,2).EQ.0 ) &
8378 is(k) = lbx + mod(nx - 1 + mod(is(k) - lbx + 1, nx), nx)
8379 IF ( mod(iclo,3).EQ.0 ) &
8380 js(k) = lby + mod(ny - 1 + mod(js(k) - lby + 1, ny), ny)
8381 IF ( iclo.EQ.
iclo_trpl .AND. js(k).GT.uby )
THEN
8382 is(k) = ubx + lbx - is(k)
8383 js(k) = 2*uby - js(k) + 1
8386 IF ( is(k).LT.lbx .OR. is(k).GT.ubx ) cycle
8387 IF ( js(k).LT.lby .OR. js(k).GT.uby ) cycle
8389 r2 = (pr - ii)**2 + (qr - jj)**2
8394 gij = exp( sfac*r2 )
8401 WHERE ( ls ) cs = cs/gsum
8404 END SUBROUTINE getgfc
8408 #define DXYDP_SINGLE_POINT_WIDE_CHANNEL_ERROR
8409 #undef DXYDP_SINGLE_POINT_WIDE_CHANNEL_WARNING
8410 SUBROUTINE dxydp( N, K, C, IJG, LLG, ICLO, &
8411 PTILED, QTILED, PRANGE, QRANGE, &
8412 LB, UB, P, Q, DXDP, DYDP, MASK, &
8413 X4, Y4, X8, Y8, RC )
8415 INTEGER,
INTENT(IN) :: n
8416 INTEGER,
INTENT(IN) :: k(0:n,0:n,1:n)
8417 REAL(8),
INTENT(IN) :: c(0:n,0:n,1:n)
8418 LOGICAL,
INTENT(IN) :: ijg
8419 LOGICAL,
INTENT(IN) :: llg
8420 INTEGER,
INTENT(IN) :: iclo
8421 LOGICAL,
INTENT(IN) :: ptiled, qtiled
8422 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
8423 INTEGER,
INTENT(IN) :: lb(2), ub(2)
8424 INTEGER,
INTENT(IN) :: p, q
8425 REAL(8),
INTENT(OUT) :: dxdp, dydp
8426 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lb(1):ub(1),lb(2):ub(2))
8427 REAL(4),
INTENT(IN),
OPTIONAL :: x4(lb(1):ub(1),lb(2):ub(2))
8428 REAL(4),
INTENT(IN),
OPTIONAL :: y4(lb(1):ub(1),lb(2):ub(2))
8429 REAL(8),
INTENT(IN),
OPTIONAL :: x8(lb(1):ub(1),lb(2):ub(2))
8430 REAL(8),
INTENT(IN),
OPTIONAL :: y8(lb(1):ub(1),lb(2):ub(2))
8431 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
8434 INTEGER,
PARAMETER :: m = 1
8435 LOGICAL,
PARAMETER :: debug = .false.
8436 CHARACTER(64) :: fstr
8437 LOGICAL :: comp_m, type_r4, type_r8
8439 INTEGER :: np, nq, lbp, lbq, ubp, ubq, p0, q0
8440 INTEGER :: istat=0, i, l, ii, ni, ii0, iin
8441 INTEGER :: kp(0:n), kq(0:n)
8447 REAL(8) :: x0, y0, lon0, lat0, c0
8448 REAL(8) :: d1dp, d2dp
8453 IF (
PRESENT(rc) ) rc = 0
8455 type_r4 =
PRESENT(x4).AND.
PRESENT(y4)
8456 type_r8 =
PRESENT(x8).AND.
PRESENT(y8)
8457 IF ( .NOT.type_r4.AND..NOT.type_r8 )
THEN
8458 WRITE(0,
'(/1A,1A/)')
'DXYDP ERROR -- ', &
8459 'no input grid coordinates specified'
8461 IF (
PRESENT(rc) )
THEN
8469 np = prange(2) - prange(1) + 1
8470 nq = qrange(2) - qrange(1) + 1
8473 lbp = lb(1); lbq = lb(2)
8474 ubp = ub(1); ubq = ub(2)
8476 lbp = lb(2); lbq = lb(1)
8477 ubp = ub(2); ubq = ub(1)
8480 IF ( p.LT.lbp .OR. p.GT.ubp .OR. q.LT.lbq .OR. q.GT.ubq )
THEN
8481 WRITE(0,
'(/1A,/1A,1L2,5I6,/1A,1L2,5I6/)')
'DXYDP ERROR -- '// &
8482 'input index coordinates outside input array bounds', &
8483 'DXYDP ERROR -- PTILED,PRANGE,P,LBP,UBP:',ptiled,prange,p,lbp,ubp, &
8484 'DXYDP ERROR -- QTILED,QRANGE,Q,LBQ,UBQ:',qtiled,qrange,q,lbq,ubq
8486 IF (
PRESENT(rc) )
THEN
8496 IF ( mod(iclo,2).EQ.0 ) &
8497 p0 = prange(1) + mod(np - 1 + mod(p0 - prange(1) + 1, np), np)
8498 IF ( mod(iclo,3).EQ.0 ) &
8499 q0 = qrange(1) + mod(nq - 1 + mod(q0 - qrange(1) + 1, nq), nq)
8500 IF ( iclo.EQ.
iclo_trpl .AND. q0.GT.qrange(2) )
THEN
8501 p0 = prange(2) + prange(1) - p0
8502 q0 = 2*qrange(2) - q0 + 1
8504 IF ( p0.LT.prange(1) .OR. p0.GT.prange(2) .OR. &
8505 q0.LT.qrange(1) .OR. q0.GT.qrange(2) )
THEN
8506 WRITE(0,
'(/1A,/1A,4I6,/1A,4I6/)')
'DXYDP ERROR -- '// &
8507 'shifted input index coordinates outside allowed range', &
8508 'DXYDP ERROR -- PRANGE,P,P0:',prange,p,p0, &
8509 'DXYDP ERROR -- QRANGE,Q,Q0:',qrange,q,q0
8511 IF (
PRESENT(rc) )
THEN
8521 comp_m =
PRESENT(mask)
8524 IF ( mask(p0,q0) )
RETURN
8526 IF ( mask(q0,p0) )
RETURN
8533 IF ( mod(iclo,2).EQ.0 )
THEN
8536 IF (p0-prange(1).LT.n/2)
THEN
8538 ELSE IF (prange(2)-p0.LT.n/2)
THEN
8539 i = n + p0 - prange(2)
8545 kp(:) = p + k(:,i,n)
8547 IF ( .NOT.ptiled )
THEN
8548 IF ( mod(iclo,2).EQ.0 )
THEN
8549 kp = prange(1) + mod(np - 1 + mod(kp - prange(1) + 1, np), np)
8553 IF ( minval(kp).LT.lbp .OR. maxval(kp).GT.ubp .OR. &
8554 minval(kq).LT.lbq .OR. maxval(kq).GT.ubq )
THEN
8555 WRITE(0,
'(/1A,/1A,1L2,8I6,/1A,1L2,8I6/)')
'DXYDP ERROR -- '// &
8556 'stencil index coordinates outside array bounds', &
8557 'DXYDP ERROR -- PTILED,PRANGE,P,P0,LBP,UBP,PMIN,PMAX:', &
8558 ptiled,prange,p,p0,lbp,ubp,minval(kp),maxval(kp), &
8559 'DXYDP ERROR -- QTILED,QRANGE,Q,Q0,LBQ,UBQ,QMIN,QMAX:', &
8560 qtiled,qrange,q,q0,lbq,ubq,minval(kq),maxval(kq)
8562 IF (
PRESENT(rc) )
THEN
8572 IF ( comp_m ) mp(l) = mask(kp(l),kq(l))
8574 xp(l) = x4(kp(l),kq(l))
8575 yp(l) = y4(kp(l),kq(l))
8577 xp(l) = x8(kp(l),kq(l))
8578 yp(l) = y8(kp(l),kq(l))
8581 IF ( comp_m ) mp(l) = mask(kq(l),kp(l))
8583 xp(l) = x4(kq(l),kp(l))
8584 yp(l) = y4(kq(l),kp(l))
8586 xp(l) = x8(kq(l),kp(l))
8587 yp(l) = y8(kq(l),kp(l))
8609 ii = count(.NOT.mp(0:i)) - 1
8610 ni = count(.NOT.mp(0:n)) - 1
8614 #ifdef DXYDP_SINGLE_POINT_WIDE_CHANNEL_ERROR
8616 WRITE(0,
'(/1A,1A,4I6/)')
'DXYDP ERROR -- ', &
8617 'single point wide channel not allowed',p,q,p0,q0
8619 IF (
PRESENT(rc) )
THEN
8630 #define DXYDP_USE_SPLX
8631 #ifdef DXYDP_USE_SPLX
8634 x0 = x4(p,q); y0 = y4(p,q);
8636 x0 = x8(p,q); y0 = y8(p,q);
8640 x0 = x4(q,p); y0 = y4(q,p);
8642 x0 = x8(q,p); y0 = y8(q,p);
8645 ihem = 1;
IF (maxval(yp(ii0:iin)).LT.zero) ihem = -1;
8646 lon0 = zero; lat0 = sign(d90,real(ihem,8));
8648 CALL w3splx(lon0,lat0,c0,xp(ii0:iin),yp(ii0:iin), &
8649 up(ii0:iin),vp(ii0:iin))
8650 d1dp = dot_product(c(0:ni,ii,ni),up(ii0:iin))
8651 d2dp = dot_product(c(0:ni,ii,ni),vp(ii0:iin))
8652 CALL spddp(lon0,c0,ihem,x0,y0,d1dp,d2dp,dxdp,dydp)
8654 dxdp = dot_product(c(0:ni,ii,ni),xp(ii0:iin))
8655 dydp = dot_product(c(0:ni,ii,ni),yp(ii0:iin))
8658 dxdp = dot_product(c(0:ni,ii,ni),xp(ii0:iin))
8659 dydp = dot_product(c(0:ni,ii,ni),yp(ii0:iin))
8662 WRITE(fstr,
'(A,I0,A,I0,A)') &
8663 '(/1A,12I8,5(/1A,2E16.8),/1A,', &
8664 ni+1,
'I16,3(/1A,',ni+1,
'E16.8))'
8665 WRITE(*,trim(fstr)) &
8666 'DXYDP -- PRANGE,QRANGE,P,Q,P0,Q0,NI,II,II0,IIN:',&
8667 prange,qrange,p,q,p0,q0,ni,ii,ii0,iin, &
8668 'DXYDP -- X0, Y0:',x0,y0, &
8669 'DXYDP -- LON0,LAT0:',lon0,lat0, &
8670 'DXYDP -- C0,IHEM:',c0,real(ihem), &
8671 'DXYDP -- D1DP,D2DP:',d1dp,d2dp, &
8672 'DXYDP -- DXDP,DYDP:',dxdp,dydp, &
8673 'DXYDP -- K:', k(0:ni,ii,ni), &
8674 'DXYDP -- C:', c(0:ni,ii,ni), &
8675 'DXYDP -- XP:',xp(ii0:iin), &
8676 'DXYDP -- YP:',yp(ii0:iin)
8679 #ifdef DXYDP_SINGLE_POINT_WIDE_CHANNEL_WARNING
8680 WRITE(0,
'(/1A,1A,4I6/)')
'DXYDP WARNING -- ', &
8681 'single point wide channel, DXDP & DYDP set to zero:',p,q,p0,q0
8687 END SUBROUTINE dxydp
8691 #define DXYDQ_SINGLE_POINT_WIDE_CHANNEL_ERROR
8692 #undef DXYDQ_SINGLE_POINT_WIDE_CHANNEL_WARNING
8693 SUBROUTINE dxydq( N, K, C, IJG, LLG, ICLO, &
8694 PTILED, QTILED, PRANGE, QRANGE, &
8695 LB, UB, P, Q, DXDQ, DYDQ, MASK, &
8696 X4, Y4, X8, Y8, RC )
8698 INTEGER,
INTENT(IN) :: n
8699 INTEGER,
INTENT(IN) :: k(0:n,0:n,1:n)
8700 REAL(8),
INTENT(IN) :: c(0:n,0:n,1:n)
8701 LOGICAL,
INTENT(IN) :: ijg
8702 LOGICAL,
INTENT(IN) :: llg
8703 INTEGER,
INTENT(IN) :: iclo
8704 LOGICAL,
INTENT(IN) :: ptiled, qtiled
8705 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
8706 INTEGER,
INTENT(IN) :: lb(2), ub(2)
8707 INTEGER,
INTENT(IN) :: p, q
8708 REAL(8),
INTENT(OUT) :: dxdq, dydq
8709 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lb(1):ub(1),lb(2):ub(2))
8710 REAL(4),
INTENT(IN),
OPTIONAL :: x4(lb(1):ub(1),lb(2):ub(2))
8711 REAL(4),
INTENT(IN),
OPTIONAL :: y4(lb(1):ub(1),lb(2):ub(2))
8712 REAL(8),
INTENT(IN),
OPTIONAL :: x8(lb(1):ub(1),lb(2):ub(2))
8713 REAL(8),
INTENT(IN),
OPTIONAL :: y8(lb(1):ub(1),lb(2):ub(2))
8714 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
8717 INTEGER,
PARAMETER :: m = 1
8718 LOGICAL,
PARAMETER :: debug = .false.
8719 CHARACTER(64) :: fstr
8720 LOGICAL :: comp_m, type_r4, type_r8
8722 INTEGER :: np, nq, lbp, lbq, ubp, ubq, p0, q0
8723 INTEGER :: istat=0, j, l, jj, nj, jj0, jjn
8724 INTEGER :: kp(0:n), kq(0:n)
8730 REAL(8) :: x0, y0, lon0, lat0, c0
8731 REAL(8) :: d1dq, d2dq
8736 IF (
PRESENT(rc) ) rc = 0
8738 type_r4 =
PRESENT(x4).AND.
PRESENT(y4)
8739 type_r8 =
PRESENT(x8).AND.
PRESENT(y8)
8740 IF ( .NOT.type_r4.AND..NOT.type_r8 )
THEN
8741 WRITE(0,
'(/1A,1A/)')
'DXYDQ ERROR -- ', &
8742 'no input grid coordinates specified'
8744 IF (
PRESENT(rc) )
THEN
8752 np = prange(2) - prange(1) + 1
8753 nq = qrange(2) - qrange(1) + 1
8756 lbp = lb(1); lbq = lb(2)
8757 ubp = ub(1); ubq = ub(2)
8759 lbp = lb(2); lbq = lb(1)
8760 ubp = ub(2); ubq = ub(1)
8763 IF ( p.LT.lbp .OR. p.GT.ubp .OR. q.LT.lbq .OR. q.GT.ubq )
THEN
8764 WRITE(0,
'(/1A,/1A,1L2,5I6,/1A,1L2,5I6/)')
'DXYDQ ERROR -- '// &
8765 'input index coordinates outside input array bounds', &
8766 'DXYDQ ERROR -- PTILED,PRANGE,P,LBP,UBP:',ptiled,prange,p,lbp,ubp, &
8767 'DXYDQ ERROR -- QTILED,QRANGE,Q,LBQ,UBQ:',qtiled,qrange,q,lbq,ubq
8769 IF (
PRESENT(rc) )
THEN
8779 IF ( mod(iclo,2).EQ.0 ) &
8780 p0 = prange(1) + mod(np - 1 + mod(p0 - prange(1) + 1, np), np)
8781 IF ( mod(iclo,3).EQ.0 ) &
8782 q0 = qrange(1) + mod(nq - 1 + mod(q0 - qrange(1) + 1, nq), nq)
8783 IF ( iclo.EQ.
iclo_trpl .AND. q0.GT.qrange(2) )
THEN
8784 p0 = prange(2) + prange(1) - p0
8785 q0 = 2*qrange(2) - q0 + 1
8787 IF ( p0.LT.prange(1) .OR. p0.GT.prange(2) .OR. &
8788 q0.LT.qrange(1) .OR. q0.GT.qrange(2) )
THEN
8789 WRITE(0,
'(/1A,/1A,4I6,/1A,4I6/)')
'DXYDQ ERROR -- '// &
8790 'shifted input index coordinates outside allowed range', &
8791 'DXYDQ ERROR -- PRANGE,P,P0:',prange,p,p0, &
8792 'DXYDQ ERROR -- QRANGE,Q,Q0:',qrange,q,q0
8794 IF (
PRESENT(rc) )
THEN
8804 comp_m =
PRESENT(mask)
8807 IF ( mask(p0,q0) )
RETURN
8809 IF ( mask(q0,p0) )
RETURN
8816 IF ( mod(iclo,3).EQ.0 )
THEN
8819 IF (q0-qrange(1).LT.n/2)
THEN
8825 IF (q0-qrange(1).LT.n/2)
THEN
8827 ELSE IF (qrange(2)-q0.LT.n/2)
THEN
8828 j = n + q0 - qrange(2)
8835 kq(:) = q + k(:,j,n)
8836 IF ( .NOT.qtiled )
THEN
8837 IF ( mod(iclo,3).EQ.0 )
THEN
8838 kq = qrange(1) + mod(nq - 1 + mod(kq - qrange(1) + 1, nq), nq)
8840 IF ( iclo.EQ.
iclo_trpl .AND. .NOT.ptiled )
THEN
8841 WHERE ( kq.GT.qrange(2) )
8842 kp = prange(2) + prange(1) - kp
8843 kq = 2*qrange(2) - kq + 1
8848 IF ( minval(kp).LT.lbp .OR. maxval(kp).GT.ubp .OR. &
8849 minval(kq).LT.lbq .OR. maxval(kq).GT.ubq )
THEN
8850 WRITE(0,
'(/1A,/1A,1L2,8I6,/1A,1L2,8I6/)')
'DXYDQ ERROR -- '// &
8851 'stencil index coordinates outside array bounds', &
8852 'DXYDQ ERROR -- PTILED,PRANGE,P,P0,LBP,UBP,PMIN,PMAX:', &
8853 ptiled,prange,p,p0,lbp,ubp,minval(kp),maxval(kp), &
8854 'DXYDQ ERROR -- QTILED,QRANGE,Q,Q0,LBQ,UBQ,QMIN,QMAX:', &
8855 qtiled,qrange,q,q0,lbq,ubq,minval(kq),maxval(kq)
8857 IF (
PRESENT(rc) )
THEN
8867 IF ( comp_m ) mq(l) = mask(kp(l),kq(l))
8869 xq(l) = x4(kp(l),kq(l))
8870 yq(l) = y4(kp(l),kq(l))
8872 xq(l) = x8(kp(l),kq(l))
8873 yq(l) = y8(kp(l),kq(l))
8876 IF ( comp_m ) mq(l) = mask(kq(l),kp(l))
8878 xq(l) = x4(kq(l),kp(l))
8879 yq(l) = y4(kq(l),kp(l))
8881 xq(l) = x8(kq(l),kp(l))
8882 yq(l) = y8(kq(l),kp(l))
8904 jj = count(.NOT.mq(0:j)) - 1
8905 nj = count(.NOT.mq(0:n)) - 1
8909 #ifdef DXYDQ_SINGLE_POINT_WIDE_CHANNEL_ERROR
8911 WRITE(0,
'(/1A,1A,4I6/)')
'DXYDQ ERROR -- ', &
8912 'single point wide channel not allowed',p,q,p0,q0
8914 IF (
PRESENT(rc) )
THEN
8925 #define DXYDQ_USE_SPLX
8926 #ifdef DXYDQ_USE_SPLX
8929 x0 = x4(p,q); y0 = y4(p,q);
8931 x0 = x8(p,q); y0 = y8(p,q);
8935 x0 = x4(q,p); y0 = y4(q,p);
8937 x0 = x8(q,p); y0 = y8(q,p);
8940 ihem = 1;
IF (maxval(yq(jj0:jjn)).LT.zero) ihem = -1;
8941 lon0 = zero; lat0 = sign(d90,real(ihem,8));
8943 CALL w3splx(lon0,lat0,c0,xq(jj0:jjn),yq(jj0:jjn), &
8944 uq(jj0:jjn),vq(jj0:jjn))
8945 d1dq = dot_product(c(0:nj,jj,nj),uq(jj0:jjn))
8946 d2dq = dot_product(c(0:nj,jj,nj),vq(jj0:jjn))
8947 CALL spddq(lon0,c0,ihem,x0,y0,d1dq,d2dq,dxdq,dydq)
8949 dxdq = dot_product(c(0:nj,jj,nj),xq(jj0:jjn))
8950 dydq = dot_product(c(0:nj,jj,nj),yq(jj0:jjn))
8953 dxdq = dot_product(c(0:nj,jj,nj),xq(jj0:jjn))
8954 dydq = dot_product(c(0:nj,jj,nj),yq(jj0:jjn))
8957 WRITE(fstr,
'(A,I0,A,I0,A)') &
8958 '(/1A,12I8,5(/1A,2E16.8),/1A,', &
8959 nj+1,
'I16,3(/1A,',nj+1,
'E16.8))'
8960 WRITE(*,trim(fstr)) &
8961 'DXYDQ -- PRANGE,QRANGE,P,Q,P0,Q0,NJ,JJ,JJ0,JJN:',&
8962 prange,qrange,p,q,p0,q0,nj,jj,jj0,jjn, &
8963 'DXYDQ -- X0, Y0:',x0,y0, &
8964 'DXYDQ -- LON0,LAT0:',lon0,lat0, &
8965 'DXYDQ -- C0,IHEM:',c0,real(ihem), &
8966 'DXYDQ -- D1DQ,D1DQ:',d1dq,d1dq, &
8967 'DXYDQ -- DXDQ,DYDQ:',dxdq,dydq, &
8968 'DXYDQ -- K:', k(0:nj,jj,nj), &
8969 'DXYDQ -- C:', c(0:nj,jj,nj), &
8970 'DXYDQ -- XQ:',xq(jj0:jjn), &
8971 'DXYDQ -- YQ:',yq(jj0:jjn)
8974 #ifdef DXYDQ_SINGLE_POINT_WIDE_CHANNEL_WARNING
8975 WRITE(0,
'(/1A,1A,4I6/)')
'DXYDQ WARNING -- ', &
8976 'single point wide channel, DXDQ & DYDQ set to zero:',p,q,p0,q0
8982 END SUBROUTINE dxydq
8986 SUBROUTINE spddp( LAM0, C0, IHEM, LAM, PHI, DXDP, DYDP, &
9006 REAL(8),
INTENT(IN) :: lam0, c0
9007 INTEGER,
INTENT(IN) :: ihem
9008 REAL(8),
INTENT(IN) :: lam, phi
9009 REAL(8),
INTENT(IN) :: dxdp, dydp
9010 REAL(8),
INTENT(OUT):: dlamdp, dphidp
9013 REAL(8),
PARAMETER :: small = 1d-6
9014 REAL(8) :: k0, a, mu, nu, fac
9015 REAL(8) :: cosmu, sinmu, cosnu2, cotnu
9016 REAL(8) :: dlamdx, dlamdy, dphidx, dphidy
9018 k0 = cos(half*c0*d2r)**2
9020 a = sign(one,real(ihem,8))
9021 nu = pio4 - a*half*phi*d2r
9022 nu = sign(max(small,abs(nu)),nu)
9023 fac = r2d*half/rearth/k0
9030 dlamdx = fac*cotnu*cosmu
9031 dlamdy = a*fac*cotnu*sinmu
9032 dphidx = -a*two*fac*cosnu2*sinmu
9033 dphidy = two*fac*cosnu2*cosmu
9035 dlamdp = dxdp*dlamdx + dydp*dlamdy
9036 dphidp = dxdp*dphidx + dydp*dphidy
9038 END SUBROUTINE spddp
9042 SUBROUTINE spddq( LAM0, C0, IHEM, LAM, PHI, DXDQ, DYDQ, &
9062 REAL(8),
INTENT(IN) :: lam0, c0
9063 INTEGER,
INTENT(IN) :: ihem
9064 REAL(8),
INTENT(IN) :: lam, phi
9065 REAL(8),
INTENT(IN) :: dxdq, dydq
9066 REAL(8),
INTENT(OUT):: dlamdq, dphidq
9069 REAL(8),
PARAMETER :: small = 1d-6
9070 REAL(8) :: k0, a, mu, nu, fac
9071 REAL(8) :: cosmu, sinmu, cosnu2, cotnu
9072 REAL(8) :: dlamdx, dlamdy, dphidx, dphidy
9074 k0 = cos(half*c0*d2r)**2
9076 a = sign(one,real(ihem,8))
9077 nu = pio4 - a*half*phi*d2r
9078 nu = sign(max(small,abs(nu)),nu)
9079 fac = r2d*half/rearth/k0
9086 dlamdx = fac*cotnu*cosmu
9087 dlamdy = a*fac*cotnu*sinmu
9088 dphidx = -a*two*fac*cosnu2*sinmu
9089 dphidy = two*fac*cosnu2*cosmu
9091 dlamdq = dxdq*dlamdx + dydq*dlamdy
9092 dphidq = dxdq*dphidx + dydq*dphidy
9094 END SUBROUTINE spddq
9098 #undef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_ERROR
9099 #undef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_WARNING
9100 SUBROUTINE dfdpq( N, K, C, IJG, ICLO, &
9104 F4, F8, DFDP, DFDQ, &
9105 G4, G8, DGDP, DGDQ, &
9106 H4, H8, DHDP, DHDQ, &
9107 NSDP, ISDP, JSDP, CSDP, &
9108 NSDQ, ISDQ, JSDQ, CSDQ, &
9111 INTEGER,
INTENT(IN) :: n
9112 INTEGER,
INTENT(IN) :: k(0:n,0:n,1:n)
9113 REAL(8),
INTENT(IN) :: c(0:n,0:n,1:n)
9114 LOGICAL,
INTENT(IN) :: ijg
9115 INTEGER,
INTENT(IN) :: iclo
9116 LOGICAL,
INTENT(IN) :: ptiled, qtiled
9117 INTEGER,
INTENT(IN) :: prange(2), qrange(2)
9118 INTEGER,
INTENT(IN) :: lb(2), ub(2)
9119 INTEGER,
INTENT(IN) :: p, q
9120 REAL(4),
INTENT(IN),
OPTIONAL :: f4(lb(1):ub(1),lb(2):ub(2))
9121 REAL(8),
INTENT(IN),
OPTIONAL :: f8(lb(1):ub(1),lb(2):ub(2))
9122 REAL(8),
INTENT(OUT),
OPTIONAL :: dfdp, dfdq
9123 REAL(4),
INTENT(IN),
OPTIONAL :: g4(lb(1):ub(1),lb(2):ub(2))
9124 REAL(8),
INTENT(IN),
OPTIONAL :: g8(lb(1):ub(1),lb(2):ub(2))
9125 REAL(8),
INTENT(OUT),
OPTIONAL :: dgdp, dgdq
9126 REAL(4),
INTENT(IN),
OPTIONAL :: h4(lb(1):ub(1),lb(2):ub(2))
9127 REAL(8),
INTENT(IN),
OPTIONAL :: h8(lb(1):ub(1),lb(2):ub(2))
9128 REAL(8),
INTENT(OUT),
OPTIONAL :: dhdp, dhdq
9129 INTEGER,
INTENT(OUT),
OPTIONAL :: nsdp
9130 INTEGER,
POINTER,
OPTIONAL :: isdp(:), jsdp(:)
9131 REAL(8),
POINTER,
OPTIONAL :: csdp(:)
9132 INTEGER,
INTENT(OUT),
OPTIONAL :: nsdq
9133 INTEGER,
POINTER,
OPTIONAL :: isdq(:), jsdq(:)
9134 REAL(8),
POINTER,
OPTIONAL :: csdq(:)
9135 LOGICAL,
INTENT(IN),
OPTIONAL :: mask(lb(1):ub(1),lb(2):ub(2))
9136 INTEGER,
INTENT(OUT),
OPTIONAL :: rc
9139 INTEGER,
PARAMETER :: m = 1
9140 LOGICAL,
PARAMETER :: debug = .false.
9141 CHARACTER(64) :: fstr
9142 LOGICAL :: comp_m, comp_f, comp_g, comp_h, type_r4
9143 LOGICAL :: comp_cp, comp_cq
9144 INTEGER :: np, nq, lbp, lbq, ubp, ubq, p0, q0
9145 INTEGER :: istat=0, i, j, l, ii, jj, ni, nj, ii0, iin, jj0, jjn
9146 INTEGER :: kp(0:n), kq(0:n)
9147 LOGICAL :: mp(0:n), mq(0:n)
9148 REAL(8) :: fp(0:n), fq(0:n)
9149 REAL(8) :: gp(0:n), gq(0:n)
9150 REAL(8) :: hp(0:n), hq(0:n)
9151 INTEGER :: ip(0:n), iq(0:n)
9152 INTEGER :: jp(0:n), jq(0:n)
9157 IF (
PRESENT(rc) ) rc = 0
9159 comp_f = (
PRESENT(f4) .OR.
PRESENT(f8) ) .AND. &
9160 PRESENT(dfdp) .AND.
PRESENT(dfdq)
9161 comp_g = (
PRESENT(g4) .OR.
PRESENT(g8) ) .AND. &
9162 PRESENT(dgdp) .AND.
PRESENT(dgdq)
9163 comp_h = (
PRESENT(h4) .OR.
PRESENT(h8) ) .AND. &
9164 PRESENT(dhdp) .AND.
PRESENT(dhdq)
9165 comp_cp =
PRESENT(nsdp) .AND.
PRESENT(isdp) .AND. &
9166 PRESENT(jsdp) .AND.
PRESENT(csdp)
9167 comp_cq =
PRESENT(nsdq) .AND.
PRESENT(isdq) .AND. &
9168 PRESENT(jsdq) .AND.
PRESENT(csdq)
9169 IF ( .NOT.comp_f.AND..NOT.comp_g.AND..NOT.comp_h.AND. &
9170 .NOT.comp_cp.AND..NOT.comp_cq )
RETURN
9173 type_r4 =
PRESENT(f4)
9174 ELSE IF ( comp_g )
THEN
9175 type_r4 =
PRESENT(g4)
9176 ELSE IF ( comp_h )
THEN
9177 type_r4 =
PRESENT(h4)
9180 np = prange(2) - prange(1) + 1
9181 nq = qrange(2) - qrange(1) + 1
9184 lbp = lb(1); lbq = lb(2)
9185 ubp = ub(1); ubq = ub(2)
9187 lbp = lb(2); lbq = lb(1)
9188 ubp = ub(2); ubq = ub(1)
9191 IF ( p.LT.lbp .OR. p.GT.ubp .OR. q.LT.lbq .OR. q.GT.ubq )
THEN
9192 WRITE(0,
'(/1A,/1A,1L2,5I6,/1A,1L2,5I6/)')
'DFDPQ ERROR -- '// &
9193 'input index coordinates outside input array bounds', &
9194 'DFDPQ ERROR -- PTILED,PRANGE,P,LBP,UBP:',ptiled,prange,p,lbp,ubp, &
9195 'DFDPQ ERROR -- QTILED,QRANGE,Q,LBQ,UBQ:',qtiled,qrange,q,lbq,ubq
9197 IF (
PRESENT(rc) )
THEN
9207 IF ( mod(iclo,2).EQ.0 ) &
9208 p0 = prange(1) + mod(np - 1 + mod(p0 - prange(1) + 1, np), np)
9209 IF ( mod(iclo,3).EQ.0 ) &
9210 q0 = qrange(1) + mod(nq - 1 + mod(q0 - qrange(1) + 1, nq), nq)
9211 IF ( iclo.EQ.
iclo_trpl .AND. q0.GT.qrange(2) )
THEN
9212 p0 = prange(2) + prange(1) - p0
9213 q0 = 2*qrange(2) - q0 + 1
9215 IF ( p0.LT.prange(1) .OR. p0.GT.prange(2) .OR. &
9216 q0.LT.qrange(1) .OR. q0.GT.qrange(2) )
THEN
9217 WRITE(0,
'(/1A,/1A,4I6,/1A,4I6/)')
'DFDPQ ERROR -- '// &
9218 'shifted input index coordinates outside allowed range', &
9219 'DFDPQ ERROR -- PRANGE,P,P0:',prange,p,p0, &
9220 'DFDPQ ERROR -- QRANGE,Q,Q0:',qrange,q,q0
9222 IF (
PRESENT(rc) )
THEN
9230 comp_m =
PRESENT(mask)
9233 IF ( mask(p0,q0) )
RETURN
9235 IF ( mask(q0,p0) )
RETURN
9242 IF ( comp_f.OR.comp_g.OR.comp_h.OR.comp_cp )
THEN
9244 IF ( mod(iclo,2).EQ.0 )
THEN
9247 IF (p0-prange(1).LT.n/2)
THEN
9249 ELSE IF (prange(2)-p0.LT.n/2)
THEN
9250 i = n + p0 - prange(2)
9256 kp(:) = p + k(:,i,n)
9258 IF ( .NOT.ptiled )
THEN
9259 IF ( mod(iclo,2).EQ.0 )
THEN
9260 kp = prange(1) + mod(np - 1 + mod(kp - prange(1) + 1, np), np)
9264 IF ( minval(kp).LT.lbp .OR. maxval(kp).GT.ubp .OR. &
9265 minval(kq).LT.lbq .OR. maxval(kq).GT.ubq )
THEN
9266 WRITE(0,
'(/1A,/1A,1L2,8I6,/1A,1L2,8I6/)')
'DFDPQ ERROR -- '// &
9267 'stencil index coordinates outside array bounds', &
9268 'DFDPQ ERROR -- PTILED,PRANGE,P,P0,LBP,UBP,PMIN,PMAX:', &
9269 ptiled,prange,p,p0,lbp,ubp,minval(kp),maxval(kp), &
9270 'DFDPQ ERROR -- QTILED,QRANGE,Q,Q0,LBQ,UBQ,QMIN,QMAX:', &
9271 qtiled,qrange,q,q0,lbq,ubq,minval(kq),maxval(kq)
9273 IF (
PRESENT(rc) )
THEN
9282 ip(:) = p0 + k(:,i,n)
9284 IF ( mod(iclo,2).EQ.0 )
THEN
9285 ip = prange(1) + mod(np - 1 + mod(ip - prange(1) + 1, np), np)
9291 IF ( comp_m ) mp(l) = mask(kp(l),kq(l))
9293 IF ( comp_f ) fp(l) = f4(kp(l),kq(l))
9294 IF ( comp_g ) gp(l) = g4(kp(l),kq(l))
9295 IF ( comp_h ) hp(l) = h4(kp(l),kq(l))
9297 IF ( comp_f ) fp(l) = f8(kp(l),kq(l))
9298 IF ( comp_g ) gp(l) = g8(kp(l),kq(l))
9299 IF ( comp_h ) hp(l) = h8(kp(l),kq(l))
9302 IF ( comp_m ) mp(l) = mask(kq(l),kp(l))
9304 IF ( comp_f ) fp(l) = f4(kq(l),kp(l))
9305 IF ( comp_g ) gp(l) = g4(kq(l),kp(l))
9306 IF ( comp_h ) hp(l) = h4(kq(l),kp(l))
9308 IF ( comp_f ) fp(l) = f8(kq(l),kp(l))
9309 IF ( comp_g ) gp(l) = g8(kq(l),kp(l))
9310 IF ( comp_h ) hp(l) = h8(kq(l),kp(l))
9332 ii = count(.NOT.mp(0:i)) - 1
9333 ni = count(.NOT.mp(0:n)) - 1
9337 #ifdef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_ERROR
9339 WRITE(0,
'(/1A,1A,4I6/)')
'DFDPQ ERROR -- ', &
9340 'DFDP -- single point wide channel not allowed',p,q,p0,q0
9342 IF (
PRESENT(rc) )
THEN
9352 IF ( comp_f ) dfdp = dot_product(c(0:ni,ii,ni),fp(ii0:iin))
9353 IF ( comp_g ) dgdp = dot_product(c(0:ni,ii,ni),gp(ii0:iin))
9354 IF ( comp_h ) dhdp = dot_product(c(0:ni,ii,ni),hp(ii0:iin))
9356 IF (
ASSOCIATED(isdp) )
DEALLOCATE(isdp)
9357 IF (
ASSOCIATED(jsdp) )
DEALLOCATE(jsdp)
9358 IF (
ASSOCIATED(csdp) )
DEALLOCATE(csdp)
9360 ALLOCATE(isdp(nsdp),jsdp(nsdp),csdp(nsdp))
9361 isdp(1:nsdp) = ip(ii0:iin)
9362 jsdp(1:nsdp) = jp(ii0:iin)
9363 csdp(1:nsdp) = c(0:ni,ii,ni)
9365 IF ( debug .AND. comp_f )
THEN
9366 WRITE(fstr,
'(A,I0,A,I0,A,I0,A)')
'(/1A,8I6,E16.8,/1A,',&
9367 ni+1,
'I16,/1A,',ni+1,
'E16.8,/1A,',ni+1,
'E16.8)'
9368 WRITE(*,trim(fstr)) &
9369 'DFDPQ -- DFDP -- P,Q,P0,Q0,NI,II,II0,IIN,DFDP:',&
9370 p,q,p0,q0,ni,ii,ii0,iin,dfdp, &
9371 'DFDPQ -- DFDP -- K:', k(0:ni,ii,ni), &
9372 'DFDPQ -- DFDP -- C:', c(0:ni,ii,ni), &
9373 'DFDPQ -- DFDP -- FP:', fp(ii0:iin)
9376 #ifdef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_WARNING
9377 WRITE(0,
'(/1A,1A,4I6/)')
'DFDPQ WARNING -- ', &
9378 'single point wide channel, DFDP set to zero:',p,q,p0,q0
9380 IF ( comp_f ) dfdp = zero
9381 IF ( comp_g ) dgdp = zero
9382 IF ( comp_h ) dhdp = zero
9383 IF ( comp_cp ) nsdp = 0
9391 IF ( comp_f.OR.comp_g.OR.comp_h.OR.comp_cq )
THEN
9393 IF ( mod(iclo,3).EQ.0 )
THEN
9396 IF (q0-qrange(1).LT.n/2)
THEN
9402 IF (q0-qrange(1).LT.n/2)
THEN
9404 ELSE IF (qrange(2)-q0.LT.n/2)
THEN
9405 j = n + q0 - qrange(2)
9412 kq(:) = q + k(:,j,n)
9413 IF ( .NOT.qtiled )
THEN
9414 IF ( mod(iclo,3).EQ.0 )
THEN
9415 kq = qrange(1) + mod(nq - 1 + mod(kq - qrange(1) + 1, nq), nq)
9417 IF ( iclo.EQ.
iclo_trpl .AND. .NOT.ptiled )
THEN
9418 WHERE ( kq.GT.qrange(2) )
9419 kp = prange(2) + prange(1) - kp
9420 kq = 2*qrange(2) - kq + 1
9425 IF ( minval(kp).LT.lbp .OR. maxval(kp).GT.ubp .OR. &
9426 minval(kq).LT.lbq .OR. maxval(kq).GT.ubq )
THEN
9427 WRITE(0,
'(/1A,/1A,1L2,8I6,/1A,1L2,8I6/)')
'DFDPQ ERROR -- '// &
9428 'stencil index coordinates outside array bounds', &
9429 'DFDPQ ERROR -- PTILED,PRANGE,P,P0,LBP,UBP,PMIN,PMAX:', &
9430 ptiled,prange,p,p0,lbp,ubp,minval(kp),maxval(kp), &
9431 'DFDPQ ERROR -- QTILED,QRANGE,Q,Q0,LBQ,UBQ,QMIN,QMAX:', &
9432 qtiled,qrange,q,q0,lbq,ubq,minval(kq),maxval(kq)
9434 IF (
PRESENT(rc) )
THEN
9444 jq(:) = q0 + k(:,j,n)
9445 IF ( mod(iclo,3).EQ.0 )
THEN
9446 jq = qrange(1) + mod(nq - 1 + mod(jq - qrange(1) + 1, nq), nq)
9449 WHERE ( jq.GT.qrange(2) )
9450 iq = prange(2) + prange(1) - iq
9451 jq = 2*qrange(2) - jq + 1
9458 IF ( comp_m ) mq(l) = mask(kp(l),kq(l))
9460 IF ( comp_f ) fq(l) = f4(kp(l),kq(l))
9461 IF ( comp_g ) gq(l) = g4(kp(l),kq(l))
9462 IF ( comp_h ) hq(l) = h4(kp(l),kq(l))
9464 IF ( comp_f ) fq(l) = f8(kp(l),kq(l))
9465 IF ( comp_g ) gq(l) = g8(kp(l),kq(l))
9466 IF ( comp_h ) hq(l) = h8(kp(l),kq(l))
9469 IF ( comp_m ) mq(l) = mask(kq(l),kp(l))
9471 IF ( comp_f ) fq(l) = f4(kq(l),kp(l))
9472 IF ( comp_g ) gq(l) = g4(kq(l),kp(l))
9473 IF ( comp_h ) hq(l) = h4(kq(l),kp(l))
9475 IF ( comp_f ) fq(l) = f8(kq(l),kp(l))
9476 IF ( comp_g ) gq(l) = g8(kq(l),kp(l))
9477 IF ( comp_h ) hq(l) = h8(kq(l),kp(l))
9499 jj = count(.NOT.mq(0:j)) - 1
9500 nj = count(.NOT.mq(0:n)) - 1
9504 #ifdef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_ERROR
9506 WRITE(0,
'(/1A,1A,4I6/)')
'DFDPQ ERROR -- ', &
9507 'DFDQ -- single point wide channel not allowed',p,q,p0,q0
9509 IF (
PRESENT(rc) )
THEN
9519 IF ( comp_f ) dfdq = dot_product(c(0:nj,jj,nj),fq(jj0:jjn))
9520 IF ( comp_g ) dgdq = dot_product(c(0:nj,jj,nj),gq(jj0:jjn))
9521 IF ( comp_h ) dhdq = dot_product(c(0:nj,jj,nj),hq(jj0:jjn))
9523 IF (
ASSOCIATED(isdq) )
DEALLOCATE(isdq)
9524 IF (
ASSOCIATED(jsdq) )
DEALLOCATE(jsdq)
9525 IF (
ASSOCIATED(csdq) )
DEALLOCATE(csdq)
9527 ALLOCATE(isdq(nsdq),jsdq(nsdq),csdq(nsdq))
9528 isdq(1:nsdq) = iq(jj0:jjn)
9529 jsdq(1:nsdq) = jq(jj0:jjn)
9530 csdq(1:nsdq) = c(0:nj,jj,nj)
9532 IF ( debug .AND. comp_f )
THEN
9533 WRITE(fstr,
'(A,I0,A,I0,A,I0,A)')
'(/1A,8I6,E16.8,/1A,',&
9534 nj+1,
'I16,/1A,',nj+1,
'E16.8,/1A,',nj+1,
'E16.8)'
9535 WRITE(*,trim(fstr)) &
9536 'DFDPQ -- DFDQ -- P,Q,P0,Q0,NJ,JJ,JJ0,JJN,DFDQ:',&
9537 p,q,p0,q0,nj,jj,jj0,jjn,dfdq, &
9538 'DFDPQ -- DFDQ -- K:', k(0:nj,jj,nj), &
9539 'DFDPQ -- DFDQ -- C:', c(0:nj,jj,nj), &
9540 'DFDPQ -- DFDQ -- FQ:', fq(jj0:jjn)
9543 #ifdef DFDPQ_SINGLE_POINT_WIDE_CHANNEL_WARNING
9544 WRITE(0,
'(/1A,1A,4I6/)')
'DFDPQ WARNING -- ', &
9545 'single point wide channel, DFDQ set to zero:',p,q,p0,q0
9547 IF ( comp_f ) dfdq = zero
9548 IF ( comp_g ) dgdq = zero
9549 IF ( comp_h ) dhdq = zero
9550 IF ( comp_cq ) nsdq = 0
9555 END SUBROUTINE dfdpq
9559 SUBROUTINE get_fdw2( N, M, K, C )
9561 INTEGER,
INTENT(IN) :: n, m
9562 INTEGER,
INTENT(OUT):: k(0:n,0:n)
9563 REAL(8),
INTENT(OUT):: c(0:n,0:n)
9565 REAL(8) :: a(0:n), b(0:n,0:m)
9572 CALL w3fdwt( n, n, m, zero, a, b )
9578 END SUBROUTINE get_fdw2
9582 SUBROUTINE get_fdw3( N, M, K, C )
9584 INTEGER,
INTENT(IN) :: n, m
9585 INTEGER,
INTENT(OUT):: k(0:n,0:n,1:n)
9586 REAL(8),
INTENT(OUT):: c(0:n,0:n,1:n)
9588 REAL(8) :: a(0:n), b(0:n,0:m)
9597 CALL w3fdwt( l, n, m, zero, a, b )
9598 c(0:l,i,l) = b(0:l,m)
9604 END SUBROUTINE get_fdw3
9624 SUBROUTINE extcde(IEXIT)
9628 INTEGER,
INTENT(IN) :: iexit
9632 CALL mpi_initialized ( run, ierr_mpi )
9634 CALL mpi_abort ( mpi_comm_world, iexit, ierr_mpi )
9638 END SUBROUTINE extcde