diff --git a/mitgcm_code/forward_step.F b/mitgcm_code/forward_step.F index 8a7a484..95d3665 100644 --- a/mitgcm_code/forward_step.F +++ b/mitgcm_code/forward_step.F @@ -1145,6 +1145,8 @@ SUBROUTINE FORWARD_STEP( iloop, myTime, myIter, myThid ) C-- Acoustical physics CALL IHOP_SOUND_SPEED( myThid ) CALL IHOP_DRIVER( myTime,myIter,myThid ) +CML I think that this could go here instead of being called twice in the_main_loop +CML CALL IHOP_COST_INLOOP( myTime, myThid ) #endif /* ALLOW_IHOP */ #ifdef ALLOW_DIAGNOSTICS diff --git a/mitgcm_code/the_main_loop.F b/mitgcm_code/the_main_loop.F index d5fe08e..7263bc8 100644 --- a/mitgcm_code/the_main_loop.F +++ b/mitgcm_code/the_main_loop.F @@ -73,6 +73,16 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) C \ev C !USES: +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_IHOP + use arr_mod, only: arr, narr + use bdry_mod, only: bdry + use ssp_mod, only: ssp + use srpos_mod, only: pos + use angle_mod, only: angles + use ihop_mod, only: beam, ray2d +# endif +#endif IMPLICIT NONE C == Global variables == #include "SIZE.h" @@ -428,6 +438,11 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) # endif max_lev2=nTimeSteps/nchklev_1+1 +C Declaring variables as passive does not work properly +CML# ifdef ALLOW_IHOP +CML!$TAF PASSIVE angles +CML# endif + c************************************** # ifdef ALLOW_DIVIDED_ADJOINT CADJ loop = divided @@ -443,6 +458,10 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) # endif CALL AUTODIFF_STORE( myThid ) #include "checkpoint_lev4_directives.h" +# ifdef ALLOW_IHOP +!$TAF STORE arr, bdry, narr, pos, ssp = tapelev4, key = ilev_4 +!$TAF STORE angles, beam, ray2d = tapelev4, key = ilev_4 +# endif CALL AUTODIFF_RESTORE( myThid ) # ifdef ALLOW_AUTODIFF_WHTAPEIO CALL AUTODIFF_WHTAPEIO_SYNC( 4 , 1, myThid ) @@ -462,6 +481,10 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) # endif CALL AUTODIFF_STORE( myThid ) #include "checkpoint_lev3_directives.h" +# ifdef ALLOW_IHOP +!$TAF STORE arr, bdry, narr, pos, ssp = tapelev3, key = ilev_3 +!$TAF STORE angles, beam, ray2d = tapelev3, key = ilev_3 +# endif CALL AUTODIFF_RESTORE( myThid ) # ifdef ALLOW_AUTODIFF_WHTAPEIO CALL AUTODIFF_WHTAPEIO_SYNC( 3 , 1, myThid ) @@ -480,6 +503,19 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) # endif CALL AUTODIFF_STORE( myThid ) #include "checkpoint_lev2_directives.h" +# ifdef ALLOW_IHOP +CML!$TAF PASSIVE angles +CML!$TAF STORE arr, bdry, narr, pos, ssp = tapelev2, key = ilev_2 +CML!$TAF STORE angles, beam, ray2d = tapelev2, key = ilev_2 +!$TAF STORE arr, narr, ray2d = tapelev2, key = ilev_2 +!$TAF STORE angles%dalpha = tapelev2, key = ilev_2 +!$TAF STORE bdry%bot, bdry%top = tapelev2, key = ilev_2 +!$TAF STORE beam%nsteps = tapelev2, key = ilev_2 +!$TAF STORE pos%rz, pos%sz = tapelev2, key = ilev_2 +!$TAF STORE ssp%c, ssp%cmat = tapelev2, key = ilev_2 +!$TAF STORE ssp%cz, ssp%czmat = tapelev2, key = ilev_2 +!$TAF STORE ssp%rho = tapelev2, key = ilev_2 +# endif CALL AUTODIFF_RESTORE( myThid ) # ifdef ALLOW_AUTODIFF_WHTAPEIO CALL AUTODIFF_WHTAPEIO_SYNC( 2 , 1, myThid ) @@ -520,7 +556,11 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) c-- # ifdef ALLOW_IHOP ! for ssp_mod>gcmSSP; make local? -!$TAF INIT comlev1_bibj_ij_ihop = COMMON, nchklev_1*nSx*nSy*(2*OLy+sNy)*(2*OLx+sNx) +!$TAF INIT comlev1_ihop = STATIC, nchklev_1 +!$TAF INIT comlev1_ihop_nts = STATIC, nchklev_1*nts +C global tapes for gcmSSP +!$TAF INIT comlev1_ihop_iijj = STATIC, nchklev_1*nSx*nSy*sNx*sNy*IHOP_MAX_RANGE*IHOP_MAX_NC_SIZE +!$TAF INIT comlev1_ihop_iijj_k = STATIC, nchklev_1*nSx*nSy*sNx*sNy*IHOP_MAX_RANGE*IHOP_MAX_NC_SIZE*(Nr + 2) # endif /* ALLOW_IHOP */ c-- # ifdef ALLOW_MOM_COMMON @@ -698,16 +738,30 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) ENDIF #endif #ifdef ALLOW_IHOP +# ifdef ALLOW_AUTODIFF_TAMC +CML!$TAF PASSIVE angles +CML!$TAF STORE arr,narr,bdry,pos,ssp = comlev1_ihop, key = ikey_dynamics +!$TAF STORE arr, narr, ray2d = comlev1_ihop, key = ikey_dynamics +!$TAF STORE angles%dalpha = comlev1_ihop, key = ikey_dynamics +!$TAF STORE bdry%bot, bdry%top = comlev1_ihop, key = ikey_dynamics +!$TAF STORE beam%nsteps = comlev1_ihop, key = ikey_dynamics +!$TAF STORE pos%rz, pos%sz = comlev1_ihop, key = ikey_dynamics +CML!$TAF STORE ssp%c, ssp%cmat = comlev1_ihop, key = ikey_dynamics +CML!$TAF STORE ssp%cz, ssp%czmat = comlev1_ihop, key = ikey_dynamics +C for some reason TAF does not do this right, so we store ssp instead +CML!$TAF STORE ssp%rho = comlev1_ihop, key = ikey_dynamics +!$TAF STORE ssp = comlev1_ihop, key = ikey_dynamics +#endif IF (useIHOP) THEN #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_CALL('ihop_cost_inloop',myThid) #endif c-- Accumulate in-situ acoutsic travel times. -#ifdef ALLOW_AUTODIFF +# ifdef ALLOW_AUTODIFF C-- Reset the model iteration counter and the model time. myIter = nIter0 + (iloop-1) myTime = startTime + deltaTClock*(iloop-1) -#endif +# endif CALL TIMER_START('IHOP_COST_INLOOP [MAIN_DO_LOOP]', myThid) CALL IHOP_COST_INLOOP( myTime, myThid ) CALL TIMER_STOP ('IHOP_COST_INLOOP [MAIN_DO_LOOP]', myThid) @@ -771,6 +825,10 @@ SUBROUTINE THE_MAIN_LOOP( myTime, myIter, myThid ) #endif #ifdef ALLOW_IHOP IF (useIHOP) THEN +# ifdef ALLOW_AUTODIFF +!$TAF STORE arr, narr = comlev1_ihop, key = ikey_dynamics +!$TAF STORE pos%rz, pos%sz = comlev1_ihop, key = ikey_dynamics +# endif c-- Accumulate travel times CALL TIMER_START('IHOP_COST_INLOOP [THE_MAIN_LOOP]', myThid) CALL IHOP_COST_INLOOP( endtime, myThid ) diff --git a/src/ihop.F90 b/src/ihop.F90 index 1a2a597..462aed2 100644 --- a/src/ihop.F90 +++ b/src/ihop.F90 @@ -1,4 +1,7 @@ #include "IHOP_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif !BOP ! !ROUTINE: IHOP ! !INTERFACE: @@ -59,10 +62,16 @@ SUBROUTINE IHOP_MAIN ( myTime, myIter, myThid ) USE ihop_init_diag, only: initPRTFile, openOutputFiles, resetMemory USE bdry_mod, only: Bdry, writeBdry USE ssp_mod, only: setSSP - USE refCoef, only: writeRefCoef +#ifdef ALLOW_AUTODIFF_TAMC + USE ssp_mod, only: SSP +#endif + USE refCoef, only: writeRefCoef USE beampat, only: writePat USE ihop_mod, only: Beam +#ifdef ALLOW_AUTODIFF +# include "tamc.h" +#endif ! == Routine Arguments == ! myThid :: Thread number. Unused by IESCO ! msgBuf :: Used to build messages for printing. @@ -91,16 +100,19 @@ SUBROUTINE IHOP_MAIN ( myTime, myIter, myThid ) CALL writeBdry ( myThid ) ! write refCoef: OPTIONAL - CALL writeRefCoef( myThid ) + CALL writeRefCoef( myThid ) ! Source Beam Pattern: OPTIONAL, default is omni source pattern CALL writePat( myThid ) ! set SSP%cmat from gcm SSP: REQUIRED +#ifdef ALLOW_AUTODIFF_TAMC +!$TAF STORE ssp%c, ssp%cz = comlev1_ihop, key = ikey_dynamics +!$TAF STORE ssp%czmat, ssp%rho = comlev1_ihop, key = ikey_dynamics +#endif CALL setSSP( myThid ) - ! open all output files IF ( IHOP_dumpfreq .GE. 0 ) & CALL OpenOutputFiles( IHOP_fileroot, myTime, myIter, myThid ) diff --git a/src/ihop_cost_inloop.F b/src/ihop_cost_inloop.F index 3de521a..5412ed3 100644 --- a/src/ihop_cost_inloop.F +++ b/src/ihop_cost_inloop.F @@ -86,10 +86,6 @@ SUBROUTINE IHOP_COST_INLOOP( myTime, myThid ) & .AND.(ihopObs_time(num_file,m,bi,bj).LT. & (myTime+deltaTclock))) THEN -#ifdef ALLOW_AUTODIFF_TAMC -!$TAF STORE ihop_modval = comlev1_bibj, key=itdkey, byte=isbyte -#endif - ihop_modval=0 CALL ihop_cost_modval(ihop_modval, num_file, diff --git a/src/ihop_driver.F b/src/ihop_driver.F index 5e4a54d..2d5c16f 100644 --- a/src/ihop_driver.F +++ b/src/ihop_driver.F @@ -36,16 +36,18 @@ SUBROUTINE IHOP_DRIVER( myTime, myIter, myThid ) C !LOCAL VARIABLES: ==================================================== INTEGER t + LOGICAL do_ihop_main + CEOP #ifdef ALLOW_IHOP + do_ihop_main=.FALSE. DO t=1,nts - IF ( IHOP_iter(t).GE.0 ) THEN - IF ( IHOP_iter(t).EQ.myIter ) THEN - CALL IHOP_MAIN( myTime, myIter, myThid ) - ENDIF + IF ( IHOP_iter(t).GE.0 .AND. IHOP_iter(t).EQ.myIter ) THEN + do_ihop_main=.TRUE. ENDIF ENDDO + IF ( do_ihop_main ) CALL IHOP_MAIN( myTime, myIter, myThid ) #endif /* ALLOW_IHOP */ RETURN diff --git a/src/splinec_mod.F90 b/src/splinec_mod.F90 index 87d35f5..760562e 100644 --- a/src/splinec_mod.F90 +++ b/src/splinec_mod.F90 @@ -79,7 +79,9 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) L = N - 1 -!$TAF init cspline1 = static, 4*ndim*L +!$TAF init cspline1 = static, 1 +!$TAF init cspline1_ndim = static, NDIM +!ML!$TAF init cspline1 = static, 4*NDIM*L ! not sure if this correct DO M = 2,N C(3,M) = TAU(M) - TAU(M-1) @@ -88,7 +90,7 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) ! * BEGINNING BOUNDARY CONDITION SECTION * -!$TAF store C = cspline1 +!$TAF INCOMPLETE C IF (IBCBEG==0) THEN ! IBCBEG = 0 !$TAF store C = cspline1 IF (N.GT.2) THEN ! N > 2 @@ -121,17 +123,20 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) ! * RUNNING CALCULATIONS TO N-1 - LOOP IS NOT EXECUTED IF N = 2 * DO M = 2,L -!$TAF store C = cspline1 +!$TAF store C = cspline1_ndim, key = M G = -C(3,M+1) / C(4,M-1) C(2,M) = G*C(2,M-1) + 3.0*(C(3,M)*C(4,M+1) + C(3,M+1)*C(4,M)) C(4,M) = G*C(3,M-1) + 2.0*(C(3,M) + C(3,M+1)) END DO +!$TAF store G = cspline1 ! * ENDING BOUNDARY CONDITION SECTION * +!$TAF INCOMPLETE C IF (IBCEND /= 1) THEN !$TAF store C = cspline1 IF (IBCEND==0) THEN +!$TAF INCOMPLETE C IF (N==2 .AND. IBCBEG==0) THEN C(2,N) = C(4,N) ELSE IF ((N==3 .AND. IBCBEG==0) .OR. N==2) THEN @@ -154,10 +159,10 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) G = -1.0 / C(4,N-1) ELSE ! do nothing - C(2,N) = C(2,N) - C(4,N) = C(4,N) + C(2,N) = C(2,N) + C(4,N) = C(4,N) G = G - END IF + ENDIF IF ( IBCBEG.GT.0 .OR. N.GT.2) THEN !$TAF store C = cspline1 @@ -165,6 +170,7 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) C(2,N) = (G*C(2,N-1) + C(2,N)) / C(4,N) END IF END IF +!$TAF store C = cspline1 ! * RUN THE ENDING BOUNDARY EFFECT BACK THROUGH THE EQUATIONS * @@ -194,7 +200,7 @@ SUBROUTINE CSPLINE (TAU, C, N, IBCBEG, IBCEND, NDIM) C(4,N) = (0.0,0.0) DO I = 1,L ! INTEGRATE OVER THE INTERVAL -!$TAF store C = cspline1 +!$TAF store C(4,:) = cspline1_ndim, key = I DTAU = TAU(I+1) - TAU(I) C(4,N) = C(4,N) + DTAU*(C(1,I) + DTAU*(C(2,I)/2.0 + & DTAU*(C(3,I)/6.0 + DTAU*C(4,I)/24.0))) diff --git a/src/ssp_mod.F90 b/src/ssp_mod.F90 index 65de9b6..9455438 100644 --- a/src/ssp_mod.F90 +++ b/src/ssp_mod.F90 @@ -120,6 +120,9 @@ SUBROUTINE setSSP( myThid ) ! Set SSP structures for particular SSP%Type routine +#ifdef ALLOW_AUTODIFF +# include "tamc.h" +#endif ! == Routine Arguments == ! myThid :: Thread number. Unused by IESCO ! msgBuf :: Used to build messages for printing. @@ -129,14 +132,16 @@ SUBROUTINE setSSP( myThid ) ! == Local Variables == INTEGER :: ir, iz -!$TAF init setssp1 = 'ssp_mod_setssp' +!ML!$TAF init setssp1 = 'ssp_mod_setssp' ! IESCO24: Write derived type with allocatable memory by type: SSP from ssp_mod ! Scalar components ! Fixed arrays -!$TAF store ssp%c = setssp1 +!ML!$TAF store ssp%c = setssp1 ! Allocatable arrays -!$TAF store ssp%cmat,ssp%czmat,ssp%seg%r = setssp1 +!ML!$TAF store ssp%cmat,ssp%czmat,ssp%seg%r = setssp1 +!$TAF STORE ssp%c = comlev1_ihop, key = ikey_dynamics +!$TAF STORE ssp%cz, ssp%rho = comlev1_ihop, key = ikey_dynamics ! init defaults for ssp_mod scoped arrays n2 = (-1.,-1.) @@ -844,9 +849,12 @@ SUBROUTINE gcmSSP( myThid ) ! parameter in a different way after ssp_mod is split btwn fixed and varia REAL (KIND=_RL90) :: bPower, fT #ifdef ALLOW_AUTODIFF_TAMC - INTEGER tkey, ijkey, hkey, lockey + INTEGER tkey, ijkey, hkey, kkey +!ML INTEGER nloctape -!$TAF init loctape_ihop_gcmssp_bibj_ij_iijj_k = STATIC, nSx*nSy*sNx*sNy*IHOP_MAX_RANGE*IHOP_MAX_NC_SIZE*(Nr + 2) +!ML nloctape = nSx*nSy*sNx*sNy*IHOP_MAX_RANGE*IHOP_MAX_NC_SIZE +!ML!$TAF init loctape_ihop_iijj = STATIC, nloctape +!ML!$TAF init loctape_ihop_iijj_k = STATIC, nloctape*(Nr + 2) #endif ! IESCO24 fT init @@ -877,10 +885,17 @@ SUBROUTINE gcmSSP( myThid ) tmpSSP = 0.0 _d 0 globSSP = 0.0 _d 0 +#ifdef ALLOW_AUTODIFF_TAMC + ! This is weird. The alternative is to copy SPP%Nr and SPP%Nz to + ! local variables. +!$TAF INIT loctape_ihop_ssp = COMMON, 1 +!$TAF STORE ssp%nr, ssp%nz = loctape_ihop_ssp, key = 1 +#endif ! interpolate SSP with adaptive IDW from gcm grid to ihop grid DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) #ifdef ALLOW_AUTODIFF_TAMC +!ML tkey = bi + (bj-1)*nSx tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy #endif @@ -889,13 +904,15 @@ SUBROUTINE gcmSSP( myThid ) DO i=1,sNx #ifdef ALLOW_AUTODIFF_TAMC ijkey = i + (j-1)*sNx + (tkey-1)*sNx*sNy -!$TAF store njj(ii) = comlev1_bibj_ij_ihop, key=ijkey, kind=isbyte #endif DO ii=1,IHOP_npts_range interp_finished = .FALSE. DO jj=1,IHOP_npts_idw #ifdef ALLOW_AUTODIFF_TAMC -!$TAF STORE interp_finished = comlev1_bibj_ij_ihop, key=ijkey, kind=isbyte + hkey = jj + (ii-1)*IHOP_npts_idw & + + (ijkey-1)*IHOP_npts_idw*IHOP_npts_range +!ML!$TAF STORE interp_finished = loctape_ihop_iijj, key=hkey +!$TAF STORE interp_finished = comlev1_ihop_iijj, key=hkey #endif ! Interpolate from GCM grid cell centers @@ -904,14 +921,11 @@ SUBROUTINE gcmSSP( myThid ) .NOT. interp_finished) THEN njj(ii) = njj(ii) + 1 -#ifdef ALLOW_AUTODIFF_TAMC -#endif DO iz = 1, SSP%Nz - 1 #ifdef ALLOW_AUTODIFF_TAMC - hkey = jj + (ii-1)*IHOP_npts_idw + (ijkey-1)*sNy*sNx*nSy*nSx - lockey = iz + (hkey-1)*(SSP%Nz-1)*IHOP_npts_idw*IHOP_npts_range*sNx*sNy*nSx*nSy -! + ((jj-1) + ((ii-1) + (ijkey-1)*IHOP_npts_range)*IHOP_npts_idw)*(SSP%Nz-1) -!$TAF store njj(ii) = loctape_ihop_gcmssp_bibj_ij_iijj_k, key=lockey, kind=isbyte + kkey = iz + (hkey-1)*(SSP%Nz-1) +!ML!$TAF store njj(ii) = loctape_ihop_iijj_k, key=kkey +!$TAF store njj(ii) = comlev1_ihop_iijj_k, key=kkey #endif IF (iz .EQ. 1) THEN @@ -922,7 +936,7 @@ SUBROUTINE gcmSSP( myThid ) ELSE ! 2:(SSP%Nz-1) ! Middle depth layers, only when not already underground IF (ihop_sumweights(ii, iz-1) .GT. 0.0) THEN - ! isolate njj increments for TAF + ! isolate njj increments for TAF IF (ihop_idw_weights(ii, jj) .EQ. 0.0) THEN njj(ii) = IHOP_npts_idw + 1 @@ -990,6 +1004,9 @@ SUBROUTINE gcmSSP( myThid ) END DO END DO END DO +#ifdef ALLOW_AUTODIFF_TAMC +!$TAF STORE tmpssp = comlev1_ihop, key = ikey_dynamics +#endif !IESCO c68s: IF ((nPx.GT.1) .OR. (nPy.GT.1)) THEN !IESCO c68s: CALL GLOBAL_VEC_SUM_R8(SSP%Nz*SSP%Nr,SSP%Nz*SSP%Nr,tileSSP,myThid) @@ -1004,6 +1021,9 @@ SUBROUTINE gcmSSP( myThid ) k = k + 1 END DO END DO +#ifdef ALLOW_AUTODIFF_TAMC +!$TAF STORE ssp%cMat = comlev1_ihop, key = ikey_dynamics +#endif ! IESCO24: END MITgcm checkpoint69a uses a new global sum subroutine... IF(ALLOCATED(tileSSP)) DEALLOCATE(tileSSP)