122 real*8,
ALLOCATABLE ::
nm(:,:,:),
dtsi(:)
123 INTEGER,
ALLOCATABLE ::
iter(:)
221 INTEGER,
SAVE :: IENT = 0
228 INTEGER :: I, J, IBND_MAP, ISEA, IP, IX, JSEA, nb
230 INTEGER :: myrank, ierr, iproc
231 INTEGER,
ALLOCATABLE :: NSEAL_arr(:)
234 INTEGER,
INTENT(in) :: IMOD
236 INTEGER IK0, ISP0, ITH
238 REAL,
PARAMETER :: COEF4 = 5.0e-7
240 CALL strace (ient,
'PDLIB_INIT')
242 #ifdef W3_DEBUGSOLVER
243 WRITE(740+
iaproc,*)
'PDLIB_INIT, IMOD (no print)'
255 #ifdef W3_DEBUGSOLVER
258 WRITE(740+
iaproc,*)
'PDLIB_INIT, myrank=', myrank
268 #ifdef W3_DEBUGSOLVER
269 WRITE(740+
iaproc,*)
'After initFromGridDim'
280 IF (isea .gt. 0) pdlib_nseal = pdlib_nseal + 1
282 #ifdef W3_DEBUGSOLVER
283 WRITE(740+
iaproc,*)
'npa is augmented domain over NX'
284 WRITE(740+
iaproc,*)
'PDLIB_NSEAL is basicall npa but only over the wet points'
285 WRITE(740+
iaproc,*)
'NSEAL is set to PDLIB_NSEAL'
286 WRITE(740+
iaproc,*)
'PDLIB_NSEAL=', pdlib_nseal
287 WRITE(740+
iaproc,*)
'NSEAL =', nseal,
'NP =',
np,
'NPA =',
npa
291 #ifdef W3_DEBUGSOLVER
292 WRITE(740+
iaproc,*)
'ISEA_TO_JSEA ALLOCATEd'
302 IF (isea .gt. 0)
THEN
309 #ifdef W3_DEBUGSOLVER
310 WRITE(740+
iaproc,*)
'After JX_TO_JSEA, ISEA_TO_JSEA and friend computation'
318 IF (mapfs(1,ix) .gt. 0) nb = nb + 1
321 IF (nb .ne. nsea)
THEN
322 WRITE(*,*)
'Logical error in computation of NSEA / nb'
323 WRITE(*,*)
'nb=', nb,
' NSEA=', nsea
326 #ifdef W3_DEBUGSOLVER
327 WRITE(740+
iaproc,*)
'nb / NSEA consistency check'
332 IF ((flcx .eqv. .true.).and.(flcy .eqv. .true.))
THEN
340 ALLOCATE(nseal_arr(
naproc))
341 nseal_arr(1)=pdlib_nseal
344 nseal_arr(iproc)=iscal(1)
346 pdlib_nsealm=maxval(nseal_arr)
347 DEALLOCATE(nseal_arr)
350 CALL mpi_send(iscal,1,mpi_int, 0, 23,
mpi_comm_wave, ierr_mpi)
355 iscal(1)=pdlib_nsealm
357 CALL mpi_send(iscal,1,mpi_int, iproc-1, 24,
mpi_comm_wave, ierr_mpi)
361 pdlib_nsealm=iscal(1)
366 WRITE(740+
iaproc,*)
'PDLIB_NSEALM=', pdlib_nsealm
385 DO jsea=1, pdlib_nseal
388 isea = mapfs(1,ip_glob)
389 IF (isea .ne. ip_glob)
THEN
390 WRITE(*,*) jsea, pdlib_nseal, ip, ip_glob, isea
391 WRITE(*,*) .ne.
'ISEA IP_glob'
466 INTEGER,
SAVE :: IENT = 0
471 INTEGER :: IBND_MAP, ISEA, JSEA, IX, IP, IP_glob
472 INTEGER,
INTENT(in) :: IMOD
473 INTEGER :: Status(NX), istat
476 CALL strace (ient,
'PDLIB_MAPSTA_INIT')
479 WRITE(*,*)
'Passing by PDLIB_MAPSTA_INIT IAPROC=',
iaproc
481 IF (
iaproc .gt. naproc)
THEN
485 ALLOCATE(grids(imod)%MAPSTA_LOC(
npa), stat=istat)
497 IF ((mapsta(1,ix) .lt. 1).and.(status(ix).gt.0))
THEN
502 ALLOCATE(grids(imod)%INDEX_MAP(
nbnd_map), stat=istat)
507 IF ((mapsta(1,ix) .lt. 1).and.(status(ix).gt.0))
THEN
508 ibnd_map = ibnd_map + 1
582 INTEGER,
SAVE :: IENT = 0
587 INTEGER :: IBND_MAP, ISEA, JSEA, IX, IP, IP_glob
588 INTEGER,
INTENT(in) :: IMOD
589 INTEGER :: Status(NX), istat
592 CALL strace (ient,
'PDLIB_MAPSTA_INIT')
595 WRITE(*,*)
'Passing by PDLIB_MAPSTA_INIT IAPROC=',
iaproc
597 IF (
iaproc .gt. naproc)
THEN
601 ALLOCATE(grids(imod)%IOBP_LOC(
npa), stat=istat)
603 ALLOCATE(grids(imod)%IOBPD_LOC(nth,
npa), stat=istat)
605 ALLOCATE(grids(imod)%IOBDP_LOC(
npa), stat=istat)
607 ALLOCATE(grids(imod)%IOBPA_LOC(
npa), stat=istat)
624 DEALLOCATE(grids(imod)%IOBP,grids(imod)%IOBPD)
631 SUBROUTINE pdlib_w3xypug ( ISP, FACX, FACY, DTG, VGX, VGY, LCALC )
713 INTEGER,
INTENT(IN) :: ISP
714 REAL,
INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
715 LOGICAL,
INTENT(IN) :: LCALC
721 INTEGER :: ITH, IK, ISEA
722 INTEGER :: I, J, IE, IBND_MAP
724 REAL :: CCOS, CSIN, CCURX, CCURY, WN1, CG1
730 REAL :: VLCFLX(npa), VLCFLY(npa)
732 REAL :: AC_MAP(NBND_MAP)
742 #ifdef W3_DEBUGSOLVER
743 WRITE(740+
iaproc,*)
'Begin of PDLIB_W3XYPUG'
746 ith = 1 + mod(isp-1,nth)
748 ccos = facx * ecos(ith)
749 csin = facy * esin(ith)
764 isea = mapfs(1,ip_glob)
767 ac(ip) =
va(isp,jsea) / cg1 * clats(isea)
768 vlcflx(ip) = ccos * cg1 / clats(isea)
769 vlcfly(ip) = csin * cg1
771 ac(ip) =
va(isp,jsea) /
cg(ik,isea) * clats(isea)
772 vlcflx(ip) = ccos *
cg(ik,isea) / clats(isea)
773 vlcfly(ip) = csin *
cg(ik,isea)
776 vlcflx(ip) = vlcflx(ip) - ccurx*vgx/clats(isea)
777 vlcfly(ip) = vlcfly(ip) - ccury*vgy
781 #ifdef W3_DEBUGSOLVER
782 WRITE(740+
iaproc,*)
'ISP=', isp,
' ITH=', ith,
' IK=', ik
783 WRITE(740+
iaproc,*)
'1: maxval(VLCFLX)=', maxval(vlcflx)
784 WRITE(740+
iaproc,*)
'1: maxval(VLCFLY)=', maxval(vlcfly)
792 isea = mapfs(1,ip_glob)
796 IF (iobp_loc(ip) .GT. 0)
THEN
797 vlcflx(ip) = vlcflx(ip) + ccurx*
cx(isea)/clats(isea)
798 vlcfly(ip) = vlcfly(ip) + ccury*
cy(isea)
803 c(:,1) = vlcflx(:) * iobdp_loc
804 c(:,2) = vlcfly(:) * iobdp_loc
810 rd2 =
dsec21( tbpi0, tbpin )
819 DO ibnd_map=1,nbnd_map
820 ip=index_map(ibnd_map)
821 ac_map(ibnd_map) = ac(ip)
827 #ifdef W3_DEBUGSOLVER
828 WRITE(740+
iaproc,*)
'maxval(C)=', maxval(c)
837 ELSE IF (fsnimp)
THEN
838 stop
'For PDLIB and FSNIMP, no function has been programmed yet'
842 DO ibnd_map=1,nbnd_map
843 ip=index_map(ibnd_map)
844 ac(ip) = ac_map(ibnd_map)
847 #ifdef W3_DEBUGSOLVER
848 WRITE(740+
iaproc,*)
'After solutioning'
857 isea=mapfs(1,ip_glob)
858 va(isp,jsea) = max( 0. ,
cg(ik,isea)/clats(isea)*ac(ip) )
860 #ifdef W3_DEBUGSOLVER
861 WRITE(740+
iaproc,*)
'Leaving PDLIB_W3XYPUG'
934 USE mpi,
only : mpi_min
939 INTEGER,
INTENT(IN) :: ISP
941 REAL,
INTENT(IN) :: DT
944 REAL,
INTENT(IN) :: C(npa,2)
946 REAL,
INTENT(INOUT) :: AC(npa)
948 REAL,
INTENT(IN) :: RD10, RD20
951 LOGICAL,
INTENT(IN) :: LCALC
954 INTEGER,
SAVE :: IENT = 0
957 INTEGER(KIND=1) :: IOBPDR(NX)
959 INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
960 INTEGER :: IBI, NI(3)
972 REAL :: FL11, FL12, FL21, FL22, FL31, FL32
973 REAL :: FL111, FL112, FL211, FL212, FL311, FL312
974 REAL :: DTSI(npa), U(npa)
975 REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST
976 REAL :: LAMBDA(2), KTMP(3)
977 REAL :: KELEM(3,NE), FLALL(3,NE)
978 REAL :: KKSUM(npa), ST(npa)
980 INTEGER :: ISPROC, JSEA, IP_glob, ierr, IX
981 REAL :: eSumAC, sumAC, sumBPI0, sumBPIN, sumCG, sumCLATS
983 REAL :: FIN(1), FOUT(1)
985 CALL strace (ient,
'W3XYPFSN')
987 #ifdef W3_DEBUGSOLVER
988 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 1'
993 ith = 1 + mod(isp-1,nth)
999 iobpdr(:)=(1-iobp_loc(:))*(1-iobpd_loc(ith,:))
1002 #ifdef W3_DEBUGSOLVER
1003 WRITE(740+
iaproc,*)
'NX=', nx
1004 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 2'
1015 lambda(1) =
onesixth *(c(i1,1)+c(i2,1)+c(i3,1))
1016 lambda(2) =
onesixth *(c(i1,2)+c(i2,2)+c(i3,2))
1021 nm(ie) = - 1.d0/min(-
thr,sum(min(
zero,ktmp)))
1022 kelem(:,ie) = max(
zero,ktmp)
1029 fl111 = 2.d0*fl11+fl12
1030 fl112 = 2.d0*fl12+fl11
1031 fl211 = 2.d0*fl21+fl22
1032 fl212 = 2.d0*fl22+fl21
1033 fl311 = 2.d0*fl31+fl32
1034 fl312 = 2.d0*fl32+fl31
1035 flall(1,ie) = (fl311 + fl212) *
onesixth + kelem(1,ie)
1036 flall(2,ie) = (fl111 + fl312) *
onesixth + kelem(2,ie)
1037 flall(3,ie) = (fl211 + fl112) *
onesixth + kelem(3,ie)
1039 #ifdef W3_DEBUGSOLVER
1040 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 3'
1047 kksum(ni) = kksum(ni) + kelem(:,ie)
1052 IF (iobp_loc(ip) .EQ. 1 .OR. fsbccfl)
THEN
1053 dtmaxexp =
pdlib_si(ip)/max(dble(10.e-10),kksum(ip)*iobdp_loc(ip))
1054 dtmax = min( dtmax, dtmaxexp)
1056 cflxymax(ip) = max(cflxymax(ip),dble(dt)/dtmaxexp)
1061 cflxy = dble(dt)/dtmax_gl
1062 rest = abs(mod(cflxy,1.0d0))
1063 IF (rest .LT.
thr)
THEN
1064 iter(ik,ith) = abs(nint(cflxy))
1065 ELSE IF (rest .GT.
thr .AND. rest .LT. 0.5d0)
THEN
1066 iter(ik,ith) = abs(nint(cflxy)) + 1
1068 iter(ik,ith) = abs(nint(cflxy))
1072 #ifdef W3_DEBUGSOLVER
1073 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 4'
1079 #ifdef W3_DEBUGSOLVER
1080 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 4.1'
1083 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 4.2'
1086 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 5'
1087 WRITE(740+
iaproc,*)
'IK=', ik,
' ITH=', ith
1093 DO it = 1,
iter(ik,ith)
1094 #ifdef W3_DEBUGSOLVER
1095 WRITE(740+
iaproc,*)
'IK=', ik,
' ITH=', ith
1096 WRITE(740+
iaproc,*)
'IT=', it,
' ITER=',
iter(ik,ith)
1099 WRITE(740+
iaproc,*)
'IT=', it
1107 utilde = nm(ie) * (dot_product(flall(:,ie),u(ni)))
1108 st(ni) = st(ni) + kelem(:,ie) * (u(ni) - utilde)
1110 #ifdef W3_DEBUGSOLVER
1120 u(ip) = max(
zero,u(ip)-dtsi(ip)*st(ip)*(1-iobpa_loc(ip)))*dble(iobpd_loc(ith,ip))*iobdp_loc(ip)
1122 IF (
refpars(3).LT.0.5.AND.iobpd_loc(ith,ip).EQ.0.AND.iobpa_loc(ip).EQ.0) u(ip) = ac(ip)
1125 #ifdef W3_DEBUGSOLVER
1132 #ifdef W3_DEBUGSOLVER
1142 rd1=rd10 - dt * real(
iter(ik,ith)-it)/real(
iter(ik,ith))
1144 IF ( rd2 .GT. 0.001 )
THEN
1145 rd2 = min(1.,max(0.,rd1/rd2))
1151 #ifdef W3_DEBUGSOLVER
1159 ip_glob = mapsf(isbpi(ibi),1)
1162 ac(jx) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) ) &
1163 / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
1164 #ifdef W3_DEBUGSOLVER
1165 sumac=sumac + ac(jx)
1166 sumbpi0=sumbpi0 + bbpi0(isp,ibi)
1167 sumbpin=sumbpin + bbpin(isp,ibi)
1168 sumcg=sumcg + cg(ik,isbpi(ibi))
1169 sumclats=sumclats + clats(isbpi(ibi))
1175 #ifdef W3_DEBUGSOLVER
1176 WRITE(740+
iaproc,*)
'NBI=', nbi
1177 WRITE(740+
iaproc,*)
'RD1=', rd1,
' RD2=', rd2
1178 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumAC=', sumac
1179 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPI0=', sumbpi0
1180 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPIN=', sumbpin
1181 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCG=', sumcg
1182 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCLATS=', sumclats
1187 #ifdef W3_DEBUGSOLVER
1193 #ifdef W3_DEBUGSOLVER
1194 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 6'
1262 USE mpi,
only : mpi_min
1267 INTEGER,
INTENT(IN) :: ISP
1269 REAL,
INTENT(IN) :: DT
1272 REAL,
INTENT(IN) :: C(npa,2)
1274 REAL,
INTENT(INOUT) :: AC(npa)
1276 REAL,
INTENT(IN) :: RD10, RD20
1279 LOGICAL,
INTENT(IN) :: LCALC
1282 INTEGER,
SAVE :: IENT = 0
1285 INTEGER(KIND=1) :: IOBPDR(NX)
1287 INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
1288 INTEGER :: IBI, NI(3), JX
1289 INTEGER :: ISPROC, IP_glob, JSEA, ierr
1293 REAL :: FL1, FL2, FL3
1295 REAL :: FL11, FL12, FL21, FL22, FL31, FL32
1296 REAL :: FL111, FL112, FL211, FL212, FL311, FL312
1297 REAL :: DTSI(npa), U(npa)
1298 REAL :: DTMAX, DTMAX_GL, DTMAXEXP, REST
1299 REAL :: LAMBDA(2), KTMP(3), TMP(3)
1300 REAL :: THETA_L(3), BET1(3), BETAHAT(3)
1301 REAL :: KELEM(3,NE), FLALL(3,NE)
1302 REAL :: KKSUM(npa), ST(npa)
1303 REAL :: NM(NE), FIN(1), FOUT(1)
1305 CALL strace (ient,
'W3XYPFSN')
1307 #ifdef W3_DEBUGSOLVER
1308 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 1'
1313 ith = 1 + mod(isp-1,nth)
1314 ik = 1 + (isp-1)/nth
1315 dtmax = dble(10.e10)
1318 iobpdr(:)=(1-iobp_loc(:))*(1-iobpd_loc(ith,:))
1321 #ifdef W3_DEBUGSOLVER
1322 WRITE(740+
iaproc,*)
'NX=', nx
1323 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 2'
1335 lambda(1) =
onesixth *(c(i1,1)+c(i2,1)+c(i3,1))
1336 lambda(2) =
onesixth *(c(i1,2)+c(i2,2)+c(i3,2))
1341 nm(ie) = - 1.d0/min(-
thr,sum(min(
zero,ktmp)))
1342 kelem(:,ie) = max(
zero,ktmp)
1349 fl111 = 2.d0*fl11+fl12
1350 fl112 = 2.d0*fl12+fl11
1351 fl211 = 2.d0*fl21+fl22
1352 fl212 = 2.d0*fl22+fl21
1353 fl311 = 2.d0*fl31+fl32
1354 fl312 = 2.d0*fl32+fl31
1355 flall(1,ie) = (fl311 + fl212)
1356 flall(2,ie) = (fl111 + fl312)
1357 flall(3,ie) = (fl211 + fl112)
1364 kksum(ni) = kksum(ni) + kelem(:,ie)
1369 IF (iobp_loc(ip) .EQ. 1 .OR. fsbccfl)
THEN
1370 dtmaxexp =
pdlib_si(ip)/max(dble(10.e-10),kksum(ip)*iobdp_loc(ip))
1371 dtmax = min( dtmax, dtmaxexp)
1373 cflxymax(ip) = max(cflxymax(ip),dble(dt)/dtmaxexp)
1378 cflxy = dble(dt)/dtmax_gl
1379 rest = abs(mod(cflxy,1.0d0))
1380 IF (rest .LT.
thr)
THEN
1381 iter(ik,ith) = abs(nint(cflxy))
1382 ELSE IF (rest .GT.
thr .AND. rest .LT. 0.5d0)
THEN
1383 iter(ik,ith) = abs(nint(cflxy)) + 1
1385 iter(ik,ith) = abs(nint(cflxy))
1393 DO it = 1,
iter(ik,ith)
1400 ft = -
onesixth*dot_product(u(ni),flall(:,ie))
1401 utilde = nm(ie) * ( dot_product(kelem(:,ie),u(ni)) - ft )
1402 theta_l(:) = kelem(:,ie) * (u(ni) - utilde)
1403 IF (abs(ft) .GT. 0.0d0)
THEN
1404 bet1(:) = theta_l(:)/ft
1405 IF (any( bet1 .LT. 0.0d0) )
THEN
1406 betahat(1) = bet1(1) + 0.5d0 * bet1(2)
1407 betahat(2) = bet1(2) + 0.5d0 * bet1(3)
1408 betahat(3) = bet1(3) + 0.5d0 * bet1(1)
1409 bet1(1) = max(
zero,min(betahat(1),1.d0-betahat(2),1.d0))
1410 bet1(2) = max(
zero,min(betahat(2),1.d0-betahat(3),1.d0))
1411 bet1(3) = max(
zero,min(betahat(3),1.d0-betahat(1),1.d0))
1412 theta_l(:) = ft * bet1
1417 st(ni) = st(ni) + theta_l
1420 #ifdef W3_DEBUGSOLVER
1427 u(ip) = max(
zero,u(ip)-dtsi(ip)*st(ip)*(1-iobpa_loc(ip)))*iobpd_loc(ith,ip)*iobdp_loc(ip)
1429 IF (
refpars(3).LT.0.5.AND.iobpd_loc(ith,ip).EQ.0.AND.iobpa_loc(ip).EQ.0) u(ip) = ac(ip)
1437 rd1=rd10 - dt * real(
iter(ik,ith)-it)/real(
iter(ik,ith))
1439 IF ( rd2 .GT. 0.001 )
THEN
1440 rd2 = min(1.,max(0.,rd1/rd2))
1453 ip_glob = mapsf(isbpi(ibi),1)
1456 ac(jx) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) ) &
1457 / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
1458 #ifdef W3_DEBUGSOLVER
1459 sumac=sumac + ac(jx)
1460 sumbpi0=sumbpi0 + bbpi0(isp,ibi)
1461 sumbpin=sumbpin + bbpin(isp,ibi)
1462 sumcg=sumcg + cg(ik,isbpi(ibi))
1463 sumclats=sumclats + clats(isbpi(ibi))
1469 #ifdef W3_DEBUGSOLVER
1470 WRITE(740+
iaproc,*)
'NBI=', nbi
1471 WRITE(740+
iaproc,*)
'RD1=', rd1,
' RD2=', rd2
1472 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumAC=', sumac
1473 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPI0=', sumbpi0
1474 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPIN=', sumbpin
1475 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCG=', sumcg
1476 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCLATS=', sumclats
1481 #ifdef W3_DEBUGSOLVER
1488 #ifdef W3_DEBUGSOLVER
1489 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 6'
1558 USE mpi,
only : mpi_min
1564 INTEGER,
INTENT(IN) :: ISP
1566 REAL,
INTENT(IN) :: DT
1569 REAL,
INTENT(IN) :: C(npa,2)
1571 REAL,
INTENT(INOUT) :: AC(npa)
1573 REAL,
INTENT(IN) :: RD10, RD20
1576 LOGICAL,
INTENT(IN) :: LCALC
1579 INTEGER,
SAVE :: IENT = 0
1582 INTEGER(KIND=1) :: IOBPDR(NX)
1584 INTEGER :: IP, IE, POS, IT, I1, I2, I3, I, J, ITH, IK
1585 INTEGER :: IBI, NI(3)
1594 REAL :: SUMTHETA, CFLXY
1595 real*8 :: ft, utilde
1596 real*8 :: fl11, fl12, fl21, fl22, fl31, fl32
1597 real*8 :: fl111, fl112, fl211, fl212, fl311, fl312
1598 REAL :: DTSI(npa), U(npa), UL(npa)
1599 REAL :: DTMAX_GL, DTMAX, DTMAXEXP, REST
1600 real*8 :: lambda(2), ktmp(3)
1601 real*8 :: kelem(3,
ne), flall(3,
ne)
1602 real*8 :: kksum(npa), st(npa)
1603 real*8 ::
nm(
ne), bet1(3), betahat(3), tmp(3), tmp1
1604 INTEGER :: ISPROC, JSEA, IP_glob, ierr, IX
1605 REAL :: eSumAC, sumAC, sumBPI0, sumBPIN, sumCG, sumCLATS
1606 LOGICAL :: testWrite
1607 REAL :: FIN(1), FOUT(1)
1608 REAL :: UIP(NE), UIPIP(NPA), UIMIP(NPA), U3(3)
1609 real*8 :: theta_h(3), theta_ace(3,ne), theta_l(3,ne)
1610 real*8 :: pm(npa), pp(npa), uim(ne), wii(2,npa)
1611 REAL :: USTARI(2,NPA)
1614 CALL strace (ient,
'W3XYPFSN')
1616 #ifdef W3_DEBUGSOLVER
1617 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 1'
1622 ith = 1 + mod(isp-1,nth)
1623 ik = 1 + (isp-1)/nth
1624 dtmax = dble(10.e10)
1627 iobpdr(:)=(1-iobp_loc(:))*(1-iobpd_loc(ith,:))
1630 #ifdef W3_DEBUGSOLVER
1631 WRITE(740+
iaproc,*)
'NX=', nx
1632 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 2'
1644 lambda(1) =
onesixth *(c(i1,1)+c(i2,1)+c(i3,1))
1645 lambda(2) =
onesixth *(c(i1,2)+c(i2,2)+c(i3,2))
1650 nm(ie) = - 1.d0/min(-
thr,sum(min(
zero,ktmp)))
1651 kelem(:,ie) = max(
zero,ktmp)
1658 fl111 = 2.d0*fl11+fl12
1659 fl112 = 2.d0*fl12+fl11
1660 fl211 = 2.d0*fl21+fl22
1661 fl212 = 2.d0*fl22+fl21
1662 fl311 = 2.d0*fl31+fl32
1663 fl312 = 2.d0*fl32+fl31
1664 flall(1,ie) = (fl311 + fl212)
1665 flall(2,ie) = (fl111 + fl312)
1666 flall(3,ie) = (fl211 + fl112)
1673 kksum(ni) = kksum(ni) + kelem(:,ie)
1678 IF (iobp_loc(ip) .EQ. 1 .OR. fsbccfl)
THEN
1679 dtmaxexp =
pdlib_si(ip)/max(dble(10.e-10),kksum(ip)*iobdp_loc(ip))
1680 dtmax = min( dtmax, dtmaxexp)
1682 cflxymax(ip) = max(cflxymax(ip),dble(dt)/dtmaxexp)
1687 cflxy = dble(dt)/dtmax_gl
1688 rest = abs(mod(cflxy,1.0d0))
1689 IF (rest .LT.
thr)
THEN
1690 iter(ik,ith) = abs(nint(cflxy))
1691 ELSE IF (rest .GT.
thr .AND. rest .LT. 0.5d0)
THEN
1692 iter(ik,ith) = abs(nint(cflxy)) + 1
1694 iter(ik,ith) = abs(nint(cflxy))
1702 DO it = 1,
iter(ik,ith)
1710 ft = -
onesixth*dot_product(u(ni),flall(:,ie))
1711 utilde =
nm(ie) * ( dot_product(kelem(:,ie),u(ni)) - ft )
1712 theta_l(:,ie) = kelem(:,ie) * (u(ni) - utilde)
1713 IF (abs(ft) .GT. 0.0d0)
THEN
1714 bet1(:) = theta_l(:,ie)/ft
1715 IF (any( bet1 .LT. 0.0d0) )
THEN
1716 betahat(1) = bet1(1) + 0.5d0 * bet1(2)
1717 betahat(2) = bet1(2) + 0.5d0 * bet1(3)
1718 betahat(3) = bet1(3) + 0.5d0 * bet1(1)
1719 bet1(1) = max(
zero,min(betahat(1),1.d0-betahat(2),1.d0))
1720 bet1(2) = max(
zero,min(betahat(2),1.d0-betahat(3),1.d0))
1721 bet1(3) = max(
zero,min(betahat(3),1.d0-betahat(1),1.d0))
1722 theta_l(:,ie) = ft * bet1
1725 st(ni) = st(ni) + theta_l(:,ie)
1726 theta_h = (1./3.+dt/(2.*
pdlib_tria(ie)) * kelem(:,ie) ) * ft
1728 theta_ace(:,ie) = theta_h-theta_l(:,ie)
1729 pp(ni) = pp(ni) + max(
zero, -theta_ace(:,ie)) * dtsi(ni)
1730 pm(ni) = pm(ni) + min(
zero, -theta_ace(:,ie)) * dtsi(ni)
1733 #ifdef W3_DEBUGSOLVER
1740 ul(ip) = max(
zero,u(ip)-dtsi(ip)*st(ip)*(1-iobpa_loc(ip)))*dble(iobpd_loc(ith,ip))*iobdp_loc(ip)
1743 #ifdef MPI_PARALL_GRID
1747 ustari(1,:) = max(ul,u)
1748 ustari(2,:) = min(ul,u)
1754 uip(ni) = max(uip(ni), maxval( ustari(1,ni) ))
1755 uim(ni) = min(uim(ni), minval( ustari(2,ni) ))
1758 wii(1,:) = min(1.0d0,(uip-ul)/max(
thr,pp))
1759 wii(2,:) = min(1.0d0,(uim-ul)/min(-
thr,pm))
1766 IF (theta_ace(1,ie) .LT.
zero)
THEN
1771 IF (theta_ace(2,ie) .LT.
zero)
THEN
1776 IF (theta_ace(3,ie) .LT.
zero)
THEN
1782 st(i1) = st(i1) + theta_ace(1,ie) * tmp1
1783 st(i2) = st(i2) + theta_ace(2,ie) * tmp1
1784 st(i3) = st(i3) + theta_ace(3,ie) * tmp1
1788 u(ip) = max(
zero,ul(ip)-dtsi(ip)*st(ip)*(1-iobpa_loc(ip)))*dble(iobpd_loc(ith,ip))*iobdp_loc(ip)
1790 IF (
refpars(3).LT.0.5.AND.iobpd_loc(ith,ip).EQ.0.AND.iobpa_loc(ip).EQ.0) u(ip) = ac(ip)
1796 #ifdef W3_DEBUGSOLVER
1806 rd1=rd10 - dt * real(
iter(ik,ith)-it)/real(
iter(ik,ith))
1808 IF ( rd2 .GT. 0.001 )
THEN
1809 rd2 = min(1.,max(0.,rd1/rd2))
1815 #ifdef W3_DEBUGSOLVER
1823 ip_glob = mapsf(isbpi(ibi),1)
1826 ac(jx) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) ) &
1827 / cg(ik,isbpi(ibi)) * clats(isbpi(ibi))
1828 #ifdef W3_DEBUGSOLVER
1829 sumac=sumac + ac(jx)
1830 sumbpi0=sumbpi0 + bbpi0(isp,ibi)
1831 sumbpin=sumbpin + bbpin(isp,ibi)
1832 sumcg=sumcg + cg(ik,isbpi(ibi))
1833 sumclats=sumclats + clats(isbpi(ibi))
1839 #ifdef W3_DEBUGSOLVER
1840 WRITE(740+
iaproc,*)
'NBI=', nbi
1841 WRITE(740+
iaproc,*)
'RD1=', rd1,
' RD2=', rd2
1842 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumAC=', sumac
1843 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPI0=', sumbpi0
1844 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPIN=', sumbpin
1845 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCG=', sumcg
1846 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCLATS=', sumclats
1851 #ifdef W3_DEBUGSOLVER
1858 #ifdef W3_DEBUGSOLVER
1859 WRITE(740+
iaproc,*)
'PDLIB_W3XYPFSN2, step 6'
1918 CHARACTER(*),
INTENT(in) :: string
1922 WRITE(740+
iaproc,*)
'TEST_MPI_STATUS, at string=', string
1927 WRITE(740+
iaproc,*)
'After status settings'
1937 CALL mpi_send(vcollexp,1,mpi_real, 0, 37,
mpi_comm_wcmp, ierr)
1939 WRITE(740+
iaproc,*)
'Leaving the TEST_MPI_STATUS'
2002 real*8,
INTENT(in) :: v(nseal)
2003 CHARACTER(*),
INTENT(in) :: string
2004 INTEGER,
INTENT(IN) :: maxidx
2005 LOGICAL,
INTENT(in) :: CheckUncovered
2006 LOGICAL,
INTENT(in) :: PrintFullValue
2008 real*8,
allocatable :: vcoll(:)
2009 INTEGER,
allocatable :: Status(:)
2010 real*8,
allocatable :: listval(:)
2011 INTEGER,
allocatable :: ListIdx(:)
2013 REAL CoherencyError, eVal1, eVal2, eErr
2014 INTEGER NSEAL_dist, maxidx_dist
2015 INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
2016 INTEGER nbIncorr, idx
2030 allocate(vcoll(nx), status(nx))
2036 isea=mapfs(1,ip_glob)
2037 vcoll(ip_glob)=v(jsea)
2042 nseal_dist = singv(1)
2043 maxidx_dist = singv(2)
2044 allocate(listval(nseal_dist), listidx(nseal_dist))
2047 DO idx=1,maxidx_dist
2048 ip_glob = listidx(idx)
2049 eval1 = vcoll(ip_glob)
2050 eval2 = listval(idx)
2051 vcoll(ip_glob) = eval2
2052 IF (status(ip_glob) .eq. 1)
THEN
2053 eerr=abs(eval1 - eval2)
2054 coherencyerror = coherencyerror + eerr
2058 deallocate(listval, listidx)
2060 WRITE(740+
iaproc,
'(a,f14.7,f14.7,a,a)')
'sum,coh=', sum(vcoll), coherencyerror,
' ', trim(string)
2064 IF (isea .gt. 0)
THEN
2065 IF (status(ix) .eq. 0)
THEN
2070 IF (checkuncovered)
THEN
2071 IF (nbincorr .gt. 0)
THEN
2072 WRITE(*,*)
' nbIncorr=', nbincorr
2073 WRITE(*,*)
' NX=', nx
2074 WRITE(*,*)
' NSEAL=', nseal
2075 WRITE(*,*)
' npa=',
npa
2079 IF (printfullvalue)
THEN
2080 WRITE(740+
iaproc,*)
'Value of V at nodes'
2082 WRITE(740+
iaproc,*)
'IX=', ix,
' V=', vcoll(ix)
2086 deallocate(vcoll, status)
2090 CALL mpi_send(singv,2,mpi_integer, 0, 360,
mpi_comm_wcmp, ierr)
2091 allocate(listval(nseal), listidx(nseal))
2095 isea=mapfs(1,ip_glob)
2096 listval(jsea) = v(jsea)
2097 listidx(jsea) = ip_glob
2099 CALL mpi_send(listval, nseal, mpi_real8, 0, 370,
mpi_comm_wcmp, ierr)
2100 CALL mpi_send(listidx, nseal, mpi_integer, 0, 430,
mpi_comm_wcmp, ierr)
2101 deallocate(listval, listidx)
2153 real*8,
INTENT(in) :: v(
nseal)
2154 CHARACTER(*),
INTENT(in) :: string
2156 LOGICAL :: CheckUncovered = .false.
2157 LOGICAL :: PrintFullValue = .false.
2210 REAL,
INTENT(in) :: V(NSEAL)
2211 CHARACTER(*),
INTENT(in) :: string
2212 LOGICAL :: CheckUncovered = .false.
2213 LOGICAL :: PrintFullValue = .false.
2271 CHARACTER(*),
INTENT(in) :: string
2272 INTEGER,
INTENT(in) :: choice
2273 REAL :: FIELD(NSPEC,NSEAL)
2274 INTEGER ISPEC, JSEA, maxidx
2275 LOGICAL :: PrintMinISP = .false.
2276 LOGICAL :: LocalizeMaximum = .false.
2279 field(ispec,jsea) =
vaold(ispec,jsea)
2282 IF (choice .eq. 1)
THEN
2343 INTEGER,
INTENT(in) :: IMOD
2344 CHARACTER(*),
INTENT(in) :: string
2345 INTEGER,
INTENT(in) :: choice
2346 REAL :: FIELD(NSPEC,NSEAL)
2347 INTEGER ISPEC, JSEA, IP_glob, maxidx
2348 LOGICAL :: PrintMinISP = .false.
2349 LOGICAL :: LocalizeMaximum = .false.
2350 INTEGER :: TEST_IP = 46
2351 INTEGER :: TEST_ISP = 370
2358 WRITE(740+
iaproc,*)
'Entering ALL_INTEGRAL_PRINT, NSEAL=', nseal
2360 IF (nseal .ne.
npa)
THEN
2361 print *,
'NSEAL=', nseal,
" npa=",
npa
2367 field(ispec,jsea) =
va(ispec,jsea)
2368 IF ((ip_glob .eq. test_ip).and.(ispec .eq. test_isp))
THEN
2369 WRITE(740+
iaproc,*)
'ASS TEST_IP=', test_ip,
' TEST_ISP=', test_isp,
' val=',
va(ispec,jsea)
2373 WRITE(740+
iaproc,*)
'Before call to ALL_FIELD_INTEGRAL'
2374 WRITE(740+
iaproc,*)
'NSPEC=', nspec,
' NX=', nx
2376 IF (choice .eq. 1)
THEN
2382 WRITE(740+
iaproc,*)
'After call to ALL_FIELD_INTEGRAL'
2438 REAL,
INTENT(in) :: FIELD(NSPEC,NSEAL)
2439 CHARACTER(*),
INTENT(in) :: string
2440 LOGICAL :: PrintMinISP = .false.
2441 LOGICAL :: LocalizeMaximum = .false.
2505 CHARACTER(*),
INTENT(in) :: string
2506 INTEGER,
INTENT(in) :: maxidx
2507 REAL,
INTENT(in) :: TheARR(NSPEC, npa)
2508 LOGICAL,
INTENT(in) :: PrintMinISP, LocalizeMaximum
2510 REAL Vcoll(NSPEC,NX), VcollExp(NSPEC*NX), rVect(NSPEC*NX)
2511 REAL CoherencyError_Max, CoherencyError_Sum
2512 REAL eVal1, eVal2, eErr
2513 INTEGER LocateMax_I, LocateMax_ISP
2514 INTEGER rStatus(NX), Status(NX)
2515 INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
2516 REAL :: mval, eVal, eSum
2517 REAL :: TheMax, TheSum, TheNb, TheAvg
2518 REAL :: eFact, Threshold
2520 INTEGER nbIncorr, n_control
2522 INTEGER :: TEST_IP = 46
2523 INTEGER :: TEST_ISP = 370
2530 WRITE(740+
iaproc,*)
'CHECK_ARRAY_INTEGRAL NSEAL=', nseal,
' npa=', npa,
' maxidx=', maxidx
2536 vcollexp(isp+nspec*(ip_glob-1)) = thearr(isp,ip)
2537 IF ((ip_glob .eq. test_ip).and.(isp .eq. test_isp))
THEN
2538 WRITE(740+
iaproc,*)
'TEST_IP=', test_ip,
' TEST_ISP=', test_isp,
' val=', thearr(isp,ip)
2546 coherencyerror_max = 0
2547 coherencyerror_sum = 0
2557 IF (rstatus(i) .eq. 1)
THEN
2559 eval1 = vcollexp(isp+nspec*(i-1))
2560 eval2 = rvect(isp+nspec*(i-1))
2561 IF (status(i) .eq. 1)
THEN
2562 eerr=abs(eval1 - eval2)
2563 coherencyerror_sum = coherencyerror_sum + eerr
2564 IF (eerr .gt. coherencyerror_max)
THEN
2565 coherencyerror_max = eerr
2569 IF (isp .eq. 1)
THEN
2570 n_control = n_control + 1
2573 vcollexp(isp+nspec*(i-1))=eval2
2581 CALL mpi_send(vcollexp,nspec*nx,mpi_real , 0, 37,
mpi_comm_wcmp, ierr)
2582 CALL mpi_send(status ,nx ,mpi_integer, 0, 43,
mpi_comm_wcmp, ierr)
2587 vcoll(isp,i)=vcollexp(isp + nspec*(i-1))
2593 IF (isea .gt. 0)
THEN
2594 IF (status(ix) .eq. 0)
THEN
2599 IF (nbincorr .gt. 0)
THEN
2600 WRITE(*,*)
' nbIncorr=', nbincorr
2601 WRITE(*,*)
' NX=', nx
2602 WRITE(*,*)
' npa=', npa
2605 WRITE(740+
iaproc,*)
'CHECK_ARRAY_INTEGRAL n_control=', n_control
2606 WRITE(740+
iaproc,*)
'ARRAY_NX sum,coh=', sum(vcoll), coherencyerror_sum, trim(string)
2607 WRITE(740+
iaproc,*)
'ARRAY_NX max,loc=', coherencyerror_max,locatemax_i,locatemax_isp, trim(string)
2608 IF (printminisp)
THEN
2613 eval=abs(vcoll(isp, ip))
2615 IF (isfirst.eqv. .true.)
then
2618 IF (eval .lt. mval)
THEN
2624 WRITE(740+
iaproc,*)
'ISP=', isp,
' mval/sum=', mval, esum
2628 IF (localizemaximum)
THEN
2634 eval = abs(vcoll(isp, ip))
2635 thesum = thesum + eval
2637 IF (eval .gt. themax)
THEN
2642 theavg = thesum / thenb
2643 WRITE(740+
iaproc,*)
'TheAvg/TheMax=', theavg, themax
2645 threshold=efact * themax
2648 eval = abs(vcoll(isp, ip))
2649 IF (eval .gt. threshold)
THEN
2650 WRITE(740+
iaproc,*)
'ISP/IP/val=', isp, ip, eval
2708 CHARACTER(*),
INTENT(in) :: string
2709 INTEGER,
INTENT(in) :: maxidx
2710 REAL,
INTENT(in) :: TheARR(NSPEC, npa)
2711 real*8 :: thearr_red(npa)
2717 LOGICAL :: FULL_NSPEC = .true.
2718 LOGICAL :: PrintMinISP = .true.
2719 LOGICAL :: LocalizeMaximum = .true.
2720 LOGICAL :: CheckUncovered = .true.
2721 LOGICAL :: PrintFullValue = .true.
2724 IF (full_nspec)
THEN
2728 thearr_red(ip) = sum(abs(thearr(:,ip)))
2788 LOGICAL,
INTENT(IN) :: LCALC
2789 INTEGER,
INTENT(IN) :: IMOD
2790 REAL,
INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
2791 #ifdef W3_DEBUGSOLVER
2799 WRITE(*,*)
'Error: You need to use with JGS_USE_JACOBI'
2800 stop
'Correct your implicit solver options'
2858 LOGICAL,
INTENT(IN) :: LCALC
2859 INTEGER,
INTENT(IN) :: IMOD
2860 REAL,
INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
2924 CHARACTER(*),
INTENT(in) :: string
2925 REAL TotalSumDMM, eDMM, sumDMM
2926 INTEGER IP, IK, ISEA
2927 WRITE(740+
iaproc,*)
'PRINT_WN_STATISTIC'
2932 edmm =
wn(ik+1,isea) -
wn(ik,isea)
2933 sumdmm=sumdmm + abs(edmm)
2935 IF (isea .eq. 1)
THEN
2936 WRITE(740+
iaproc,*)
'ISEA=', isea
2937 WRITE(740+
iaproc,*)
'sumDMM=', sumdmm
2939 totalsumdmm = totalsumdmm + sumdmm
2941 WRITE(740+
iaproc,*)
'string=', string
2942 WRITE(740+
iaproc,*)
'TotalSumDMM=', totalsumdmm
3009 CHARACTER(*),
INTENT(in) :: eFile
3010 REAL,
INTENT(in) :: TheARR(NSPEC, npa)
3012 REAL Vcoll(NSPEC,NX), VcollExp(NSPEC*NX), rVect(NSPEC*NX)
3013 REAL CoherencyError, eVal1, eVal2, eErr
3014 INTEGER rStatus(NX), Status(NX)
3015 INTEGER JSEA, ISEA, iProc, I, IX, ierr, ISP, IP, IP_glob
3031 vcollexp(isp+nspec*(ip_glob-1))=thearr(isp,ip)
3044 IF (rstatus(i) .eq. 1)
THEN
3046 eval1=vcollexp(isp+nspec*(i-1))
3047 eval2=rvect(isp+nspec*(i-1))
3048 vcollexp(isp+nspec*(i-1))=rvect(isp+nspec*(i-1))
3049 IF (status(i) .eq. 1)
THEN
3050 eerr=abs(eval1 - eval2)
3051 coherencyerror = coherencyerror + eerr
3053 vcollexp(isp+nspec*(i-1))=eval2
3061 CALL mpi_send(vcollexp,nspec*nx,mpi_double , 0, 37,
mpi_comm_wcmp, ierr)
3062 CALL mpi_send(status ,nx ,mpi_integer, 0, 43,
mpi_comm_wcmp, ierr)
3067 vcoll(isp,i)=vcollexp(isp + nspec*(i-1))
3070 OPEN(fhndl,
file=efile)
3072 esum=sum(vcoll(:,ix))
3073 WRITE(fhndl,*)
'IX=', ix,
'eSum=', esum
3135 CHARACTER(*),
INTENT(in) :: string
3136 INTEGER J, IP, JP, I, ISP
3137 REAL TheSum1, TheSum2
3155 IF (jp .ne. ip)
THEN
3162 WRITE(740+
iaproc,
'(a,f14.7,f14.7,a,a)')
'TheSum12=', thesum1, thesum2,
' ', string
3227 REAL,
INTENT(IN) :: A(NTH,NK), CG(NK), WN(NK)
3228 REAL,
INTENT(OUT) :: EMEAN, FMEAN, WNMEAN, AMAX
3231 INTEGER,
SAVE :: IENT = 0
3233 REAL :: EB(NK), EBAND
3235 CALL strace (ient,
'W3SPR0')
3248 eb(ik) = eb(ik) + a(ith,ik)
3249 amax = max( amax , a(ith,ik) )
3256 eb(ik) = eb(ik) *
dden(ik) / cg(ik)
3257 emean = emean + eb(ik)
3258 fmean = fmean + eb(ik) /
sig(ik)
3259 wnmean = wnmean + eb(ik) / sqrt(wn(ik))
3265 eband = eb(nk) /
dden(nk)
3266 emean = emean + eband *
fte
3267 fmean = fmean + eband *
ftf
3268 wnmean = wnmean + eband *
ftwn
3272 fmean =
tpiinv * emean / max( 1.e-7 , fmean )
3273 wnmean = ( emean / max( 1.e-7 , wnmean ) )**2
3276 WRITE (
ndst,9000) emean, fmean, wnmean
3284 9000
FORMAT (
' TEST W3SPR0 : E,F,WN MEAN ',3e10.3)
3384 REAL,
INTENT(in) :: DTG, FACX, FACY, VGX, VGY
3385 INTEGER :: IP, ISP, ISEA, IP_glob
3387 INTEGER :: I, J, ITH, IK, J2
3388 INTEGER :: IE, POS, JSEA
3389 INTEGER :: I1, I2, I3, NI(3)
3399 REAL :: CRFS(3), CXY(3,2)
3400 REAL :: KP(3,NSPEC,NE)
3402 REAL :: K1, eSI, eVS, eVD
3403 REAL :: eVal1, eVal2, eVal3
3404 REAL :: DELTAL(3,NSPEC,NE)
3405 REAL :: NM(NSPEC,NE)
3406 REAL :: TRIA03, SIDT, CCOS, CSIN
3407 REAL :: SPEC(NSPEC), DEPTH
3409 #ifdef W3_DEBUGSOLVER
3410 WRITE(740+
iaproc,*)
'calcARRAY_JACOBI, begin'
3424 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 0')
3432 ith = 1 + mod(is-1,nth)
3434 ccos = facx * ecos(ith)
3435 csin = facy * esin(ith)
3436 cxy(:,1) = ccos *
cg(ik,ni) / clats(ni)
3437 cxy(:,2) = csin *
cg(ik,ni)
3439 cxy(:,1) = cxy(:,1) + facx *
cx(ni)/clats(ni)
3440 cxy(:,2) = cxy(:,2) + facy *
cy(ni)
3443 cxy(:,1) = cxy(:,1) - ccurx*vgx/clats(isea)
3444 cxy(:,2) = cxy(:,2) - ccury*vgy
3452 crfs(1) = - onesixth * (2.0d0 *fl31 + fl32 + fl21 + 2.0d0 * fl22 )
3453 crfs(2) = - onesixth * (2.0d0 *fl32 + 2.0d0 * fl11 + fl12 + fl31 )
3454 crfs(3) = - onesixth * (2.0d0 *fl12 + 2.0d0 * fl21 + fl22 + fl11 )
3455 lambda(1) = onesixth * sum(cxy(:,1))
3456 lambda(2) = onesixth * sum(cxy(:,2))
3460 kp(:,is,ie) = max(
zero,k(:))
3461 deltal(:,is,ie) = crfs(:) - kp(:,is,ie)
3462 km(:) = min(
zero,k(:))
3463 nm(is,ie) = 1.d0/min(-thr,sum(km))
3481 ith = 1 + mod(isp-1,nth)
3482 ik = 1 + (isp-1)/nth
3485 eiobpdr=(1-iobp_loc(ip_glob))*(1-iobpd_loc(ith,ip_glob))
3486 IF (eiobpdr .eq. 1)
THEN
3491 dtk = k1 * dtg * iobdp_loc(ip) * (1-iobpa_loc(ip)) * iobpd_loc(ith,ip)
3492 b_jac(isp,ip) =
b_jac(isp,ip) + tria03 *
va(isp,ip) * iobdp_loc(ip) * (1-iobpa_loc(ip)) * iobpd_loc(ith,ip)
3493 tmp3 = dtk * nm(isp,ie)
3504 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 1')
3505 #ifdef W3_DEBUGSOLVER
3506 WRITE(740+
iaproc,*)
'sum(VA)=', sum(
va)
3608 REAL,
INTENT(in) :: DTG, FACX, FACY, VGX, VGY
3609 INTEGER :: IP, ISP, ISEA, IP_glob
3611 INTEGER :: I, J, ITH, IK, J2
3612 INTEGER :: IE, POS, JSEA
3613 INTEGER :: I1, I2, I3, NI(3)
3614 INTEGER :: counter, IB1, IB2, IBR
3616 REAL :: LAMBDA(2), CXYY(2,3), CXY(2,NPA)
3620 REAL :: CRFS(3), K(3)
3622 REAL :: KM(3), DELTAL(3,NE)
3623 REAL :: K1, eSI, eVS, eVD
3624 REAL :: eVal1, eVal2, eVal3
3626 REAL :: TRIA03, SIDT, CCOS, CSIN
3627 REAL :: SPEC(NSPEC), DEPTH, CCOSA(NTH), CSINA(NTH)
3628 INTEGER :: IOBPTH1(NTH), IOBPTH2(NTH)
3630 #ifdef W3_DEBUGSOLVER
3631 WRITE(740+
iaproc,*)
'calcARRAY_JACOBI, begin'
3645 ccosa = facx * ecos(1:nth)
3646 csina = facx * esin(1:nth)
3647 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 0')
3651 ith = 1 + mod(isp-1,nth)
3652 ik = 1 + (isp-1)/nth
3662 cg1 =
cg(ik,ip_glob)
3664 cxy(1,ip) = ccos * cg1/clats(ip_glob)
3665 cxy(2,ip) = csin * cg1
3667 cxy(1,ip) = cxy(1,ip) + facx *
cx(ip_glob)/clats(ip_glob)*iobdp_loc(ip)
3668 cxy(2,ip) = cxy(2,ip) + facy *
cy(ip_glob)*iobdp_loc(ip)
3671 cxy(1,ip) = cxy(1,ip) - ccurx*vgx/clats(isea)
3672 cxy(2,ip) = cxy(2,ip) - ccury*vgy
3678 cxyy(1,:) = cxy(1,ni)
3679 cxyy(2,:) = cxy(2,ni)
3686 crfs(1) = - onesixth * (2.0d0 *fl31 + fl32 + fl21 + 2.0d0 * fl22 )
3687 crfs(2) = - onesixth * (2.0d0 *fl32 + 2.0d0 * fl11 + fl12 + fl31 )
3688 crfs(3) = - onesixth * (2.0d0 *fl12 + 2.0d0 * fl21 + fl22 + fl11 )
3689 lambda(1) = onesixth * sum(cxyy(1,:))
3690 lambda(2) = onesixth * sum(cxyy(2,:))
3694 kp(1:3,ie) = max(
zero,k(1:3))
3695 deltal(1:3,ie) = (crfs(1:3) - kp(1:3,ie)) * 1.d0/min(-thr,sum(min(
zero,k(1:3))))
3700 ib1 = (1-iobpa_loc(ip)) * iobpd_loc(ith,ip)
3701 ib2 = iobpd_loc(ith,ip)
3703 ibr = (1-iobp_loc(ip)) * (1-iobpd_loc(ith,ip)) * (1-iobpa_loc(ip))
3705 IF (iobdp_loc(ip) .eq. 1)
THEN
3716 dtk = kp(pos,ie) * dtg
3719 dtk = kp(pos,ie) * dtg * ib1
3723 dtk = kp(pos,ie) * dtg * ib1
3750 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 1')
3751 #ifdef W3_DEBUGSOLVER
3752 WRITE(740+
iaproc,*)
'sum(VA)=', sum(
va)
3853 REAL,
INTENT(in) :: DTG, FACX, FACY, VGX, VGY
3854 INTEGER :: IP, ISP, ISEA, IP_glob
3856 INTEGER :: I, J, ITH, IK, J2
3857 INTEGER :: IE, POS, JSEA
3858 INTEGER :: I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
3863 INTEGER :: IP1, IP2, IPP1, IPP2
3869 REAL :: CRFS(3), K(3)
3871 REAL :: KM(3), CXY(3,2)
3872 REAL :: K1, eSI, eVS, eVD
3873 REAL :: eVal1, eVal2, eVal3
3875 REAL :: NM, TRIA03, SIDT
3876 REAL :: IEN_LOCAL(6), CG2(NK,NTH)
3878 REAL :: SPEC(NSPEC), DEPTH
3882 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 0')
3887 isea=mapfs(1,ip_glob)
3902 ni_isea = mapfs(1,ni_glob)
3904 ith = 1 + mod(isp-1,nth)
3905 ik = 1 + (isp-1)/nth
3906 ccos = facx * ecos(ith)
3907 csin = facy * esin(ith)
3908 cxy(:,1) = ccos *
cg(ik,ni_isea) / clats(ni_isea)
3909 cxy(:,2) = csin *
cg(ik,ni_isea)
3911 cxy(:,1) = cxy(:,1) + facx *
cx(ni_isea)/clats(ni_isea)
3912 cxy(:,2) = cxy(:,2) + facy *
cy(ni_isea)
3915 cxy(:,1) = cxy(:,1) - ccurx*vgx/clats(isea)
3916 cxy(:,2) = cxy(:,2) - ccury*vgy
3918 fl11 = cxy(2,1)*ien_local(1)+cxy(2,2)*ien_local(2)
3919 fl12 = cxy(3,1)*ien_local(1)+cxy(3,2)*ien_local(2)
3920 fl21 = cxy(3,1)*ien_local(3)+cxy(3,2)*ien_local(4)
3921 fl22 = cxy(1,1)*ien_local(3)+cxy(1,2)*ien_local(4)
3922 fl31 = cxy(1,1)*ien_local(5)+cxy(1,2)*ien_local(6)
3923 fl32 = cxy(2,1)*ien_local(5)+cxy(2,2)*ien_local(6)
3924 crfs(1) = -
onesixth * (2.0d0 *fl31 + fl32 + fl21 + 2.0d0 * fl22 )
3925 crfs(2) = -
onesixth * (2.0d0 *fl32 + 2.0d0 * fl11 + fl12 + fl31 )
3926 crfs(3) = -
onesixth * (2.0d0 *fl12 + 2.0d0 * fl21 + fl22 + fl11 )
3927 lambda(1) =
onesixth * sum(cxy(:,1))
3928 lambda(2) =
onesixth * sum(cxy(:,2))
3929 k(1) = lambda(1) * ien_local(1) + lambda(2) * ien_local(2)
3930 k(2) = lambda(1) * ien_local(3) + lambda(2) * ien_local(4)
3931 k(3) = lambda(1) * ien_local(5) + lambda(2) * ien_local(6)
3932 kp(:) = max(
zero,k(:))
3933 deltal(:) = crfs(:) - kp(:)
3934 km(:) = min(
zero,k(:))
3935 nm = 1.d0/min(-
thr,sum(km))
3938 eiobpdr=(1-iobp_loc(ip))*(1-iobpd_loc(ith,ip))
3939 IF (eiobpdr .eq. 1)
THEN
3944 dtk = k1 * dtg * iobdp_loc(ip) * iobpd_loc(ith,ip) * (1-iobpa_loc(ip))
3953 b_jac(isp,ip) =
b_jac(isp,ip) + tria03 *
va(isp,ip) * iobdp_loc(ip) * iobpd_loc(ith,ip)
3957 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_JACOBI SECTION 1')
3963 SUBROUTINE calcarray_jacobi3(IP,J,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
4056 INTEGER,
INTENT(IN) :: IP
4057 INTEGER,
INTENT(INOUT) :: J
4058 REAL,
INTENT(in) :: DTG, FACX, FACY, VGX, VGY
4059 REAL,
INTENT(out) :: ASPAR_DIAG_LOCAL(NSPEC), B_JAC_LOCAL(NSPEC), ASPAR_OFF_DIAG_LOCAL(NSPEC)
4060 INTEGER :: ISP, ISEA, IP_glob, IPP1, IPP2
4061 INTEGER :: idx, IS, IP1, IP2
4062 INTEGER :: I, ITH, IK, J2
4063 INTEGER :: IE, POS, JSEA
4064 INTEGER :: I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
4074 REAL :: CRFS(3), K(3)
4076 REAL :: KM(3), CXY(3,2)
4077 REAL :: K1, eSI, eVS, eVD
4078 REAL :: eVal1, eVal2, eVal3
4079 REAL :: ien_local(6)
4082 REAL :: TRIA03, SIDT, CCOS, CSIN
4085 aspar_diag_local = 0.d0
4087 aspar_off_diag_local = 0.d0
4104 ni_isea = mapfs(1,ni_glob)
4107 ith = 1 + mod(isp-1,nth)
4108 ik = 1 + (isp-1)/nth
4109 ccos = facx * ecos(ith)
4110 csin = facy * esin(ith)
4111 cxy(:,1) = ccos *
cg(ik,ni_isea) / clats(ni_isea)
4112 cxy(:,2) = csin *
cg(ik,ni_isea)
4114 cxy(:,1) = cxy(:,1) + facx *
cx(ni_isea)/clats(ni_isea)
4115 cxy(:,2) = cxy(:,2) + facy *
cy(ni_isea)
4119 cxy(:,1) = cxy(:,1) - ccurx*vgx/clats(isea)
4120 cxy(:,2) = cxy(:,2) - ccury*vgy
4122 fl11 = cxy(2,1)*ien_local(1)+cxy(2,2)*ien_local(2)
4123 fl12 = cxy(3,1)*ien_local(1)+cxy(3,2)*ien_local(2)
4124 fl21 = cxy(3,1)*ien_local(3)+cxy(3,2)*ien_local(4)
4125 fl22 = cxy(1,1)*ien_local(3)+cxy(1,2)*ien_local(4)
4126 fl31 = cxy(1,1)*ien_local(5)+cxy(1,2)*ien_local(6)
4127 fl32 = cxy(2,1)*ien_local(5)+cxy(2,2)*ien_local(6)
4128 crfs(1) = -
onesixth * (2.0d0 *fl31 + fl32 + fl21 + 2.0d0 * fl22 )
4129 crfs(2) = -
onesixth * (2.0d0 *fl32 + 2.0d0 * fl11 + fl12 + fl31 )
4130 crfs(3) = -
onesixth * (2.0d0 *fl12 + 2.0d0 * fl21 + fl22 + fl11 )
4131 lambda(1) =
onesixth * sum(cxy(:,1))
4132 lambda(2) =
onesixth * sum(cxy(:,2))
4133 k(1) = lambda(1) * ien_local(1) + lambda(2) * ien_local(2)
4134 k(2) = lambda(1) * ien_local(3) + lambda(2) * ien_local(4)
4135 k(3) = lambda(1) * ien_local(5) + lambda(2) * ien_local(6)
4136 kp(:) = max(
zero,k(:))
4137 deltal(:) = crfs(:) - kp(:)
4138 km(:) = min(
zero,k(:))
4139 nm = 1.d0/min(-
thr,sum(km))
4141 eiobpdr=(1-iobp_loc(ip))*(1-iobpd_loc(ith,ip))
4142 IF (eiobpdr .eq. 1)
THEN
4147 dtk = kp(pos) * dble(dtg) * iobdp_loc(ip) * iobpd_loc(ith,ip) * (1-iobpa_loc(ip))
4150 aspar_diag_local(isp) = aspar_diag_local(isp) + tria03 + dtk - tmp3*deltal(pos)
4151 aspar_off_diag_local(isp) = aspar_off_diag_local(isp) - tmp3*deltal(ipp1)*
va(isp,ip1)
4152 aspar_off_diag_local(isp) = aspar_off_diag_local(isp) - tmp3*deltal(ipp2)*
va(isp,ip2)
4154 aspar_diag_local(isp) = aspar_diag_local(isp) + tria03
4156 b_jac_local(isp) = b_jac_local(isp) + tria03 *
vaold(isp,ip) * iobdp_loc(ip) * iobpd_loc(ith,ip)
4164 SUBROUTINE calcarray_jacobi4(IP,DTG,FACX,FACY,VGX,VGY,ASPAR_DIAG_LOCAL,ASPAR_OFF_DIAG_LOCAL,B_JAC_LOCAL)
4256 INTEGER,
INTENT(IN) :: IP
4257 REAL,
INTENT(in) :: DTG, FACX, FACY, VGX, VGY
4258 REAL,
INTENT(out) :: ASPAR_DIAG_LOCAL(NSPEC), B_JAC_LOCAL(NSPEC), ASPAR_OFF_DIAG_LOCAL(NSPEC)
4262 INTEGER :: IE, POS, JSEA
4263 INTEGER :: I, I1, I2, I3, NI(3), NI_GLOB(3), NI_ISEA(3)
4264 INTEGER :: ISP, IP_glob, IPP1, IPP2, IOBPTH1(NTH), IOBPTH2(NTH)
4269 REAL :: DTK, TMP3, D1, D2
4271 REAL :: CRFS(3), K(3)
4272 REAL :: KP(3), UV_CUR(3,2)
4273 REAL :: KM(3), CSX(3), CSY(3)
4274 REAL :: K1, eSI, eVS, eVD
4275 REAL :: eVal1, eVal2, eVal3
4276 REAL :: ien_local(6)
4277 REAL :: DELTAL(3), K_X(3,NK), K_Y(3,NK), K_U(3)
4278 REAL :: CRFS_X(3,NK), CRFS_Y(3,NK), CRFS_U(3)
4279 REAL :: NM, CGFAK(3,NK), CSINA(NTH), CCOSA(NTH)
4280 REAL :: TRIA03, SIDT, CCOS, CSIN
4281 REAL :: FL11_X, FL12_X, FL21_X, FL22_X, FL31_X, FL32_X
4282 REAL :: FL11_Y, FL12_Y, FL21_Y, FL22_Y, FL31_Y, FL32_Y
4283 REAL :: FL11_U, FL12_U, FL21_U, FL22_U, FL31_U, FL32_U
4286 aspar_diag_local =
zero
4288 aspar_off_diag_local =
zero
4291 ccosa(ith) = facx * ecos(ith)
4292 csina(ith) = facx * esin(ith)
4293 iobpth1(ith) = iobdp_loc(ip) * (1-iobpa_loc(ip)) * iobpd_loc(ith,ip)
4294 iobpth2(ith) = iobdp_loc(ip) * iobpd_loc(ith,ip)
4309 ni_isea = mapfs(1,ni_glob)
4314 uv_cur(1:3,1) = facx *
cx(ni_isea) / clats(ni_isea)
4315 uv_cur(1:3,2) = facy *
cy(ni_isea)
4317 lambda(1) =
onesixth*(uv_cur(1,1)+uv_cur(2,1)+uv_cur(3,1))
4318 lambda(2) =
onesixth*(uv_cur(1,2)+uv_cur(2,2)+uv_cur(3,2))
4320 k_u(1) = lambda(1) * ien_local(1) + lambda(2) * ien_local(2)
4321 k_u(2) = lambda(1) * ien_local(3) + lambda(2) * ien_local(4)
4322 k_u(3) = lambda(1) * ien_local(5) + lambda(2) * ien_local(6)
4324 fl11_u = uv_cur(2,1)*ien_local(1)+uv_cur(2,2)*ien_local(2)
4325 fl12_u = uv_cur(3,1)*ien_local(1)+uv_cur(3,2)*ien_local(2)
4326 fl21_u = uv_cur(3,1)*ien_local(3)+uv_cur(3,2)*ien_local(4)
4327 fl22_u = uv_cur(1,1)*ien_local(3)+uv_cur(1,2)*ien_local(4)
4328 fl31_u = uv_cur(1,1)*ien_local(5)+uv_cur(1,2)*ien_local(6)
4329 fl32_u = uv_cur(2,1)*ien_local(5)+uv_cur(2,2)*ien_local(6)
4331 crfs_u(1) = -
onesixth*(2.d0 *fl31_u + fl32_u + fl21_u + 2.d0 * fl22_u)
4332 crfs_u(2) = -
onesixth*(2.d0 *fl32_u + 2.d0 * fl11_u + fl12_u + fl31_u)
4333 crfs_u(3) = -
onesixth*(2.d0 *fl12_u + 2.d0 * fl21_u + fl22_u + fl11_u)
4338 csx =
cg(ik,ni_isea) / clats(ni_isea)
4339 csy =
cg(ik,ni_isea)
4340 lambda(1) =
onesixth * (csx(1) + csx(2) + csx(3))
4341 lambda(2) =
onesixth * (csy(1) + csy(2) + csy(3))
4342 k_x(1,ik) = lambda(1) * ien_local(1)
4343 k_x(2,ik) = lambda(1) * ien_local(3)
4344 k_x(3,ik) = lambda(1) * ien_local(5)
4345 k_y(1,ik) = lambda(2) * ien_local(2)
4346 k_y(2,ik) = lambda(2) * ien_local(4)
4347 k_y(3,ik) = lambda(2) * ien_local(6)
4348 fl11_x = csx(2) * ien_local(1)
4349 fl12_x = csx(3) * ien_local(1)
4350 fl21_x = csx(3) * ien_local(3)
4351 fl22_x = csx(1) * ien_local(3)
4352 fl31_x = csx(1) * ien_local(5)
4353 fl32_x = csx(2) * ien_local(5)
4354 fl11_y = csy(2) * ien_local(2)
4355 fl12_y = csy(3) * ien_local(2)
4356 fl21_y = csy(3) * ien_local(4)
4357 fl22_y = csy(1) * ien_local(4)
4358 fl31_y = csy(1) * ien_local(6)
4359 fl32_y = csy(2) * ien_local(6)
4360 crfs_x(1,ik) = -
onesixth*(2.d0*fl31_x + fl32_x + fl21_x + 2.d0 * fl22_x)
4361 crfs_x(2,ik) = -
onesixth*(2.d0*fl32_x + 2.d0 * fl11_x + fl12_x + fl31_x)
4362 crfs_x(3,ik) = -
onesixth*(2.d0*fl12_x + 2.d0 * fl21_x + fl22_x + fl11_x)
4363 crfs_y(1,ik) = -
onesixth*(2.d0*fl31_y + fl32_y + fl21_y + 2.d0 * fl22_y)
4364 crfs_y(2,ik) = -
onesixth*(2.d0*fl32_y + 2.d0 * fl11_y + fl12_y + fl31_y)
4365 crfs_y(3,ik) = -
onesixth*(2.d0*fl12_y + 2.d0 * fl21_y + fl22_y + fl11_y)
4369 ith = 1 + mod(isp-1,nth)
4370 ik = 1 + (isp-1)/nth
4371 k(1) = k_x(1,ik) * ccosa(ith) + k_y(1,ik) * csina(ith) + k_u(1)
4372 k(2) = k_x(2,ik) * ccosa(ith) + k_y(2,ik) * csina(ith) + k_u(2)
4373 k(3) = k_x(3,ik) * ccosa(ith) + k_y(3,ik) * csina(ith) + k_u(3)
4374 crfs(1) = crfs_x(1,ik) * ccosa(ith) + crfs_y(1,ik) * csina(ith) + crfs_u(1)
4375 crfs(2) = crfs_x(2,ik) * ccosa(ith) + crfs_y(2,ik) * csina(ith) + crfs_u(2)
4376 crfs(3) = crfs_x(3,ik) * ccosa(ith) + crfs_y(3,ik) * csina(ith) + crfs_u(3)
4378 kp(1:3) = max(
zero,k(1:3))
4379 deltal(1:3) = crfs(1:3) - kp(1:3)
4381 dtk = kp(pos) * dtg * iobpth1(ith)
4382 tmp3 = dtk * 1.d0/min(-
thr,sum(min(
zero,k(1:3))))
4384 aspar_diag_local(isp) = aspar_diag_local(isp) + tria03 + dtk - tmp3*deltal(pos)
4385 d1 = deltal(ipp1)*
va(isp,ip1)
4386 d2 = deltal(ipp2)*
va(isp,ip2)
4387 aspar_off_diag_local(isp) = aspar_off_diag_local(isp) - ( tmp3 * ( d1 + d2 ) )
4390 aspar_diag_local(isp) = aspar_diag_local(isp) + tria03
4392 b_jac_local(isp) = b_jac_local(isp) + tria03 *
vaold(isp,ip) * iobpth2(ith)
4455 REAL,
INTENT(in) :: DTG
4456 INTEGER IP, IP_glob, ITH, IK
4459 REAL :: B_SIG(NSPEC), B_THE(NSPEC)
4460 REAL :: CP_SIG(NSPEC), CM_SIG(NSPEC)
4461 REAL :: CP_THE(NSPEC), CM_THE(NSPEC)
4462 REAL :: CAD(NSPEC), CAS(NSPEC)
4463 REAL :: DMM(0:NK2), eVal
4464 REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
4465 LOGICAL :: DoLimiterRefraction = .false.
4466 LOGICAL :: DoLimiterFreqShit = .false.
4469 LOGICAL :: LSIG = .false.
4475 isea=
mapfs(1,ip_glob)
4477 IF (fsfreqshift .AND. lsig)
THEN
4479 IF (iobp_loc(ip).eq.1.and.iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4481 cp_sig = max(
zero,cas)
4482 cm_sig = min(
zero,cas)
4486 isp=ith + (ik-1)*nth
4487 b_sig(isp)= cp_sig(isp)/dmm(ik-1) - cm_sig(isp)/dmm(ik)
4489 isp = ith + (nk-1)*nth
4490 b_sig(isp)= b_sig(isp) + cm_sig(isp)/dmm(nk) * fachfa
4498 IF (iobp_loc(ip).eq.1.and.iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4500 #ifdef W3_DEBUGFREQSHIFT
4501 WRITE(740+
iaproc,*)
'sum(CWNB_M2)=', sum(cwnb_m2)
4505 isp = ith + (ik-1)*nth
4506 eval = dwni_m2(ik) * ( min(cwnb_m2(isp - nth),
zero) - max(cwnb_m2(isp),
zero) )
4509 eval = dwni_m2(nk) * min(cwnb_m2(ith + (nk-1)*nth),
zero) * fachfa
4522 IF (fsrefraction)
THEN
4523 IF (iobp_loc(ip) .eq. 1 .and. iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4530 #ifdef W3_DEBUGREFRACTION
4531 WRITE(740+
iaproc,*)
'refraction IP=', ip,
' ISEA=', isea
4532 WRITE(740+
iaproc,*)
'sum(abs(CAD))=', sum(abs(cad))
4535 cp_the = dtg*max(
zero,cad)
4536 cm_the = dtg*min(
zero,cad)
4537 b_the(:) = cp_the(:) - cm_the(:)
4604 REAL,
INTENT(in) :: DTG
4605 REAL,
INTENT(inout) :: ASPAR_DIAG_LOCAL(nspec,NSEAL)
4606 INTEGER IP, IP_glob, ITH, IK
4609 REAL :: B_SIG(NSPEC), B_THE(NSPEC)
4610 REAL :: CP_SIG(NSPEC), CM_SIG(NSPEC)
4611 REAL :: CP_THE(NSPEC), CM_THE(NSPEC)
4612 REAL :: CAD(NSPEC), CAS(NSPEC)
4613 REAL :: DMM(0:NK2), eVal
4614 REAL :: DWNI_M2(NK), CWNB_M2(1-NTH:NSPEC)
4615 LOGICAL :: DoLimiterRefraction = .false.
4616 LOGICAL :: DoLimiterFreqShit = .false.
4619 LOGICAL :: LSIG = .false.
4625 isea=
mapfs(1,ip_glob)
4630 IF (fsfreqshift .AND. lsig)
THEN
4632 IF (iobp_loc(ip).eq.1.and.iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4634 cp_sig = max(
zero,cas)
4635 cm_sig = min(
zero,cas)
4639 isp=ith + (ik-1)*nth
4640 b_sig(isp)= cp_sig(isp)/dmm(ik-1) - cm_sig(isp)/dmm(ik)
4642 isp = ith + (nk-1)*nth
4643 b_sig(isp)= b_sig(isp) + cm_sig(isp)/dmm(nk) * fachfa
4645 aspar_diag_local(:,ip) = aspar_diag_local(:,ip) + b_sig * esi
4653 IF (iobp_loc(ip).eq.1)
THEN
4655 #ifdef W3_DEBUGFREQSHIFT
4656 WRITE(740+
iaproc,*)
'sum(CWNB_M2)=', sum(cwnb_m2)
4660 isp = ith + (ik-1)*nth
4661 eval = dwni_m2(ik) * ( min(cwnb_m2(isp - nth),
zero) - max(cwnb_m2(isp),
zero) )
4664 ELSE IF (
imem == 2)
THEN
4665 aspar_diag_local(isp,ip) = aspar_diag_local(isp,ip) - esi * eval
4668 eval = dwni_m2(nk) * min(cwnb_m2(ith + (nk-1)*nth),
zero) * fachfa
4670 aspar_diag_local(ith0 + ith,ip) = aspar_diag_local(ith0 + ith,ip) + esi * eval
4679 IF (fsrefraction)
THEN
4680 IF (iobp_loc(ip) .eq. 1.and.iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4687 #ifdef W3_DEBUGREFRACTION
4688 WRITE(740+
iaproc,*)
'refraction IP=', ip,
' ISEA=', isea
4689 WRITE(740+
iaproc,*)
'sum(abs(CAD))=', sum(abs(cad))
4692 cp_the = dtg*max(
zero,cad)
4693 cm_the = dtg*min(
zero,cad)
4694 b_the(:) = cp_the(:) - cm_the(:)
4695 aspar_diag_local(:,ip) = aspar_diag_local(:,ip) + b_the(:)*esi
4769 REAL,
INTENT(in) :: DTG
4770 REAL,
PARAMETER :: COEF4 = 5.0e-07
4771 REAL,
PARAMETER :: FACDAM = 1
4772 INTEGER JSEA, IP, IP_glob, ISEA
4773 INTEGER IK, ITH, ISP, IS0
4775 REAL :: eSI, eVS, eVD, SIDT
4776 REAL :: DEPTH, DAM(NSPEC), RATIO, MAXDAC, VSDB(NSPEC), VDDB(NSPEC)
4777 REAL :: PreVS, eDam, DVS, FREQ, EMEAN, FMEAN, WNMEAN, AMAX, CG1(NK),WN1(NK),SPEC_VA(NSPEC)
4784 isea = mapfs(1,ip_glob)
4786 IF ((iobp_loc(ip).eq.1..or.iobp_loc(jsea).eq. 3).and.iobdp_loc(ip).eq.1.and.iobpa_loc(ip).eq.0)
THEN
4789 dam(1+(ik-1)*nth) = facp / ( sig(ik) *
wn(ik,isea)**3 )
4794 dam(ith+is0) = dam(1+is0)
4808 isp=ith + (ik-1)*nth
4809 spec_va(isp) =
va(isp,jsea) *
cg(ik,isea) / clats(isea)
4813 SELECT CASE (nint(
sdbsc))
4815 CALL w3sdb1 ( jsea, spec_va, depth, emean, fmean, wnmean, cg1, lbreak, vsdb, vddb )
4827 isp=ith + (ik-1)*nth
4828 spec_va(isp) =
va(isp,jsea) *
cg(ik,isea) / clats(isea)
4832 CALL w3sdb2 ( jsea, spec_va, depth, emean, fmean, cg1, lbreak, vsdb, vddb )
4836 isp=ith + (ik-1)*nth
4838 maxdac = facdam * dam(isp)
4839 thefactor = dtg / max( 1. , (1.-dtg*
vdtot(isp,jsea)))
4840 dvs =
vstot(isp,jsea) * thefactor
4841 dvs = sign(min(maxdac,abs(dvs)),dvs)
4842 prevs = dvs / thefactor
4844 prevs =
vstot(isp,jsea)
4846 evs = prevs * clats(isea) /
cg(ik,isea)
4847 evd = dble(
vdtot(isp,jsea))
4849 evs = evs + dble(vsdb(isp)) /
cg(ik,isea) * clats(isea)
4850 evd = evd + dble(vddb(isp))
4853 evs = evs + dble(vsdb(isp)) /
cg(ik,isea) * clats(isea)
4854 evd = evd + dble(vddb(isp))
4856 b_jac(isp,ip) =
b_jac(isp,ip) + sidt * (evs - evd*
va(isp,jsea))
4931 REAL,
INTENT(in) :: DTG
4932 REAL,
INTENT(inout) :: ASPAR_DIAG_LOCAL(:,:)
4933 REAL,
PARAMETER :: COEF4 = 5.0e-07
4934 REAL,
PARAMETER :: FACDAM = 1
4935 INTEGER JSEA, IP, IP_glob, ISEA
4936 INTEGER IK, ITH, ISP, IS0
4938 REAL :: eSI, eVS, eVD, SIDT
4939 REAL :: DEPTH, DAM(NSPEC), RATIO, MAXDAC, VSDB(NSPEC), VDDB(NSPEC)
4940 REAL :: PreVS, eDam, DVS, FREQ, EMEAN, FMEAN, WNMEAN, AMAX, CG1(NK),WN1(NK),SPEC_VA(NSPEC)
4947 isea = mapfs(1,ip_glob)
4949 IF (iobp(ip_glob).eq.1..and.iobdp(ip_glob).eq.1.and.iobpa(ip_glob).eq.0)
THEN
4951 dam(1+(ik-1)*nth) = facp / ( sig(ik) *
wn(ik,isea)**3 )
4956 dam(ith+is0) = dam(1+is0)
4969 isp=ith + (ik-1)*nth
4970 spec_va(isp) =
va(isp,jsea) *
cg(ik,isea) / clats(isea)
4974 SELECT CASE (nint(
sdbsc))
4976 CALL w3sdb1 ( jsea, spec_va, depth, emean, fmean, wnmean, cg1, lbreak, vsdb, vddb )
4988 isp=ith + (ik-1)*nth
4989 spec_va(isp) =
va(isp,jsea) *
cg(ik,isea) / clats(isea)
4993 CALL w3sdb2 ( jsea, spec_va, depth, emean, fmean, cg1, lbreak, vsdb, vddb )
4997 isp=ith + (ik-1)*nth
4999 maxdac = facdam * dam(isp)
5000 thefactor = dtg / max( 1. , (1.-dtg*
vdtot(isp,jsea)))
5001 dvs =
vstot(isp,jsea) * thefactor
5002 dvs = sign(min(maxdac,abs(dvs)),dvs)
5003 prevs = dvs / thefactor
5005 prevs =
vstot(isp,jsea)
5007 evs = prevs /
cg(ik,isea) * clats(isea)
5008 evd = dble(
vdtot(isp,jsea))
5010 evs = evs + dble(vsdb(isp)) /
cg(ik,isea) * clats(isea)
5011 evd = evd + dble(vddb(isp))
5014 evs = evs + dble(vsdb(isp)) /
cg(ik,isea) * clats(isea)
5015 evd = evd + dble(vddb(isp))
5017 b_jac(isp,ip) =
b_jac(isp,ip) + sidt * (evs - evd*
va(isp,jsea))
5018 aspar_diag_local(isp,ip) = aspar_diag_local(isp,ip) - sidt * evd
5091 INTEGER,
SAVE :: IENT = 0
5096 REAL :: RD1, RD2, RD10, RD20
5098 INTEGER :: IK, ITH, ISEA, JSEA
5101 CALL strace (ient,
'APPLY_BOUNDARY_CONDITION_VA')
5103 IF (gtype .eq. ungtype)
THEN
5114 IF ( rd2 .GT. 0.001 )
THEN
5115 rd2 = min(1.,max(0.,rd1/rd2))
5124 IF (jsea .gt. 0)
THEN
5127 isp=ith + (ik-1)*
nth
5128 eac = ( rd1*
bbpi0(isp,ibi) + rd2*
bbpin(isp,ibi) ) &
5130 eva = max(0.,
cg(ik,isea)/clats(isea)*eac)
5202 INTEGER,
INTENT(IN) :: IMOD
5211 INTEGER,
SAVE :: IENT = 0
5216 #ifdef W3_DEBUGSOLVER
5217 real*8 :: sumac(nspec)
5218 REAL :: sumBPI0(NSPEC), sumBPIN(NSPEC), sumCG, sumCLATS
5221 REAL :: ETOT, HSIG_bound, eVA, eAC, FACTOR
5223 REAL :: RD1, RD2, RD10, RD20
5224 INTEGER :: IK, ITH, ISEA
5225 INTEGER :: IBI, IP_glob, ISP, JX
5227 CALL strace (ient,
'APPLY_BOUNDARY_CONDITION')
5229 #ifdef W3_DEBUGSOLVERCOH
5243 IF ( rd2 .GT. 0.001 )
THEN
5244 rd2 = min(1.,max(0.,rd1/rd2))
5250 #ifdef W3_DEBUGSOLVER
5251 WRITE(740+
iaproc,*)
'Begin of APPLY_BOUNDARY_CONDITION'
5262 ip_glob = mapsf(isea,1)
5267 isp=ith + (ik-1)*nth
5268 va(isp,jx) = (( rd1*
bbpi0(isp,ibi) + rd2*
bbpin(isp,ibi) ) &
5269 /
cg(ik,
isbpi(ibi)) * clats(
isbpi(ibi))) * iobdp_loc(jx)
5276 factor =
dden(ik)/
cg(ik,isea)
5277 isp=ith + (ik-1)*nth
5278 eac=real(
va(isp,jx))
5279 eva=
cg(ik,isea)/clats(isea)*eac
5280 etot = etot + eva*factor
5283 hsig_bound=4.*sqrt(etot)
5284 WRITE(740+
iaproc,*)
'IBI=', ibi,
' HSIG=', hsig_bound
5287 #ifdef W3_DEBUGSOLVER
5288 sumac=sumac +
va(:,jx)
5289 sumbpi0=sumbpi0 +
bbpi0(:,ibi)
5290 sumbpin=sumbpin +
bbpin(:,ibi)
5291 sumcg=sumcg +
cg(ik,
isbpi(ibi))
5292 sumclats=sumclats + clats(
isbpi(ibi))
5296 #ifdef W3_DEBUGSOLVER
5297 WRITE(740+
iaproc,*)
'RD1=', rd1,
' RD2=', rd2
5299 #ifdef W3_DEBUGSOLVERALL
5301 WRITE(740+
iaproc,*)
'RD1=', rd1,
' RD2=', rd2
5302 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumAC=', sumac(isp)
5303 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPI0=', sumbpi0(isp)
5304 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumBPIN=', sumbpin(isp)
5305 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCG=', sumcg
5306 WRITE(740+
iaproc,*)
'ISP=', isp,
'sumCLATS=', sumclats
5309 #ifdef W3_DEBUGSOLVER
5310 WRITE(740+
iaproc,*)
'Begin of APPLY_BOUNDARY_CONDITION'
5313 #ifdef W3_DEBUGSOLVERCOH
5381 INTEGER,
SAVE :: IENT = 0
5386 INTEGER,
INTENT(in) :: IP
5387 REAL,
INTENT(in) :: ACOLD(NSPEC)
5388 REAL,
INTENT(inout) :: ACLOC(NSPEC)
5389 REAL,
INTENT(in) :: DTG
5390 INTEGER :: MELIM = 1
5391 REAL :: LIMFAK = 0.1
5392 REAL :: CONST, SND, eWN, eWK, eWKpow
5393 REAL :: eFact, eSPSIG
5395 REAL :: OLDAC, NEWAC, NEWDAC
5397 REAL :: dac, limac, eDam
5398 INTEGER IP_glob, ISEA
5399 INTEGER :: IK, ITH, ISP
5400 LOGICAL :: LLIMITER_WWM
5402 CALL strace (ient,
'ACTION_LIMITER_LOCAL')
5405 isea=
mapfs(1,ip_glob)
5407 const =
tpi**2*3.0*1.0e-7*dtg*espsig
5408 snd =
tpi*5.6*1.0e-3
5410 llimiter_wwm = .false.
5412 IF (llimiter_wwm)
THEN
5415 IF (melim .eq. 1)
THEN
5420 maxdac = dble(0.0081*limfak/(efact*ewkpow*
cg(ik,isea)))
5423 isp=ith + (ik-1)*
nth
5426 newdac = newac - oldac
5427 newdac = sign(min(maxdac,abs(newdac)), newdac)
5428 newval = max(0., oldac + newdac )
5434 edam=dble(
facp / (
sig(ik) *
wn(ik,isea)**3))
5436 isp = ith + (ik-1)*
nth
5437 dac = acloc(isp) - acold(isp)
5438 limac = sign(min(edam,abs(dac)),dac)
5439 acloc(isp) = max(0., acloc(isp) + limac)
5511 USE mpi,
only : mpi_sum, mpi_int
5540 LOGICAL,
INTENT(IN) :: LCALC
5541 INTEGER,
INTENT(IN) :: IMOD
5542 REAL,
INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
5544 INTEGER :: IP, ISP, ITH, IK, JSEA, ISEA, IP_glob, IS0
5546 INTEGER :: nbIter, ISPnextDir, ISPprevDir
5547 INTEGER :: ISPp1, ISPm1, JP, ICOUNT1, ICOUNT2
5549 real*8 :: ccos, csin, ccurx, ccury
5550 real*8 :: esum(nspec), frlocal
5551 real*8 :: ea_the, ec_the, ea_sig, ec_sig, esi
5552 real*8 :: cad(nspec), cas(nspec), acloc(nspec)
5553 real*8 :: cp_sig(nspec), cm_sig(nspec)
5554 real*8 :: efactm1, efactp1
5555 real*8 :: sum_prev, sum_new, p_is_converged, diffnew, prop_conv
5556 real*8 :: sum_l2, sum_l2_gl
5557 REAL :: DMM(0:NK2), DAM(NSPEC), DAM2(NSPEC), SPEC(NSPEC)
5558 real*8 :: ediff(nspec), eprod(nspec), ediffb(nspec)
5559 real*8 :: dwni_m2(nk), cwnb_m2(1-nth:nspec)
5560 REAL :: VAnew(NSPEC), VFLWN(1-NTH:NSPEC), JAC, JAC2
5561 REAL :: VAAnew(1-NTH:NSPEC+NTH), VAAacloc(1-NTH:NSPEC+NTH)
5562 REAL :: VAinput(NSPEC), VAacloc(NSPEC), ASPAR_DIAG(NSPEC)
5563 REAL :: aspar_diag_local(nspec), aspar_off_diag_local(nspec), b_jac_local(nspec)
5564 real*8 :: ediffsing, esumpart
5565 REAL :: EMEAN, FMEAN, FMEAN1, WNMEAN, AMAX, U10ABS, U10DIR, TAUA, TAUADIR
5566 REAL :: USTAR, USTDIR, TAUWX, TAUWY, CD, Z0, CHARN, FMEANWS, DLWMEAN
5567 real*8 :: eval1, eval2
5568 real*8 :: eva, evo, cg2, newdac, newac, oldac, maxdac
5569 REAL :: CG1(0:NK+1), WN1(0:NK+1)
5570 LOGICAL :: LCONVERGED(NSEAL), lexist, LLWS(NSPEC)
5572 INTEGER :: ipiter(nseal), ipitergl(np_global), ipiterout(np_global)
5575 REAL :: IntDiff, eVA_w3srce, eVAsolve, SumACout
5576 REAL :: SumVAin, SumVAout, SumVAw3srce, SumVS, SumVD, VS_w3srce
5577 REAL :: VAsolve(NSPEC)
5581 #ifdef W3_DEBUGSOLVERCOH
5582 REAL :: TheARR(NSPEC, npa)
5583 REAL :: PRE_VA(NSPEC, npa)
5584 REAL :: OffDIAG(NSPEC, npa)
5585 real*8 :: eoff(nspec)
5586 real*8 :: esum1(nspec), esum2(nspec)
5588 CHARACTER(len=128) eFile
5591 INTEGER is_converged, itmp
5593 INTEGER :: TESTNODE = 923
5595 LOGICAL :: LSIG = .false.
5600 #ifdef W3_DEBUGSOLVERCOH
5603 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 0')
5608 CALL mpi_comm_rank(mpi_comm_wcmp, myrank, ierr)
5610 #ifdef W3_DEBUGSOLVER
5611 WRITE(740+
iaproc,*)
'PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK, begin'
5612 WRITE(740+
iaproc,*)
'NX=', nx
5613 WRITE(740+
iaproc,*)
'NP=', np
5614 WRITE(740+
iaproc,*)
'NPA=', npa
5615 WRITE(740+
iaproc,*)
'NSEA=', nsea
5616 WRITE(740+
iaproc,*)
'NSEAL=', nseal
5617 WRITE(740+
iaproc,*)
'NBI=', nbi
5618 WRITE(740+
iaproc,*)
'B_JGS_TERMINATE_NORM=', b_jgs_terminate_norm
5619 WRITE(740+
iaproc,*)
'B_JGS_TERMINATE_DIFFERENCE=', b_jgs_terminate_difference
5620 WRITE(740+
iaproc,*)
'B_JGS_TERMINATE_MAXITER=', b_jgs_terminate_maxiter
5621 WRITE(740+
iaproc,*)
'B_JGS_MAXITER=', b_jgs_maxiter
5622 WRITE(740+
iaproc,*)
'B_JGS_BLOCK_GAUSS_SEIDEL=', b_jgs_block_gauss_seidel
5623 WRITE(740+
iaproc,*)
'FSREFRACTION=', fsrefraction
5624 WRITE(740+
iaproc,*)
'FSFREQSHIFT=', fsfreqshift
5625 WRITE(740+
iaproc,*)
'B_JGS_LIMITER=', b_jgs_limiter
5626 WRITE(740+
iaproc,*)
'B_JGS_BLOCK_GAUSS_SEIDEL=', b_jgs_block_gauss_seidel
5630 WRITE(740+
iaproc,*)
'optionCall=', optioncall
5633 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 1')
5638 WRITE(740+
iaproc,*)
'NSEAL =', nseal,
'NP =', np,
'NPA =', npa
5640 #ifdef W3_DEBUGSOLVERCOH
5647 isea = mapfs(1,ip_glob)
5649 ith = 1 + mod(isp-1,nth)
5650 ik = 1 + (isp-1)/nth
5652 CALL wavnu_local(sig(ik),dw(isea),wn1(ik),cg1(ik))
5654 cg1(ik) = cg(ik,isea)
5656 va(isp,jsea) =
va(isp,jsea) / cg1(ik) * clats(isea)
5663 WRITE(740+
iaproc,*)
'JSEA=', jsea
5664 WRITE(740+
iaproc,*)
'min/max/sum(VA)=', minval(
va(:,jsea)), maxval(
va(:,jsea)), sum(
va(:,jsea))
5668 #ifdef W3_DEBUGSOLVERCOH
5673 #ifdef W3_DEBUGSOLVER
5675 WRITE(740+
iaproc,*)
'JACOBI_SOLVER, step 4'
5676 WRITE(740+
iaproc,*)
'FSSOURCE=', fssource
5677 WRITE(740+
iaproc,*)
'FSREFRACTION=', fsrefraction
5678 WRITE(740+
iaproc,*)
'FSFREQSHIFT=', fsfreqshift
5680 WRITE(740+
iaproc,*)
'DTG=', dtg
5685 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 2')
5687 IF (.not. lsloc)
THEN
5690 ELSE IF (imem == 2)
THEN
5695 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 3')
5700 IF (.not. lsloc)
THEN
5703 ELSE IF (imem == 2)
THEN
5708 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 4')
5716 #ifdef W3_DEBUGSOLVER
5719 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 5')
5721 #ifdef W3_DEBUGSOLVER
5727 IF (fsfreqshift .or. fsrefraction)
THEN
5730 ELSE IF (imem == 2)
THEN
5735 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION 6')
5737 #ifdef W3_DEBUGSOLVERCOH
5740 thearr(:, ip)=real(
aspar_jac(:, pdlib_i_diag(ip)))
5746 lconverged(ip) = .false.
5756 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 1')
5761 isea = mapfs(1,ip_glob)
5762 IF (iobdp_loc(ip) .eq. 0)
THEN
5763 is_converged = is_converged + 1
5764 lconverged(ip) = .true.
5770 CALL wavnu_local(sig(ik),dw(isea),wn1(ik),cg1(ik))
5772 cg1(ik) = cg(ik,isea)
5773 wn1(ik) = wn(ik,isea)
5778 isea = mapfs(1,ip_glob)
5782 IF (.NOT. lconverged(ip))
THEN
5784 ipiter(ip) = ipiter(ip) + 1
5786 #ifdef W3_DEBUGFREQSHIFT
5787 WRITE(740+
iaproc,*)
'Begin loop'
5788 WRITE(740+
iaproc,*)
'IP/IP_glob/ISEA/JSEA=', ip, ip_glob, isea, jsea
5791 WRITE(740+
iaproc,*)
'IP=', ip,
' IP_glob=', ip_glob
5792 WRITE(740+
iaproc,*)
'sum(VA)in=', sum(
va(:,ip))
5794 #ifdef W3_DEBUGFREQSHIFT
5798 vainput(isp) = dble(cg(ik,isea)/clats(isea)) *
va(isp, ip)
5799 vaacloc(isp) = dble(cg(ik,isea)/clats(isea)) * acloc(isp)
5801 WRITE(740+
iaproc,*)
'sum(VAold/VAinput/VAacloc)=', sum(
vaold), sum(vainput), sum(vaacloc)
5804 sum_prev = sum(acloc)
5807 CALL calcarray_jacobi4(ip,dtg,facx,facy,vgx,vgy,aspar_diag_local,aspar_off_diag_local,b_jac_local)
5808 aspar_diag(1:nspec) = aspar_diag_local(1:nspec) +
aspar_diag_all(1:nspec,ip)
5809 esum = b_jac_local - aspar_off_diag_local +
b_jac(1:nspec,ip)
5810 ELSEIF (imem == 1)
THEN
5811 esum(1:nspec) =
b_jac(1:nspec,ip)
5812 aspar_diag(1:nspec) =
aspar_jac(1:nspec,pdlib_i_diag(ip))
5813 #ifdef W3_DEBUGFREQSHIFT
5814 WRITE(740+
iaproc,*)
'eSI=', esi
5815 WRITE(740+
iaproc,*)
'sum(ASPAR_DIAG)=', sum(aspar_diag)
5818 WRITE(740+
iaproc,*)
'Step 1: sum(eSum)=', sum(esum)
5820 #ifdef W3_DEBUGSOLVERCOH
5823 DO i = pdlib_ia_p(ip)+1, pdlib_ia_p(ip+1)
5825 IF (jp .ne. ip)
THEN
5826 eprod(1:nspec) =
aspar_jac(1:nspec,i) *
va(1:nspec,jp)
5827 esum(1:nspec) = esum(1:nspec) - eprod(1:nspec)
5828 #ifdef W3_DEBUGSOLVERALL
5829 WRITE(740+
iaproc,
'(A20,3I10,20E20.10)')
'OFF DIAGONAL', ip, i, jp, sum(
b_jac(:,ip)), sum(esum), sum(
aspar_jac(:,i)), sum(
va(:,jp))
5831 #ifdef W3_DEBUGSOLVERCOH
5838 #ifdef W3_DEBUGSOLVERCOH
5839 offdiag(:, ip)=real(eoff)
5841 #ifdef W3_DEBUGSOLVERCOHALL
5842 WRITE(740+
iaproc,*)
'Step 2: sum(eSum)=', sum(esum),
' eOff=', sum(eoff)
5844 IF (fsrefraction)
THEN
5845 #ifdef W3_DEBUGREFRACTION
5846 WRITE(740+
iaproc,*)
'Adding refraction terms to eSum'
5850 ispprevdir=listispprevdir(isp)
5851 ispnextdir=listispnextdir(isp)
5852 ea_the = - dtg*esi*max(zero,cad(ispprevdir))
5853 ec_the = dtg*esi*min(zero,cad(ispnextdir))
5854 esum(isp) = esum(isp) - ea_the *
va(ispprevdir,ip)
5855 esum(isp) = esum(isp) - ec_the *
va(ispnextdir,ip)
5859 WRITE(740+
iaproc,*)
'Step 3: sum(eSum)=', sum(esum)
5861 IF (fsfreqshift .and. lsig)
THEN
5864 cp_sig = max(zero,cas)
5865 cm_sig = min(zero,cas)
5867 dmm(ik+1) = dble(wn1(ik+1) - wn1(ik))
5873 isp = ith + (ik -1)*nth
5874 ispm1 = ith + (ik-1 -1)*nth
5875 efactm1 = cg1(ik-1) / cg1(ik)
5876 ea_sig = - esi * cp_sig(ispm1)/dmm(ik-1) * efactm1
5877 esum(isp) = esum(isp) - ea_sig*
va(ispm1,ip)
5880 isp = ith + (ik -1)*nth
5881 ispp1 = ith + (ik+1 -1)*nth
5882 efactp1 = cg1(ik+1) / cg1(ik)
5883 ec_sig = esi * cm_sig(ispp1)/dmm(ik) * efactp1
5884 esum(isp) = esum(isp) - ec_sig*
va(ispp1,ip)
5890 dwni_m2(ik) = dble( cg1(ik) / dsip(ik) )
5892 #ifdef W3_DEBUGFREQSHIFT
5893 WRITE(740+
iaproc,*)
'Before FreqShift oper eSum=', sum(abs(esum))
5897 isp = ith + (ik -1)*nth
5898 ispm1 = ith + (ik-1 -1)*nth
5899 efactm1 = dble( cg1(ik-1) / cg1(ik) )
5900 ea_sig = - esi * dwni_m2(ik) * max(cwnb_m2(ispm1),zero) *efactm1
5901 esum(isp) = esum(isp) - ea_sig*
va(ispm1,ip)
5904 isp = ith + (ik -1)*nth
5905 ispp1 = ith + (ik+1 -1)*nth
5906 efactp1 = dble( cg1(ik+1) / cg1(ik) )
5907 ec_sig = esi * dwni_m2(ik) * min(cwnb_m2(isp),zero) * efactp1
5908 esum(isp) = esum(isp) - ec_sig*
va(ispp1,ip)
5911 #ifdef W3_DEBUGFREQSHIFT
5912 WRITE(740+
iaproc,*)
' after FreqShift oper eSum=', sum(abs(esum))
5917 WRITE(740+
iaproc,*)
'Step 4: sum(eSum)=', sum(esum)
5919 #ifdef W3_DEBUGSOLVERCOH
5920 pre_va(:, ip)=real(esum)
5922 esum(1:nspec) = esum(1:nspec) / aspar_diag(1:nspec)
5923 #ifdef W3_DEBUGFREQSHIFT
5924 WRITE(740+
iaproc,*)
'JSEA=', jsea,
' nbIter=', nbiter
5927 vanew(isp) = dble(cg(ik,isea)/clats(isea)) * esum(isp)
5930 vaanew(isp) = vanew(isp)
5931 vaaacloc(isp) = vaacloc(isp)
5934 vaanew(ith + nspec) = fachfa * vaanew(ith + nspec - nth)
5935 vaanew(ith - nth ) = 0.
5936 vaaacloc(ith + nspec) = fachfa * vaaacloc(ith + nspec - nth)
5937 vaaacloc(ith - nth ) = 0.
5940 vflwn(isp) = max(cwnb_m2(isp),0.) * vaanew(isp) + min(cwnb_m2(isp),0.) * vaanew(isp + nth)
5943 ediff(isp) = vanew(isp) -
vaold(isp) - dwni_m2(mapwn(isp)) * (vflwn(isp-nth) - vflwn(isp) )
5944 eval1=max(cwnb_m2(isp-nth),0.) * vaaacloc(isp-nth) + min(cwnb_m2(isp-nth),0.) * vaanew(isp)
5945 eval2=max(cwnb_m2(isp),0.) * vaanew(isp) + min(cwnb_m2(isp),0.) * vaaacloc(isp + nth)
5946 ediffb(isp) = vanew(isp) -
vaold(isp) - dwni_m2(mapwn(isp)) * (eval1 - eval2)
5948 IF (isea .eq. 190)
THEN
5951 isp = ith + (ik-1)*nth
5952 WRITE(740+
iaproc,*)
'ISP/ITH/IK=', isp, ith, ik
5953 WRITE(740+
iaproc,*)
'eDiff(A/B)=', ediff(isp), ediffb(isp)
5957 WRITE(740+
iaproc,*)
'NK=', nk,
' NTH=', nth
5961 isp = ith + (ik-1)*nth
5962 esumpart = esumpart + abs(ediff(isp))
5964 IF (isea .eq. 190)
THEN
5965 WRITE(740+
iaproc,*)
'IK=', ik,
' eSumDiff=', esumpart
5968 WRITE(740+
iaproc,*)
'sum(eDiff/VAnew/VAold)=', sum(abs(ediff)), sum(abs(vanew)), sum(abs(
vaold))
5971 IF (b_jgs_block_gauss_seidel)
THEN
5972 va(1:nspec,ip) = real(esum) * iobdp_loc(ip)
5976 isp = ith + (ik-1)*nth
5977 IF (
refpars(3) .LT. 0.5 .AND. iobpd_loc(ith,ip) .EQ. 0 .AND. iobpa_loc(ip) .EQ. 0)
THEN
5978 va(isp,ip) =
vaold(isp,ip) * iobdp_loc(ip)
5984 u_jac(1:nspec,ip) = esum
5987 esum =
va(1:nspec,ip)
5990 IF (b_jgs_terminate_difference)
THEN
5992 if (sum_new .gt. 0.d0)
then
5993 diffnew = abs(sum(acloc-esum))/sum_new
5994 #ifdef W3_DEBUGFREQSHIFT
5995 WRITE(740+
iaproc,*)
'DiffNew=', diffnew,
' Sum_new=', sum_new
5997 p_is_converged = diffnew
5999 p_is_converged = zero
6001 #ifdef W3_DEBUGFREQSHIFT
6002 WRITE(740+
iaproc,*)
'p_is_converged=', p_is_converged
6004 IF (p_is_converged .lt. b_jgs_diff_thr .and. nbiter .gt. 1)
then
6005 is_converged = is_converged + 1
6006 lconverged(ip) = .true.
6008 lconverged(ip) = .false.
6012 WRITE(740+
iaproc,*)
'sum(VA)out=', sum(
va(:,ip))
6016 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 2')
6018 #ifdef W3_DEBUGSOLVERCOH
6019 WRITE (efile,40) nbiter
6020 40
FORMAT (
'PRE_VA_',i4.4,
'.txt')
6026 IF (b_jgs_block_gauss_seidel)
THEN
6032 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 3')
6036 IF (b_jgs_terminate_maxiter)
THEN
6037 IF (nbiter .gt. b_jgs_maxiter)
THEN
6038 #ifdef W3_DEBUGSOLVER
6039 WRITE(740+
iaproc,*)
'Exiting by TERMINATE_MAXITER'
6044 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 4')
6048 IF (b_jgs_terminate_difference .and. int(mod(nbiter,10)) == 0)
THEN
6049 CALL mpi_allreduce(is_converged, itmp, 1, mpi_int, mpi_sum, mpi_comm_wcmp, ierr)
6051 prop_conv = (dble(nx) - dble(is_converged))/dble(nx) * 100.
6052 #ifdef W3_DEBUGSOLVER
6053 WRITE(740+
iaproc,*)
'solver', nbiter, is_converged, prop_conv, b_jgs_pmin
6056 IF (myrank == 0)
WRITE(*,*)
'No. of solver iterations', nbiter, is_converged, prop_conv, b_jgs_pmin
6057 IF (prop_conv .le. b_jgs_pmin + tiny(1.))
THEN
6058 #ifdef W3_DEBUGFREQSHIFT
6059 WRITE(740+
iaproc,*)
'prop_conv=', prop_conv
6060 WRITE(740+
iaproc,*)
'NX=', nx
6061 WRITE(740+
iaproc,*)
'is_converged=', is_converged
6062 WRITE(740+
iaproc,*)
'Exiting by TERMINATE_DIFFERENCE'
6067 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 5')
6071 IF (b_jgs_terminate_norm)
THEN
6075 IF (iobp_loc(ip).eq.1)
THEN
6080 isea= mapfs(1,ip_glob)
6081 esum(:) = esum(:) - aspar_diag(:)*acloc
6082 DO i = pdlib_ia_p(ip)+1, pdlib_ia_p(ip+1)
6086 IF (fsrefraction)
THEN
6089 ispprevdir=listispprevdir(isp)
6090 ispnextdir=listispnextdir(isp)
6091 ea_the = - dtg*esi*max(zero,cad(ispprevdir))
6092 ec_the = dtg*esi*min(zero,cad(ispnextdir))
6093 esum(isp) = esum(isp) - ea_the*
va(ispprevdir,ip)
6094 esum(isp) = esum(isp) - ec_the*
va(ispnextdir,ip)
6097 IF (fsfreqshift)
THEN
6099 cp_sig = max(zero,cas)
6100 cm_sig = min(zero,cas)
6103 CALL wavnu_local(sig(ik),dw(isea),wn1(ik),cg1(ik))
6105 cg1(ik) = cg(ik,isea)
6106 wn1(ik) = wn(ik,isea)
6110 IF (iobpd_loc(ith,ip) .NE. 0)
THEN
6112 isp =ith + (ik -1)*nth
6113 ispm1=ith + (ik-1-1)*nth
6114 efactm1=cg(ik-1,isea) / cg1(ik)
6115 ea_sig= - esi*cp_sig(ispm1)/dmm(ik-1) * efactm1
6116 esum(isp) = esum(isp) - ea_sig*
va(ispm1,ip)
6119 isp =ith + (ik -1)*nth
6120 ispp1=ith + (ik+1-1)*nth
6121 efactp1=cg(ik+1,isea) / cg1(ik)
6122 ec_sig= esi*cm_sig(ispp1)/dmm(ik) * efactp1
6123 esum(isp) = esum(isp) - ec_sig*
va(ispp1,ip)
6128 sum_l2 = sum_l2 + sum(esum*esum)
6131 CALL mpi_allreduce(sum_l2, sum_l2_gl, 1,
rtype, mpi_sum, mpi_comm_wcmp, ierr)
6132 #ifdef W3_DEBUGSOLVER
6133 WRITE(740+
iaproc,*)
'Sum_L2_gl=', sum_l2_gl
6136 IF (sum_l2_gl .le. b_jgs_norm_thr)
THEN
6137 #ifdef W3_DEBUGFREQSHIFT
6138 WRITE(740+
iaproc,*)
'Exiting by TERMINATE_NORM'
6143 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION SOLVER LOOP 6')
6149 #ifdef W3_DEBUGSOLVER
6150 WRITE(740+
iaproc,*)
'nbIter=', nbiter,
' B_JGS_MAXITER=', b_jgs_maxiter
6156 WRITE(740+
iaproc,*)
'IOBPD loop, Before, sum(VA)=', sum(
va(:,ip))
6159 ith = 1 + mod(isp-1,nth)
6160 va(isp,ip)=max(zero,
va(isp,ip))*iobdp_loc(ip)*dble(iobpd_loc(ith,ip))
6162 IF (
refpars(3).LT.0.5.AND.iobpd_loc(ith,ip).EQ.0.AND.iobpa_loc(ip).EQ.0)
THEN
6168 WRITE(740+
iaproc,*)
'IOBPD loop, After, sum(VA)=', sum(
va(:,ip))
6171 #ifdef W3_DEBUGSOLVERCOH
6174 #ifdef W3_DEBUGSOLVER
6183 isea = mapfs(1,ip_glob)
6197 ik = 1 + (isp-1)/nth
6199 CALL wavnu_local(sig(ik),dw(isea),wn1(ik),cg1(ik))
6201 cg1(ik) = cg(ik,isea)
6203 eva = max( zero ,cg1(ik)/clats(isea)*real(
va(isp,ip)) )
6204 evo = max( zero ,cg1(ik)/clats(isea)*real(
vaold(isp,jsea)) )
6206 sumacout=sumacout + real(
va(isp,ip))
6207 vs_w3srce =
vstot(isp,jsea) * dtg / max(1., (1. - dtg*
vdtot(isp,jsea)))
6208 eva_w3srce = max(0.,
va(isp,jsea) + vs_w3srce)
6209 intdiff = intdiff + abs(eva - eva_w3srce)
6211 eb=
va(isp,jsea) + dtg*(
vstot(isp,jsea) -
vdtot(isp,jsea)*
va(isp,jsea))
6212 evasolve=max(0., cg(ik,isea)/clats(isea)*acsolve)
6213 vasolve(isp)=evasolve
6214 sumvs = sumvs + abs(
vstot(isp,jsea))
6215 sumvd = sumvd + abs(
vdtot(isp,jsea))
6216 sumvain = sumvain + abs(
va(isp,jsea))
6217 sumvaout = sumvaout + abs(eva)
6218 sumvaw3srce = sumvaw3srce + abs(eva_w3srce)
6220 vaold(isp,jsea) = evo
6224 WRITE(740+
iaproc,*)
'ISEA=', isea,
' IntDiff=', intdiff,
' DTG=', dtg
6225 IF (isea .eq. testnode)
THEN
6227 WRITE(740+
iaproc,*)
'ISP=', isp,
'VA/VAsolve=',
va(isp,jsea), vasolve(isp)
6230 WRITE(740+
iaproc,*)
'SHAVE=', shavetot(jsea)
6231 WRITE(740+
iaproc,*)
'Sum(VS/VD)=', sumvs, sumvd
6232 WRITE(740+
iaproc,*)
'min/max/sum(VS)=', minval(
vstot(:,jsea)), maxval(
vstot(:,jsea)), sum(
vstot(:,jsea))
6233 WRITE(740+
iaproc,*)
'min/max/sum(VD)=', minval(
vdtot(:,jsea)), maxval(
vdtot(:,jsea)), sum(
vdtot(:,jsea))
6234 WRITE(740+
iaproc,*)
'min/max/sum(VA)=', minval(
va(:,jsea)), maxval(
va(:,jsea)), sum(
va(:,jsea))
6235 WRITE(740+
iaproc,*)
'min/max/sum(VAsolve)=', minval(vasolve), maxval(vasolve), sum(vasolve)
6236 WRITE(740+
iaproc,*)
'SumVA(in/out/w3srce)=', sumvain, sumvaout, sumvaw3srce
6237 WRITE(740+
iaproc,*)
'SumACout=', sumacout
6241 IF (b_jgs_limiter)
THEN
6244 ik = 1 + (isp-1)/nth
6245 spec(isp) =
vaold(isp,jsea)
6248 CALL w3spr4 (spec, cg1, wn1, emean, fmean, fmean1, wnmean, &
6249 amax,
u10(isea),
u10d(isea), &
6251 taua, tauadir, dair, &
6254 tauwx, tauwy, cd, z0, charn, llws, fmeanws, dlwmean)
6259 dam(1+(ik-1)*nth) = 0.0081*0.1 / ( 2 * sig(ik) * wn(ik,isea)**3 * cg(ik,isea)) * cg1(ik) / clats(isea)
6265 dam(ith+is0) = dam(1+is0)
6271 jac2 = 1./
tpi/sig(ik)
6273 dam2(1+(ik-1)*nth) = 1e-06 *
grav/frlocal**4 * ustar * max(fmeanws,fmean) * dtg * jac2 * cg1(ik) / clats(isea)
6278 dam2(ith+is0) = dam2(1+is0)
6284 isp = ith + (ik-1)*nth
6285 newdac =
va(isp,ip) -
vaold(isp,jsea)
6286 maxdac = max(dam(isp),dam2(isp))
6287 newdac = sign(min(maxdac,abs(newdac)), newdac)
6288 va(isp,ip) = max(0.,
vaold(isp,ip) + newdac)
6296 INQUIRE (
file=
'weights.ww3', exist = lexist )
6297 if (.not. lexist)
then
6301 ipitergl(iplg(ip)) = ipiter(ip)
6303 call mpi_reduce(ipitergl,ipiterout,np_global,mpi_int,mpi_sum,0,mpi_comm_wcmp,ierr)
6304 if (myrank == 0)
tHEN
6305 OPEN(100001,
file=
'weights.ww3',form=
'FORMATTED',status=
'unknown')
6306 do ip = 1, np_global
6307 write(100001,*) ipiterout(ip)
6314 call print_memcheck(
memunit,
'memcheck_____:'//
' WW3_PROP SECTION LOOP 7')
6318 WRITE(740+
iaproc,*)
'JSEA=', jsea
6319 WRITE(740+
iaproc,*)
'min/max/sum(VA)=', minval(
va(:,jsea)), maxval(
va(:,jsea)), sum(
va(:,jsea))
6321 WRITE(740+
iaproc,*)
'min/max/sum(VAtot)=', minval(
va), maxval(
va), sum(
va)
6325 #ifdef W3_DEBUGSOLVER
6326 WRITE(740+
iaproc,*)
'PDLIB_JACOBI_GAUSS_SEIDEL_BLOCK, end'
6379 USE w3gdatmd,
ONLY:
nk,
nth,
nspec,
sig,
dth,
esin,
ecos,
nseal,
fsbccfl,
clats,
mapfs
6393 USE mpi,
only : mpi_min
6401 LOGICAL,
INTENT(IN) :: LCALC
6403 INTEGER,
INTENT(IN) :: IMOD
6405 REAL,
INTENT(IN) :: FACX, FACY, DTG, VGX, VGY
6407 REAL :: KTMP(3), UTILDE(NTH), ST(NTH,NPA)
6408 REAL :: FL11(NTH), FL12(NTH), FL21(NTH), FL22(NTH), FL31(NTH), FL32(NTH), KKSUM(NTH,NPA)
6409 REAL :: FL111(NTH), FL112(NTH), FL211(NTH), FL212(NTH), FL311(NTH), FL312(NTH)
6411 REAL :: KSIG(NPA), CGSIG(NPA), CXX(NTH,NPA), CYY(NTH,NPA)
6412 REAL :: LAMBDAX(NTH), LAMBDAY(NTH)
6413 REAL :: DTMAX(NTH), DTMAXEXP(NTH), DTMAXOUT, DTMAXGL
6414 REAL :: FIN(1), FOUT(1), REST, CFLXY, RD1, RD2, RD10, RD20
6415 REAL :: UOLD(NTH,NPA), U(NTH,NPA)
6417 REAL,
PARAMETER :: ONESIXTH = 1.0/6.0
6418 REAL,
PARAMETER :: ZERO = 0.0
6419 REAL,
PARAMETER :: THR = 1.0e-12
6421 INTEGER :: IK, ISP, ITH, IE, IP, IT, IBI, NI(3), I1, I2, I3, JX, IERR, IP_GLOB, ISEA
6435 CALL wavnu3 (sig(ik),
dw(
iplg(ip)), ksig(ip), cgsig(ip))
6441 cxx(ith,ip) = cgsig(ip) * facx * ecos(ith) / clats(isea)
6442 cyy(ith,ip) = cgsig(ip) * facy * esin(ith)
6447 IF (iobp_loc(ip) .GT. 0)
THEN
6448 cxx(ith,ip) = cxx(ith,ip) + facx *
cx(isea)/clats(isea)
6449 cyy(ith,ip) = cyy(ith,ip) + facy *
cy(isea)
6464 lambdax(ith) = onesixth *(cxx(ith,i1)+cxx(ith,i2)+cxx(ith,i3))
6465 lambday(ith) = onesixth *(cyy(ith,i1)+cyy(ith,i2)+cyy(ith,i3))
6469 ktmp(1) =
kelem1(ith,ie,ik)
6470 ktmp(2) =
kelem2(ith,ie,ik)
6471 ktmp(3) =
kelem3(ith,ie,ik)
6472 nm(ith,ie,ik) = - 1.d0/min(-thr,sum(min(zero,ktmp)))
6473 kelem1(ith,ie,ik) = max(zero,ktmp(1))
6474 kelem2(ith,ie,ik) = max(zero,ktmp(2))
6475 kelem3(ith,ie,ik) = max(zero,ktmp(3))
6485 fl111 = 2.d0 * fl11 + fl12
6486 fl112 = 2.d0 * fl12 + fl11
6487 fl211 = 2.d0 * fl21 + fl22
6488 fl212 = 2.d0 * fl22 + fl21
6489 fl311 = 2.d0 * fl31 + fl32
6490 fl312 = 2.d0 * fl32 + fl31
6492 flall1(:,ie,ik) = (fl311 + fl212) * onesixth +
kelem1(:,ie,ik)
6493 flall2(:,ie,ik) = (fl111 + fl312) * onesixth +
kelem2(:,ie,ik)
6494 flall3(:,ie,ik) = (fl211 + fl112) * onesixth +
kelem3(:,ie,ik)
6502 kksum(ith,ni(1)) = kksum(ith,ni(1)) +
kelem1(ith,ie,ik)
6503 kksum(ith,ni(2)) = kksum(ith,ni(2)) +
kelem2(ith,ie,ik)
6504 kksum(ith,ni(3)) = kksum(ith,ni(3)) +
kelem3(ith,ie,ik)
6511 IF (iobp_loc(ip) .EQ. 1 .OR. fsbccfl)
THEN
6513 dtmaxexp(ith) =
pdlib_si(ip)/max(thr,kksum(ith,ip)*iobdp_loc(ip))
6514 dtmax(ith) = min(dtmax(ith),dtmaxexp(ith))
6516 dtmaxout = minval(dtmax)
6524 cflxy = dble(dtg)/dtmaxgl
6525 rest = abs(mod(cflxy,1.0d0))
6526 IF (rest .LT. thr)
THEN
6527 iter(ik) = abs(nint(cflxy))
6528 ELSE IF (rest .GT. thr .AND. rest .LT. 0.5d0)
THEN
6529 iter(ik) = abs(nint(cflxy)) + 1
6531 iter(ik) = abs(nint(cflxy))
6544 isp = ith + (ik-1)*nth
6545 u(ith,ip) =
va(isp,ip) / cgsig(ip) * clats(
iplg(ip))
6555 utilde(ith) =
nm(ith,ie,ik) * (
flall1(ith,ie,ik)*u(ith,ni(1)) +
flall2(ith,ie,ik)*u(ith,ni(2)) +
flall3(ith,ie,ik)*u(ith,ni(3)))
6556 st(ith,ni(1)) = st(ith,ni(1)) +
kelem1(ith,ie,ik) * (u(ith,ni(1)) - utilde(ith))
6557 st(ith,ni(2)) = st(ith,ni(2)) +
kelem2(ith,ie,ik) * (u(ith,ni(2)) - utilde(ith))
6558 st(ith,ni(3)) = st(ith,ni(3)) +
kelem3(ith,ie,ik) * (u(ith,ni(3)) - utilde(ith))
6563 u(ith,ip) = max(zero,u(ith,ip)-
dtsi(ip)*st(ith,ip)*(1-iobpa_loc(ip)))*iobpd_loc(ith,ip)*iobdp_loc(ip)
6565 IF (
refpars(3).LT.0.5.AND.iobpd_loc(ith,ip).EQ.0.AND.iobpa_loc(ip).EQ.0) u(ith,ip) = uold(ith,ip)
6572 isp = ith + (ik-1) * nth
6573 rd1 = rd10 - dtg * real(
iter(ik)-it)/real(
iter(ik))
6575 IF ( rd2 .GT. 0.001 )
THEN
6576 rd2 = min(1.,max(0.,rd1/rd2))
6583 ip_glob = mapsf(isbpi(ibi),1)
6586 u(ith,jx) = ( rd1*bbpi0(isp,ibi) + rd2*bbpin(isp,ibi) ) / cgsig(isbpi(ibi)) * clats(isbpi(ibi))
6600 isp = ith + (ik-1)*nth
6601 va(isp,ip) = u(ith,ip) * cgsig(ip) / clats(
iplg(ip))
6737 INTEGER,
INTENT(IN) :: IMOD
6741 INTEGER ISP, ITH, IK, ISPprevFreq, ISPnextFreq
6742 INTEGER NewISP, JTH, istat
6752 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 1'
6758 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 2'
6762 ith = 1 + mod(isp-1,nth)
6763 ik = 1 + (isp-1)/nth
6767 ispprevfreq=ith + (ik-1 -1)*nth
6770 IF (ik .eq. nk)
THEN
6773 ispnextfreq=ith + (ik+1 -1)*nth
6777 IF (ith .eq. 1)
THEN
6782 newisp=jth + (ik-1)*nth
6783 listispprevdir(isp)=newisp
6784 IF (ith .eq. nth)
THEN
6789 newisp=jth + (ik-1)*nth
6790 listispnextdir(isp)=newisp
6793 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 3'
6798 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 4'
6803 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 5'
6808 WRITE(740+
iaproc,*)
'BLOCK_SOLVER_INIT, step 6'
6876 INTEGER,
SAVE :: IENT = 0
6882 INTEGER :: JSEA, ISEA, IX, IP, IP_glob
6883 real*8,
PARAMETER :: dthr = 10e-6
6885 CALL strace (ient,
'SETDEPTH_PDLIB')
6890 IF (
dw(ip_glob) .LT. dmin + dthr)
THEN
6965 INTEGER,
SAVE :: IENT = 0
6971 INTEGER :: JSEA, ISEA, IX, IP, IP_glob
6972 real*8,
PARAMETER :: dthr = 10e-6
6974 CALL strace (ient,
'SETDEPTH_PDLIB')
6977 ip_glob =
iplg(jsea)
6978 IF (mapsta(1,ip_glob).EQ.2)
THEN
7059 refpars, reflc, refld, &
7061 angle0, angle, nseal
7081 INTEGER :: ITH, IX, I, J, IP, IE, NDIRSUM
7082 REAL (KIND = 8) :: cossum, sinsum
7083 REAL (KIND = 8) :: dirmin, dirmax, shift, tempo, dircoast
7084 REAL (KIND = 8) :: x1, x2, y1, y2, dxp1, dxp2, dxp3
7085 REAL (KIND = 8) :: dyp1, dyp2, dyp3, edet1, edet2, evx, evy
7086 REAL(KIND=8), parameter :: thr = tiny(1.)
7087 INTEGER :: I1, I2, I3
7088 INTEGER :: ITMP(NX), NEXTVERT(NX), PREVVERT(NX)
7089 INTEGER :: MAX_IOBPD, MIN_IOBPD
7091 CHARACTER(60) :: FNAME
7093 INTEGER,
SAVE :: IENT = 0
7118 ELSE IF (i.eq.2)
THEN
7124 ELSE IF (i.eq.3)
THEN
7132 edet1 = thr-x1*evy+y1*evx
7133 edet2 = thr+x2*evy-y2*evx
7134 IF ((edet1.gt.0.).and.(edet2.gt.0.))
THEN
7153 IF (max_iobpd .gt. 1 .OR. min_iobpd .lt. 0)
THEN
7154 WRITE(*,*)
'MAX_IOBPD - MIN_IOBPD', max_iobpd, min_iobpd
7155 stop
'MAX_IOBPD ERRROR'
7158 #ifdef W3_DEBUGSETUGIOBP
7159 WRITE(740+
iaproc,*)
'Calling SETUGIOBP, step 5'
7165 #ifdef W3_DEBUGSETUGIOBP
7166 WRITE(740+
iaproc,*)
'Calling SETUGIOBP, step 7'
7179 IF (iobp(ip).EQ.0.AND.mapsta(1,ip).EQ.1)
THEN
7184 cossum=cossum+iobpd(ith,ip)*ecos(ith)
7185 sinsum=sinsum+iobpd(ith,ip)*esin(ith)
7186 ndirsum=ndirsum+iobpd(ith,ip)
7188 dircoast=atan2(sinsum, cossum)
7189 refld(1,mapfs(1,ip)) = 1+mod(nth+nint(dircoast/dth),nth)
7190 refld(2,mapfs(1,ip)) = 4-max(2,nint(4.*real(ndirsum)/real(nth)))
7191 reflc(1,mapfs(1,ip))= refpars(1)
7195 #ifdef W3_DEBUGSETUGIOBP
7196 WRITE(740+
iaproc,*)
'Calling SETUGIOBP, step 8'
7260 INTEGER,
SAVE :: IENT = 0
7267 CALL strace (ient,
'BLOCK_SOLVER_FINALIZE')
7342 INTEGER,
SAVE :: IENT = 0
7347 INTEGER,
INTENT(IN) :: IMOD
7353 grids(imod)%CROSSDIFF, &
7356 grids(imod)%ANGLE, &
7357 grids(imod)%ANGLE0, &
7359 grids(imod)%COUNTCON, &
7360 grids(imod)%INDEX_CELL, &
7361 grids(imod)%IE_CELL, &
7362 grids(imod)%POS_CELL, &
7366 grids(imod)%I_DIAG, &
7367 grids(imod)%JA_IE, &
7370 grids(imod)%IOBDP, &
7377 SUBROUTINE ergout(FHNDL, ERGNAME)
7428 INTEGER,
INTENT(IN) :: FHNDL
7429 CHARACTER(LEN=*),
INTENT(IN) :: ERGNAME
7430 REAL :: SUMVA(NSEAL)
7434 OPEN(fhndl,
file = trim(ergname), form =
'UNFORMATTED')
7441 sumva(jsea) = sum(
va(:,jsea))
7445 WRITE(fhndl) (sumva(jsea), sumva(jsea), sumva(jsea), jsea = 1, nseal)
7514 INTEGER,
SAVE :: IENT = 0
7519 INTEGER,
INTENT(IN) :: IMOD
7525 ELSE IF (
imem == 2)
THEN
7540 IF (.NOT. b_jgs_block_gauss_seidel)
THEN
7541 ALLOCATE(
u_jac(nspec,
npa), stat=istat)
7605 INTEGER,
SAVE :: IENT = 0
7611 CALL strace (ient,
'JACOBI_FINALIZE')
7615 ELSE IF (
imem == 2)
THEN