--- MITgcm/pkg/ecco/the_main_loop.F 2004/11/13 06:19:10 1.14 +++ MITgcm/pkg/ecco/the_main_loop.F 2006/12/13 22:31:19 1.48 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ecco/Attic/the_main_loop.F,v 1.14 2004/11/13 06:19:10 heimbach Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ecco/Attic/the_main_loop.F,v 1.48 2006/12/13 22:31:19 heimbach Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" @@ -72,6 +72,26 @@ #include "EEPARAMS.h" #include "PARAMS.h" +#ifdef ALLOW_MNC +#include "MNC_PARAMS.h" + EXTERNAL DIFFERENT_MULTIPLE + LOGICAL DIFFERENT_MULTIPLE +#endif + +#ifdef HAVE_SIGREG +#include "SIGREG.h" +#endif + +#ifdef ALLOW_SHAP_FILT +# include "SHAP_FILT.h" +#endif +#ifdef ALLOW_ZONAL_FILT +# include "ZONAL_FILT.h" +#endif +#ifdef COMPONENT_MODULE +# include "CPL_PARAMS.h" +#endif + c************************************** #ifdef ALLOW_AUTODIFF_TAMC @@ -83,8 +103,10 @@ # include "DYNVARS.h" # include "FFIELDS.h" # include "EOS.h" -# include "GAD.h" +# ifdef ALLOW_GENERIC_ADVDIFF +# include "GAD.h" +# endif # ifdef ALLOW_CD_CODE # include "CD_CODE_VARS.h" # endif @@ -95,7 +117,7 @@ # ifdef ALLOW_NONHYDROSTATIC # include "CG3D.h" # endif -# ifdef EXACT_CONSERV +# if (defined (EXACT_CONSERV) || defined (NONLIN_FRSURF)) # include "SURFACE.h" # endif # ifdef ALLOW_OBCS @@ -110,6 +132,7 @@ # endif /* ALLOW_EXF */ # ifdef ALLOW_SEAICE # include "SEAICE.h" +# include "SEAICE_COST.h" # endif # ifdef ALLOW_KPP # include "KPP.h" @@ -117,6 +140,12 @@ # ifdef ALLOW_GMREDI # include "GMREDI.h" # endif +# ifdef ALLOW_RBCS +# include "RBCS.h" +# endif +# ifdef ALLOW_PROFILES +# include "profiles.h" +# endif # ifdef ALLOW_DIVIDED_ADJOINT_MPI # include "mpif.h" # endif @@ -156,10 +185,45 @@ integer ilev_1 integer ilev_2 integer ilev_3 + integer ilev_4 integer max_lev2 integer max_lev3 + integer max_lev4 #endif +#if defined(TIME_PER_TIMESTEP) || defined(USE_PAPI_FLOPS) || defined(USE_PCL_FLOPS) + CHARACTER*(MAX_LEN_MBUF) msgBuf +#ifdef TIME_PER_TIMESTEP +CCE107 common block for per timestep timing +C !TIMING VARIABLES +C == Timing variables == + REAL*8 utnew, utold, stnew, stold, wtnew, wtold + DATA utnew, utold, stnew, stold, wtnew, wtold /6*0.0D0/ +#endif +#ifdef USE_PAPI_FLOPS +CCE107 common block for PAPI summary performance +#include + INTEGER*8 flpops, instr + DATA flpops, instr /2*0/ + INTEGER check + REAL*4 real_time, proc_time, mflops, ipc + DATA real_time, proc_time, mflops, ipc /4*0.0E0/ +#else +#ifdef USE_PCL_FLOPS +CCE107 common block for PCL summary performance +#include + INTEGER pcl_counter_list(5), flags, nevents, res, ipcl + INTEGER*8 i_result(5), descr + REAL*8 fp_result(5) + COMMON /pclvars/ i_result, descr, fp_result, pcl_counter_list, + $ flags, nevents + INTEGER nmaxevents + PARAMETER (nmaxevents = 61) + CHARACTER*22 pcl_counter_name(0:nmaxevents-1) + COMMON /pclnames/ pcl_counter_name +#endif +#endif +#endif c-- == end of interface == #ifndef DISABLE_DEBUGMODE @@ -171,11 +235,15 @@ c-- Initialize storage for the initialisations. CADJ INIT tapelev_ini_bibj_k = USER CADJ INIT tapelev_init = USER -#ifdef AUTODIFF_2_LEVEL_CHECKPOINT +c +#if (defined (AUTODIFF_2_LEVEL_CHECKPOINT)) CADJ INIT tapelev2 = USER +#elif (defined (AUTODIFF_4_LEVEL_CHECKPOINT)) +CADJ INIT tapelev4 = USER #else CADJ INIT tapelev3 = USER #endif +c # ifdef ALLOW_DIVIDED_ADJOINT CADJ INIT onetape = user cphCADJ INIT onetape = common, 1 @@ -193,8 +261,16 @@ CALL TIMER_START('ECCO SPIN-UP', mythid) +#ifdef ALLOW_CAL c-- Get the current date. call CAL_TIMESTAMP( myiter, mytime, mydate, mythid ) +#endif + +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +CADJ STORE hFacC = tapelev_init, key = 1 +# endif +#endif C-- Set initial conditions (variable arrays) #ifndef DISABLE_DEBUGMODE @@ -233,6 +309,48 @@ call timer_stop ('ECCO SPIN-UP', mythid) _BARRIER +#ifdef TIME_PER_TIMESTEP +CCE107 Initial call for timers + _BEGIN_MASTER( myThid ) + CALL TIMER_GET_TIME( utold, stold, wtold ) + _END_MASTER( myThid ) +#endif +#ifdef USE_PAPI_FLOPS +CCE107 Initial call for PAPI + _BEGIN_MASTER( myThid ) +#ifdef USE_FLIPS + call PAPIF_flips(real_time, proc_time, flpops, mflops, check) +#else + call PAPIF_flops(real_time, proc_time, flpops, mflops, check) +#endif + WRITE(msgBuf,'(A34,F10.6,A,F10.6)') + $ 'Mflop/s before timestepping:', mflops, ' ', mflops*proc_time + $ /(real_time + 1E-36) + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) +#ifdef PAPI_VERSION + call PAPIF_ipc(real_time, proc_time, instr, ipc, check) + WRITE(msgBuf,'(A34,F10.6,A,F10.6)') + $ 'IPC before timestepping:', ipc, ' ', ipc*proc_time + $ /(real_time + 1E-36) + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) +#endif + _END_MASTER( myThid ) +#else +#ifdef USE_PCL_FLOPS +CCE107 Initial call for PCL + _BEGIN_MASTER( myThid ) + res = PCLstop(descr, i_result, fp_result, nevents) + do ipcl = 1, nevents + WRITE(msgBuf,'(A22,A26,F10.6)'), + $ pcl_counter_name(pcl_counter_list(ipcl)), + $ 'before timestepping:', fp_result(ipcl) + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) + enddo + res = PCLstart(descr, pcl_counter_list, nevents, flags) + _END_MASTER( myThid ) +#endif +#endif + c-- Do the model integration. call timer_start('ECCO MAIN LOOP',mythid) @@ -242,6 +360,7 @@ #ifdef ALLOW_AUTODIFF_TAMC #ifdef ALLOW_TAMC_CHECKPOINTING + max_lev4=nTimeSteps/(nchklev_1*nchklev_2*nchklev_3)+1 max_lev3=nTimeSteps/(nchklev_1*nchklev_2)+1 max_lev2=nTimeSteps/nchklev_1+1 @@ -251,17 +370,24 @@ #endif c************************************** -#ifndef AUTODIFF_2_LEVEL_CHECKPOINT +#ifdef AUTODIFF_4_LEVEL_CHECKPOINT + do ilev_4 = 1,nchklev_4 + if(ilev_4.le.max_lev4) then +c************************************** +#include "checkpoint_lev4_directives.h" +c************************************** +c-- Initialise storage for the middle loop. +CADJ INIT tapelev3 = USER +#endif /* AUTODIFF_4_LEVEL_CHECKPOINT */ +#ifndef AUTODIFF_2_LEVEL_CHECKPOINT do ilev_3 = 1,nchklev_3 if(ilev_3.le.max_lev3) then c************************************** #include "checkpoint_lev3_directives.h" c************************************** - c-- Initialise storage for the middle loop. CADJ INIT tapelev2 = USER - #endif /* AUTODIFF_2_LEVEL_CHECKPOINT */ do ilev_2 = 1,nchklev_2 @@ -292,14 +418,17 @@ # ifdef ALLOW_PTRACERS CADJ INIT comlev1_bibj_ptracers = COMMON, CADJ & nchklev_1*nsx*nsy*nthreads_chkpt*PTRACERS_num +CADJ INIT comlev1_bibj_k_ptracers = COMMON, +CADJ & nchklev_1*nsx*nsy*nthreads_chkpt*PTRACERS_num*nr # endif /* ALLOW_PTRACERS */ c-- -# ifndef DISABLE_MULTIDIM_ADVECTION +cph Now also needed by seaice +cph# ifndef DISABLE_MULTIDIM_ADVECTION CADJ INIT comlev1_bibj_k_gad CADJ & = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt*maxpass CADJ INIT comlev1_bibj_k_gad_pass CADJ & = COMMON,nchklev_1*nsx*nsy*nr*nthreads_chkpt*maxpass*maxcube -# endif /* DISABLE_MULTIDIM_ADVECTION */ +cph# endif /* DISABLE_MULTIDIM_ADVECTION */ c-- # if (defined (ALLOW_EXF) && defined (ALLOW_BULKFORMULAE)) CADJ INIT comlev1_exf_1 @@ -328,10 +457,13 @@ c-- c-- Requirement: nchklev_1*nchklev_2*nchklev_3 .ge. nTimeSteps . - iloop = (ilev_2 - 1)*nchklev_1 + ilev_1 + iloop = (ilev_2 - 1)*nchklev_1 + ilev_1 #ifndef AUTODIFF_2_LEVEL_CHECKPOINT & + (ilev_3 - 1)*nchklev_2*nchklev_1 #endif +#ifdef AUTODIFF_4_LEVEL_CHECKPOINT + & + (ilev_4 - 1)*nchklev_3*nchklev_2*nchklev_1 +#endif if ( iloop .le. nTimeSteps ) then @@ -366,7 +498,7 @@ #endif /* ALLOW_AUTODIFF_TAMC */ #ifdef ALLOW_TAMC_CHECKPOINTING - nIter0 = INT( startTime/deltaTClock ) + nIter0 = NINT( (startTime-baseTime)/deltaTClock ) ikey_dynamics = ilev_1 #endif @@ -374,12 +506,36 @@ myiter = nIter0 + (iloop-1) mytime = startTime + float(iloop-1)*deltaTclock +#ifdef ALLOW_AUTODIFF_TAMC + CALL AUTODIFF_INADMODE_UNSET( myThid ) +#endif + +#ifdef ALLOW_DIAGNOSTICS +C-- State-variables diagnostics + IF ( useDiagnostics ) THEN +C-- Switch on/off diagnostics for snap-shot output: + CALL DIAGNOSTICS_SWITCH_ONOFF( myTime, myIter, myThid ) + CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + CALL DO_STATEVARS_DIAGS( myTime, 0, myIter, myThid ) + CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + ENDIF +#endif + +#ifdef ALLOW_PROFILES +c-- Accumulate in-situ time averages of temperature, salinity, and SSH. + call timer_start('PROFILES_INLOOP [ECCO MAIN]', mythid) + call profiles_inloop( mytime, mythid ) + call timer_stop ('PROFILES_INLOOP [ECCO MAIN]', mythid) +#endif + #ifdef ALLOW_COST -c-- Accumulate time averages of temperature, salinity, and SSH. +c-- Accumulate time averages of temperature, salinity call timer_start('COST_AVERAGESFIELDS [ECCO MAIN]', mythid) call cost_averagesFields( mytime, mythid ) call timer_stop ('COST_AVERAGESFIELDS [ECCO MAIN]', mythid) + + #ifdef ALLOW_COST_ATLANTIC c-- Compute meridional heat transport call timer_start('cost_atlantic [ECCO MAIN]', mythid) @@ -391,32 +547,25 @@ #ifdef ALLOW_AUTODIFF_TAMC c************************************** #include "checkpoint_lev1_directives.h" +#include "checkpoint_lev1_template.h" c************************************** #endif -#ifndef DISABLE_DEBUGMODE - IF ( debugLevel .GE. debLevB ) - & CALL DEBUG_CALL('EXF_GETFORCING',myThid) +C-- Call driver to load external forcing fields from file +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('LOAD_FIELDS_DRIVER',myThid) #endif - CALL TIMER_START('EXF_GETFORCING [FORWARD_STEP]',mythid) - CALL EXF_GETFORCING( mytime, myiter, mythid ) - CALL TIMER_STOP ('EXF_GETFORCING [FORWARD_STEP]',mythid) + CALL TIMER_START('LOAD_FIELDS_DRIVER [FORWARD_STEP]',myThid) + CALL LOAD_FIELDS_DRIVER( myTime, myIter, myThid ) + CALL TIMER_STOP ('LOAD_FIELDS_DRIVER [FORWARD_STEP]',myThid) -#ifdef ALLOW_SEAICE -cph this simple runtime flag causes a lot of recomp. -cph IF ( useSEAICE ) THEN -#ifndef DISABLE_DEBUGMODE - IF ( debugLevel .GE. debLevB ) - & CALL DEBUG_CALL('SEAICE_MODEL',myThid) -#endif - CALL TIMER_START('SEAICE_MODEL [FORWARD_STEP]',myThid) - CALL SEAICE_MODEL( myTime, myIter, myThid ) - CALL TIMER_STOP ('SEAICE_MODEL [FORWARD_STEP]',myThid) -#ifdef ALLOW_COST_ICE - CALL COST_ICE ( myTime, myIter, myThid ) + +#ifdef ALLOW_AUTODIFF_TAMC +# if (defined (ALLOW_AUTODIFF_MONITOR)) + CALL DUMMY_IN_STEPPING( myTime, myIter, myThid ) +# endif #endif -cph ENDIF -#endif /* ALLOW_SEAICE */ #ifdef ALLOW_AUTODIFF_TAMC # ifdef ALLOW_PTRACERS @@ -425,17 +574,6 @@ # endif #endif -#if (defined (ALLOW_AUTODIFF_TAMC) && \ - defined (ALLOW_AUTODIFF_MONITOR)) -C Include call to a dummy routine. Its adjoint will be -C called at the proper place in the adjoint code. -C The adjoint routine will print out adjoint values -C if requested. The location of the call is important, -C it has to be after the adjoint of the exchanges -C (DO_GTERM_BLOCKING_EXCHANGES). - CALL DUMMY_IN_STEPPING( myTime, myIter, myThid ) -#endif - #ifdef ALLOW_EBM IF ( useEBM ) THEN # ifdef ALLOW_DEBUG @@ -458,6 +596,20 @@ CALL DO_ATMOSPHERIC_PHYS( myTime, myIter, myThid ) CALL TIMER_STOP ('DO_ATMOSPHERIC_PHYS [FORWARD_STEP]',mythid) +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef EXACT_CONSERV +cphCADJ STORE empmr = comlev1, key = ikey_dynamics +cphCADJ STORE pmepr = comlev1, key = ikey_dynamics +# endif +# ifdef ALLOW_PTRACERS +CADJ STORE ptracer = comlev1, key = ikey_dynamics +# endif +# ifdef NONLIN_FRSURF +cph-test +CADJ STORE hFacC = comlev1, key = ikey_dynamics +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + #ifndef ALLOW_OFFLINE #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) @@ -466,9 +618,60 @@ CALL TIMER_START('DO_OCEANIC_PHYS [FORWARD_STEP]',mythid) CALL DO_OCEANIC_PHYS( myTime, myIter, myThid ) CALL TIMER_STOP ('DO_OCEANIC_PHYS [FORWARD_STEP]',mythid) +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE EmPmR = comlev1, key = ikey_dynamics +# ifdef EXACT_CONSERV +CADJ STORE pmepr = comlev1, key = ikey_dynamics +# endif +#endif #endif #ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +cph-test +CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics +CADJ STORE hfac_surfs = comlev1, key = ikey_dynamics +CADJ STORE hfac_surfw = comlev1, key = ikey_dynamics +CADJ STORE hFacC, hFacS, hFacW +CADJ & = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacC, recip_hFacS, recip_hFacW +CADJ & = comlev1, key = ikey_dynamics +c +CADJ STORE surfaceforcingu = comlev1, key = ikey_dynamics +CADJ STORE surfaceforcingv = comlev1, key = ikey_dynamics +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +#ifdef ALLOW_GCHEM +C GCHEM package is an interface for any bio-geochemical or +C ecosystem model you would like to include. +C If GCHEM_SEPARATE_FORCING is not defined, you are +C responsible for computing tendency terms for passive +C tracers and storing them on a 3DxNumPtracers-array called +C gchemTendency in GCHEM_CALC_TENDENCY. This tendency is then added +C to gPtr in ptracers_forcing later-on. +C If GCHEM_SEPARATE_FORCING is defined, you are reponsible for +C UPDATING ptracers directly in GCHEM_FORCING_SEP. This amounts +C to a completely separate time step that you have to implement +C yourself (Eulerian seems to be fine in most cases). +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC +C CAVEAT: Up to now, when GCHEM is turned on the field ptracerForcingSurf, +C which is needed for KPP is not set properly. ptracerForcingSurf must +C be treated differently depending on whether GCHEM_SEPARATE_FORCING +C is define or not. TBD. +CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC + IF ( useGCHEM ) THEN +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GCHEM_CALC_TENDENCY',myThid) +#endif + CALL TIMER_START('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid) + CALL GCHEM_CALC_TENDENCY( myTime, myIter, myThid ) + CALL TIMER_STOP ('GCHEM_CALC_TENDENCY [FORWARD_STEP]',myThid) + ENDIF +#endif /* ALLOW_GCHEM */ + +#ifdef ALLOW_AUTODIFF_TAMC cph needed to be moved here from do_oceanic_physics cph to be visible down the road c @@ -500,6 +703,14 @@ # endif #endif /* ALLOW_AUTODIFF_TAMC */ +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +CADJ STORE etaH = comlev1, key = ikey_dynamics +# ifdef ALLOW_CD_CODE +CADJ STORE etanm1 = comlev1, key = ikey_dynamics +# endif +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ IF ( .NOT.staggerTimeStep ) THEN #ifdef ALLOW_DEBUG @@ -512,6 +723,18 @@ C-- if not staggerTimeStep: end ENDIF +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +CADJ STORE hFacC = comlev1, key = ikey_dynamics +CADJ STORE hFacS = comlev1, key = ikey_dynamics +CADJ STORE hFacW = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics +CADJ STORE etaN = comlev1, key = ikey_dynamics +# endif +#endif + C-- Step forward fields and calculate time tendency terms. #ifndef ALLOW_OFFLINE #ifndef ALLOW_AUTODIFF_TAMC @@ -529,44 +752,92 @@ #endif #endif -#ifdef ALLOW_NONHYDROSTATIC -C-- Step forward W field in N-H algorithm - IF ( momStepping .AND. nonHydrostatic ) THEN -#ifdef ALLOW_DEBUG - IF ( debugLevel .GE. debLevB ) - & CALL DEBUG_CALL('CALC_GW',myThid) -#endif - CALL TIMER_START('CALC_GW [FORWARD_STEP]',myThid) - CALL CALC_GW(myThid) - CALL TIMER_STOP ('CALC_GW [FORWARD_STEP]',myThid) - ENDIF +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +cph-test +CADJ STORE gU, gV = comlev1, key = ikey_dynamics +# endif #endif C-- Update time-counter myIter = nIter0 + iLoop myTime = startTime + deltaTClock * float(iLoop) +#ifdef ALLOW_MNC +C Update the default next iter for MNC + IF ( useMNC ) THEN + CALL MNC_CW_CITER_SETG( 1, 1, -1, myIter , myThid ) + +C TODO: Logic should be added here so that users can specify, on +C a per-citer-group basis, when it is time to update the +C "current" (and not just the "next") iteration + +C TODO: the following is just a temporary band-aid (mostly, for +C Baylor) until someone writes a routine that better handles time +C boundaries such as weeks, months, years, etc. + IF ( mnc_filefreq .GT. 0 ) THEN + IF (DIFFERENT_MULTIPLE(mnc_filefreq,myTime,deltaTClock)) + & THEN + CALL MNC_CW_CITER_SETG( 1, 1, myIter, -1 , myThid ) + ENDIF + ENDIF + ENDIF +#endif + C-- Update geometric factors: #ifdef NONLIN_FRSURF C- update hfacC,W,S and recip_hFac according to etaH(n+1) : IF ( nonlinFreeSurf.GT.0) THEN IF ( select_rStar.GT.0 ) THEN +# ifndef DISABLE_RSTAR_CODE +# ifdef ALLOW_AUTODIFF_TAMC +cph-test +CADJ STORE hFacC = comlev1, key = ikey_dynamics +CADJ STORE hFacS = comlev1, key = ikey_dynamics +CADJ STORE hFacW = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics +# endif CALL TIMER_START('UPDATE_R_STAR [FORWARD_STEP]',myThid) CALL UPDATE_R_STAR( myTime, myIter, myThid ) CALL TIMER_STOP ('UPDATE_R_STAR [FORWARD_STEP]',myThid) +# ifdef ALLOW_AUTODIFF_TAMC +cph-test +CADJ STORE hFacC = comlev1, key = ikey_dynamics +CADJ STORE hFacS = comlev1, key = ikey_dynamics +CADJ STORE hFacW = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics +# endif +# endif /* DISABLE_RSTAR_CODE */ ELSE +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE hFac_surfC, hFac_surfS, hFac_surfW +CADJ & = comlev1, key = ikey_dynamics +#endif CALL TIMER_START('UPDATE_SURF_DR [FORWARD_STEP]',myThid) CALL UPDATE_SURF_DR( myTime, myIter, myThid ) CALL TIMER_STOP ('UPDATE_SURF_DR [FORWARD_STEP]',myThid) ENDIF ENDIF +# ifdef ALLOW_AUTODIFF_TAMC +cph-test +CADJ STORE hFacC = comlev1, key = ikey_dynamics +CADJ STORE hFacS = comlev1, key = ikey_dynamics +CADJ STORE hFacW = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacC = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacS = comlev1, key = ikey_dynamics +CADJ STORE recip_hFacW = comlev1, key = ikey_dynamics +# endif C- update also CG2D matrix (and preconditioner) IF ( momStepping .AND. nonlinFreeSurf.GT.2 ) THEN CALL TIMER_START('UPDATE_CG2D [FORWARD_STEP]',myThid) CALL UPDATE_CG2D( myTime, myIter, myThid ) CALL TIMER_STOP ('UPDATE_CG2D [FORWARD_STEP]',myThid) ENDIF -#endif +#endif /* NONLIN_FRSURF */ C-- Apply Filters to u*,v* before SOLVE_FOR_PRESSURE #ifdef ALLOW_SHAP_FILT @@ -600,6 +871,14 @@ C three-dimensional for non-hydrostatic and/or IGW scheme. #ifndef ALLOW_OFFLINE IF ( momStepping ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +# ifdef NONLIN_FRSURF +CADJ STORE uvel, vvel +CADJ & = comlev1, key = ikey_dynamics +CADJ STORE empmr,hfacs,hfacw +CADJ & = comlev1, key = ikey_dynamics +# endif +#endif CALL TIMER_START('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid) CALL SOLVE_FOR_PRESSURE(myTime, myIter, myThid) CALL TIMER_STOP ('SOLVE_FOR_PRESSURE [FORWARD_STEP]',myThid) @@ -617,6 +896,10 @@ #ifdef EXACT_CONSERV IF (exactConserv) THEN +#ifdef ALLOW_AUTODIFF_TAMC +cph-test +cphCADJ STORE etaH = comlev1, key = ikey_dynamics +#endif C-- Update etaH(n+1) : CALL TIMER_START('UPDATE_ETAH [FORWARD_STEP]',mythid) CALL UPDATE_ETAH( myTime, myIter, myThid ) @@ -626,16 +909,25 @@ #ifdef NONLIN_FRSURF IF ( select_rStar.NE.0 ) THEN +# ifndef DISABLE_RSTAR_CODE C-- r* : compute the future level thickness according to etaH(n+1) CALL TIMER_START('CALC_R_STAR [FORWARD_STEP]',mythid) CALL CALC_R_STAR(etaH, myTime, myIter, myThid ) CALL TIMER_STOP ('CALC_R_STAR [FORWARD_STEP]',mythid) +# endif /* DISABLE_RSTAR_CODE */ ELSEIF ( nonlinFreeSurf.GT.0) THEN C-- compute the future surface level thickness according to etaH(n+1) +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE etaH = comlev1, key = ikey_dynamics +# endif CALL TIMER_START('CALC_SURF_DR [FORWARD_STEP]',mythid) CALL CALC_SURF_DR(etaH, myTime, myIter, myThid ) CALL TIMER_STOP ('CALC_SURF_DR [FORWARD_STEP]',mythid) ENDIF +# ifdef ALLOW_AUTODIFF_TAMC +cph-test +CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics +# endif #endif /* NONLIN_FRSURF */ C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| @@ -649,6 +941,15 @@ CALL DO_STAGGER_FIELDS_EXCHANGES( myTime, myIter, myThid ) CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid) +#ifdef ALLOW_DIAGNOSTICS +C-- State-variables diagnostics + IF ( usediagnostics ) THEN + CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + CALL DO_STATEVARS_DIAGS( myTime, 1, myIter, myThid ) + CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + ENDIF +#endif + #ifdef ALLOW_DEBUG IF ( debugLevel .GE. debLevB ) & CALL DEBUG_CALL('THERMODYNAMICS',myThid) @@ -671,6 +972,20 @@ CALL TRACERS_CORRECTION_STEP(myTime, myIter, myThid) CALL TIMER_STOP ('TS_CORRECTION_STEP [FORWARD_STEP]',myThid) +#ifdef ALLOW_GCHEM +C Add separate timestepping of chemical/biological/forcing +C of ptracers here in GCHEM_FORCING_SEP + IF ( useGCHEM ) THEN +#ifdef ALLOW_DEBUG + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('GCHEM_FORCING_SEP',myThid) +#endif /* ALLOW_DEBUG */ + CALL TIMER_START('GCHEM_FORCING_SEP [FORWARD_STEP]',myThid) + CALL GCHEM_FORCING_SEP( myTime,myIter,myThid ) + CALL TIMER_STOP ('GCHEM_FORCING_SEP [FORWARD_STEP]',myThid) + ENDIF +#endif /* ALLOW_GCHEM */ + C-- Do "blocking" sends and receives for tendency "overlap" terms c CALL TIMER_START('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid) c CALL DO_GTERM_BLOCKING_EXCHANGES( myThid ) @@ -681,6 +996,14 @@ CALL DO_FIELDS_BLOCKING_EXCHANGES( myThid ) CALL TIMER_STOP ('BLOCKING_EXCHANGES [FORWARD_STEP]',myThid) +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + CALL DO_STATEVARS_DIAGS( myTime, 2, myIter, myThid ) + CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) + ENDIF +#endif + #ifdef ALLOW_FLT C-- Calculate float trajectories IF (useFLT) THEN @@ -690,10 +1013,16 @@ ENDIF #endif -C-- State-variables statistics (time-aver, diagnostics ...) - CALL TIMER_START('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) - CALL DO_STATEVARS_DIAGS( myTime, myIter, myThid ) - CALL TIMER_STOP ('DO_STATEVARS_DIAGS [FORWARD_STEP]',myThid) +#ifdef ALLOW_AUTODIFF_TAMC + CALL AUTODIFF_INADMODE_SET( myThid ) +#endif + +#ifdef ALLOW_TIMEAVE +C-- State-variables time-averaging + CALL TIMER_START('DO_STATEVARS_TAVE [FORWARD_STEP]',myThid) + CALL DO_STATEVARS_TAVE( myTime, myIter, myThid ) + CALL TIMER_STOP ('DO_STATEVARS_TAVE [FORWARD_STEP]',myThid) +#endif #ifndef ALLOW_OFFLINE #ifdef ALLOW_MONITOR @@ -715,31 +1044,94 @@ CALL TIMER_STOP ('DO_THE_MODEL_IO [FORWARD_STEP]',myThid) #endif +#ifndef ALLOW_DIVIDED_ADJOINT +# ifdef HAVE_SIGREG + IF ( useSIGREG ) THEN + IF ( i_got_signal .GT. 0 ) THEN + CALL PACKAGES_WRITE_PICKUP( + I .TRUE., myTime, myIter, myThid ) + CALL WRITE_PICKUP( + I .TRUE., myTime, myIter, myThid ) + STOP 'Checkpoint completed -- killed by signal handler' + ENDIF + ENDIF +# endif /* HAVE_SIGREG */ C-- Save state for restarts - CALL TIMER_START('WRITE_CHECKPOINT [FORWARD_STEP]',myThid) - CALL PACKAGES_WRITE_PICKUP( + CALL TIMER_START('DO_WRITE_PICKUP [FORWARD_STEP]',myThid) + CALL DO_WRITE_PICKUP( I .FALSE., myTime, myIter, myThid ) -#ifndef ALLOW_OFFLINE - CALL WRITE_CHECKPOINT( - I .FALSE., myTime, myIter, myThid ) + CALL TIMER_STOP ('DO_WRITE_PICKUP [FORWARD_STEP]',myThid) +#endif /* ALLOW_DIVIDED_ADJOINT */ + +#ifdef TIME_PER_TIMESTEP +CCE107 Time per timestep information + _BEGIN_MASTER( myThid ) + CALL TIMER_GET_TIME( utnew, stnew, wtnew ) + WRITE(msgBuf,'(A34,3F10.6,I8)') + $ 'User, system and wallclock time:', utnew - utold, + $ stnew - stold, wtnew - wtold, iloop + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) + utold = utnew + stold = stnew + wtold = wtnew + _END_MASTER( myThid ) +#endif +#ifdef USE_PAPI_FLOPS +CCE107 PAPI summary performance + _BEGIN_MASTER( myThid ) +#ifdef USE_FLIPS + call PAPIF_flips(real_time, proc_time, flpops, mflops, check) +#else + call PAPIF_flops(real_time, proc_time, flpops, mflops, check) +#endif + WRITE(msgBuf,'(F10.6,A,F10.6,A34,I8)') + $ mflops, ' ', mflops*proc_time/(real_time + 1E-36), + $ 'Mflop/s during timestep ', iloop + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) +#ifdef PAPI_VERSION + call PAPIF_ipc(real_time, proc_time, instr, ipc, check) + WRITE(msgBuf,'(F10.6,A,F10.6,A34,I8)') + $ ipc, ' ', ipc*proc_time/(real_time + 1E-36), + $ 'IPC during timestep ', iloop + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) +#endif + _END_MASTER( myThid ) +#else +#ifdef USE_PCL_FLOPS +CCE107 PCL summary performance + _BEGIN_MASTER( myThid ) + res = PCLstop(descr, i_result, fp_result, nevents) + do ipcl = 1, nevents + WRITE(msgBuf,'(F10.6,A2,A22,A17,I8)'), fp_result(ipcl), + $ ' ', pcl_counter_name(pcl_counter_list(ipcl)), + $ 'during timestep ', iloop + CALL PRINT_MESSAGE(msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) + enddo + res = PCLstart(descr, pcl_counter_list, nevents, flags) + _END_MASTER( myThid ) +#endif #endif - CALL TIMER_STOP ('WRITE_CHECKPOINT [FORWARD_STEP]',myThid) #ifdef ALLOW_AUTODIFF_TAMC #ifdef ALLOW_TAMC_CHECKPOINTING - endif + endif + enddo + endif enddo - endif - enddo #ifndef AUTODIFF_2_LEVEL_CHECKPOINT - endif - enddo + endif + enddo #endif -#else +#ifdef AUTODIFF_4_LEVEL_CHECKPOINT + endif enddo #endif +c +#else /* ndef ALLOW_TAMC_CHECKPOINTING */ + enddo +#endif /* ALLOW_TAMC_CHECKPOINTING */ -#else +#else /* ndef ALLOW_AUTODIFF_TAMC */ enddo #endif /* ALLOW_AUTODIFF_TAMC */ @@ -748,6 +1140,15 @@ call timer_start('ECCO SPIN-DOWN', mythid) +#ifdef ALLOW_PROFILES +#ifndef ALLOW_DIVIDED_ADJOINT +c-- Accumulate in-situ time averages of temperature, salinity, and SSH. + call timer_start('PROFILES_INLOOP [ECCO SPIN-DOWN]', mythid) + call profiles_inloop( mytime, mythid ) + call timer_stop ('PROFILES_INLOOP [ECCO SPIN-DOWN]', mythid) +#endif +#endif + #ifdef ALLOW_COST #ifdef ALLOW_DIVIDED_ADJOINT @@ -787,6 +1188,39 @@ call timer_start('cost_forcing [ECCO SPIN-DOWN]', mythid) call cost_forcing( myiter, mytime, mythid ) call timer_stop ('cost_forcing [ECCO SPIN-DOWN]', mythid) +cph( +c-- Compute cost function contribution of wind stress observations. +#ifdef ALLOW_MEAN_HFLUX_COST_CONTRIBUTION + call cost_mean_heatflux( myiter, mytime, mythid ) +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE objf_hfluxmm = tapelev_init, key = 1 +# endif +#endif + +c-- Compute cost function contribution of wind stress observations. +#ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION + call cost_mean_saltflux( myiter, mytime, mythid ) +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE objf_sfluxmm = tapelev_init, key = 1 +# endif +#endif +cph) + +c-- Compute cost function contribution of SSH. +#ifdef ALLOW_SSH_COST_CONTRIBUTION +# ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('cost_ssh',myThid) +# endif + call timer_start('cost_ssh [ECCO SPIN-DOWN]', mythid) + call cost_ssh( myiter, mytime, mythid ) + call timer_stop ('cost_ssh [ECCO SPIN-DOWN]', mythid) +# ifdef ALLOW_AUTODIFF_TAMC +# ifdef ALLOW_PROFILES +CADJ STORE prof_etan_mean = tapelev_init, key = 1 +# endif +# endif +#endif c-- Compute cost function contribution of Temperature and Salinity. #ifndef DISABLE_DEBUGMODE @@ -797,6 +1231,16 @@ call cost_hyd( myiter, mytime, mythid ) call timer_stop ('cost_hyd [ECCO SPIN-DOWN]', mythid) +#ifdef ALLOW_SEAICE +#ifndef DISABLE_DEBUGMODE + IF ( debugLevel .GE. debLevB ) + & CALL DEBUG_CALL('seaice_cost_driver',myThid) +#endif + call timer_start('seaice_cost_driver [ECCO SPIN-DOWN]', mythid) + call seaice_cost_driver( myiter, mytime, mythid ) + call timer_stop ('seaice_cost_driver [ECCO SPIN-DOWN]', mythid) +#endif + #ifdef ALLOW_OBCS_COST_CONTRIBUTION #ifndef DISABLE_DEBUGMODE IF ( debugLevel .GE. debLevB ) @@ -817,17 +1261,6 @@ call timer_stop ('cost_curmtr [ECCO SPIN-DOWN]', mythid) #endif -c-- Compute cost function contribution of SSH. -#ifdef ALLOW_SSH_COST_CONTRIBUTION -#ifndef DISABLE_DEBUGMODE - IF ( debugLevel .GE. debLevB ) - & CALL DEBUG_CALL('cost_ssh',myThid) -#endif - call timer_start('cost_ssh [ECCO SPIN-DOWN]', mythid) - call cost_ssh( myiter, mytime, mythid ) - call timer_stop ('cost_ssh [ECCO SPIN-DOWN]', mythid) -#endif - c-- Compute cost function contribution of drifter's velocities. #ifdef ALLOW_DRIFTER_COST_CONTRIBUTION #ifndef DISABLE_DEBUGMODE @@ -850,20 +1283,6 @@ call timer_stop ('cost_scat [ECCO SPIN-DOWN]', mythid) #endif -c-- Compute cost function contribution of wind stress observations. -#ifdef ALLOW_MEAN_HFLUX_COST_CONTRIBUTION - call timer_start('cost_mean_heatflux [ECCO SPIN-DOWN]', mythid) - call cost_mean_heatflux( myiter, mytime, mythid ) - call timer_stop ('cost_mean_heatflux [ECCO SPIN-DOWN]', mythid) -#endif - -c-- Compute cost function contribution of wind stress observations. -#ifdef ALLOW_MEAN_SFLUX_COST_CONTRIBUTION - call timer_start('cost_mean_saltflux [ECCO SPIN-DOWN]', mythid) - call cost_mean_saltflux( myiter, mytime, mythid ) - call timer_stop ('cost_mean_saltflux [ECCO SPIN-DOWN]', mythid) -#endif - c-- Compute cost function contribution of drift between the first c and the last year. #ifdef ALLOW_DRIFT_COST_CONTRIBUTION @@ -895,13 +1314,20 @@ #endif _BARRIER +c-- Eddy stress penalty term +#ifdef ALLOW_COST_TAU_EDDY + call timer_start('cost_tau_eddy [ECCO SPIN-DOWN]', mythid) + call cost_tau_eddy( mythid ) + call timer_stop ('cost_tau_eddy [ECCO SPIN-DOWN]', mythid) +#endif + c-- Sum all cost function contributions. #ifndef DISABLE_DEBUGMODE IF ( debugLevel .GE. debLevB ) & CALL DEBUG_CALL('cost_final',myThid) #endif call timer_start('COST_FINAL [ECCO SPIN-DOWN]', mythid) - call ecco_cost_final( mythid ) + call cost_final( mythid ) call timer_stop ('COST_FINAL [ECCO SPIN-DOWN]', mythid) #endif /* ALLOW_COST */