Index: /issm/trunk-jpl/test/MITgcm/code/EESUPPORT.h =================================================================== --- /issm/trunk-jpl/test/MITgcm/code/EESUPPORT.h (revision 22532) +++ /issm/trunk-jpl/test/MITgcm/code/EESUPPORT.h (revision 22532) @@ -0,0 +1,298 @@ +C $Header: /u/gcmpack/MITgcm/eesupp/inc/EESUPPORT.h,v 1.10 2009/04/21 16:00:53 jmc Exp $ +C $Name: $ +CBOP +C !ROUTINE: EESUPPORT.h +C !INTERFACE: +C include "EESUPPORT.h" +C +C !DESCRIPTION: +C *==========================================================* +C | EESUPPORT.h | +C *==========================================================* +C | Support data structures for the MITgcm UV execution | +C | environment code. This data should be private to the | +C | execution environment routines. Data which needs to be | +C | accessed directly by a numerical model goes in | +C | EEPARAMS.h. | +C *==========================================================* +CEOP + +C ERROR_HEADER - String which prefixes error messages + CHARACTER*(*) ERROR_HEADER + PARAMETER ( ERROR_HEADER = ' *** ERROR ***' ) +C PROCESS_HEADER - String which prefixes processor number + CHARACTER*(*) PROCESS_HEADER + PARAMETER ( PROCESS_HEADER = 'PID.TID' ) + +C MAX_NUM_COMM_MODES - Maximum number of communication modes +C COMM_NONE - No edge communication +C COMM_MSG - Use messages to communicate edges +C COMM_PUT - Use put to communicate edges +C COMM_GET - Use get to communicate edges +C Note - commName holds an identifying name for each communication +C mode. The COMM_ parameters are used to index commName +C so the COMM_ parameters need to be in the range +C 1 : MAX_NUM_COMM_MODES. + INTEGER MAX_NUM_COMM_MODES + PARAMETER ( MAX_NUM_COMM_MODES = 4 ) + INTEGER COMM_NONE + PARAMETER ( COMM_NONE = 1 ) + INTEGER COMM_MSG + PARAMETER ( COMM_MSG = 2 ) + INTEGER COMM_PUT + PARAMETER ( COMM_PUT = 3 ) + INTEGER COMM_GET + PARAMETER ( COMM_GET = 4 ) + COMMON /EESUPP_COMMNAME/ commName + CHARACTER*10 commName(MAX_NUM_COMM_MODES) + +C Tile identifiers +C Tiles have a number that is unique over the global domain. +C A tile that is not there has its number set to NULL_TILE + INTEGER NULL_TILE + PARAMETER ( NULL_TILE = -1 ) + + +C-- COMMON /EESUPP_C/ Execution environment support character variables +C myProcessStr - String identifying my process number + COMMON /EESUPP_C/ myProcessStr + CHARACTER*128 myProcessStr + +C-- COMMON /EESUPP_L/ Execution environment support logical variables +C initMPError - Flag indicating error during multi-processing +C initialisation. +C finMPError - Flag indicating error during multi-processing +C termination. +C ThError - Thread detected an error. +C usingMPI - Flag controlling use of MPI routines. This flag +C allows either MPI or threads to be used in a +C shared memory environment which can be a useful +C debugging/performance analysis tool. +C usingSyncMessages - Flag that causes blocking communication to be used +C if possible. When false non-blocking EXCH routines +C will be used if possible. +C notUsingXPeriodicity - Flag indicating no X/Y boundary wrap around +C notUsingYPeriodicity This affects the communication routines but +C is generally ignored in the numerical model +C code. +C threadIsRunning, threadIsComplete - Flags used to check for correct behaviour +C of multi-threaded code. +C threadIsRunning is used to check that the +C threads we need are running. This catches the +C situation where a program eedata file has nTthreads +C greater than the setenv PARALLEL or NCPUS variable. +C threadIsComplete is used to flag that a thread has +C reached the end of the model. This is used as a check to +C trap problems that might occur if one thread "escapes" +C the main.F master loop. This should not happen +C if the multi-threading compilation tools works right. +C But (see for example KAP) this is not always the case! + COMMON /EESUPP_L/ thError, threadIsRunning, threadIsComplete, + & allMyEdgesAreSharedMemory, usingMPI, usingSyncMessages, + & notUsingXPeriodicity, notUsingYPeriodicity + LOGICAL thError(MAX_NO_THREADS) + LOGICAL threadIsRunning(MAX_NO_THREADS) + LOGICAL threadIsComplete(MAX_NO_THREADS) + LOGICAL allMyEdgesAreSharedMemory(MAX_NO_THREADS) + LOGICAL usingMPI + LOGICAL usingSyncMessages + LOGICAL notUsingXPeriodicity + LOGICAL notUsingYPeriodicity + +C-- COMMON /EESUPP_I/ Parallel support integer globals +C pidW - Process ID of neighbor to West +C pidE - ditto East +C pidN - ditto North +C pidS - ditto South +C Note: pid[XY] is not necessairily the UNIX +C process id - it is just an identifying +C number. +C myPid - My own process id +C nProcs - Number of processes +C westCommunicationMode - Mode of communication for each tile face +C eastCommunicationMode +C northCommunicationMode +C southCommunicationMode +C bi0 - Low cartesian tile index for this process +C bj0 Note - In a tile distribution with holes bi0 and bj0 +C are not useful. Neighboring tile indices must +C be derived some other way. +C tileNo - Tile identification number for my tile and +C tileNo[WENS] my N,S,E,W neighbor tiles. +C tilePid[WENS] - Process identification number for +C my N,S,E,W neighbor tiles. +C nTx, nTy - No. threads in X and Y. This assumes a simple +C cartesian gridding of the threads which is not +C required elsewhere but that makes it easier. + COMMON /EESUPP_I/ + & myPid, nProcs, pidW, pidE, pidN, pidS, + & tileCommModeW, tileCommModeE, + & tileCommModeN, tileCommModeS, + & tileNo, tileNoW, tileNoE, tileNoS, tileNoN, + & tilePidW, tilePidE, tilePidS, tilePidN, + & tileBiW, tileBiE, tileBiS, tileBiN, + & tileBjW, tileBjE, tileBjS, tileBjN, + & tileTagSendW, tileTagSendE, tileTagSendS, tileTagSendN, + & tileTagRecvW, tileTagRecvE, tileTagRecvS, tileTagRecvN + INTEGER myPid + INTEGER nProcs + INTEGER pidW + INTEGER pidE + INTEGER pidN + INTEGER pidS + INTEGER tileCommModeW ( nSx, nSy ) + INTEGER tileCommModeE ( nSx, nSy ) + INTEGER tileCommModeN ( nSx, nSy ) + INTEGER tileCommModeS ( nSx, nSy ) + INTEGER tileNo( nSx, nSy ) + INTEGER tileNoW( nSx, nSy ) + INTEGER tileNoE( nSx, nSy ) + INTEGER tileNoN( nSx, nSy ) + INTEGER tileNoS( nSx, nSy ) + INTEGER tilePidW( nSx, nSy ) + INTEGER tilePidE( nSx, nSy ) + INTEGER tilePidN( nSx, nSy ) + INTEGER tilePidS( nSx, nSy ) + INTEGER tileBiW( nSx, nSy ) + INTEGER tileBiE( nSx, nSy ) + INTEGER tileBiN( nSx, nSy ) + INTEGER tileBiS( nSx, nSy ) + INTEGER tileBjW( nSx, nSy ) + INTEGER tileBjE( nSx, nSy ) + INTEGER tileBjN( nSx, nSy ) + INTEGER tileBjS( nSx, nSy ) + INTEGER tileTagSendW( nSx, nSy ) + INTEGER tileTagSendE( nSx, nSy ) + INTEGER tileTagSendN( nSx, nSy ) + INTEGER tileTagSendS( nSx, nSy ) + INTEGER tileTagRecvW( nSx, nSy ) + INTEGER tileTagRecvE( nSx, nSy ) + INTEGER tileTagRecvN( nSx, nSy ) + INTEGER tileTagRecvS( nSx, nSy ) + +#ifdef ALLOW_USE_MPI +C-- Include MPI standard Fortran header file +#include "mpif.h" +#define _mpiTRUE_ 1 +#define _mpiFALSE_ 0 + +C-- COMMON /EESUPP_MPI_I/ MPI parallel support integer globals +C mpiPidW - MPI process id for west neighbor. +C mpiPidE - MPI process id for east neighbor. +C mpiPidN - MPI process id for north neighbor. +C mpiPidS - MPI process id for south neighbor. +C mpiPidNW - MPI process id for northwest neighbor. +C mpiPidNE - MPI process id for northeast neighbor. +C mpiPidSW - MPI process id for southwest neighbor. +C mpiPidSE - MPI process id for southeast neighbor. +C mpiPidIO - MPI process to use for IO. +C mpiNprocs - No. of MPI processes. +C mpiMyId - MPI process id of me. +C mpiComm - MPI communicator to use. +C mpiPx - My MPI proc. grid X coord +C mpiPy - My MPI proc. grid Y coord +C mpiXGlobalLo - My bottom-left (south-west) x-coordinate in +C global domain. +C mpiYGlobalLo - My bottom-left (south-west) y-coordinate in +C global domain. +C mpiTypeXFaceBlock_xy_r4 - Primitives for communicating edge +C mpiTypeXFaceBlock_xy_r8 of a block. +C mpiTypeYFaceBlock_xy_r4 XFace is used in east-west transfer +C mpiTypeYFaceBlock_xy_r8 YFace is used in nrth-south transfer +C mpiTypeXFaceBlock_xyz_r4 xy is used in two-dimensional arrays +C mpiTypeXFaceBlock_xyz_r8 xyz is used with three-dimensional arrays +C mpiTypeYFaceBlock_xyz_r4 r4 is used for real*4 data +C mpiTypeYFaceBlock_xyz_r8 r8 is used for real*8 data +C mpiTypeXFaceThread_xy_r4 - Composites of the above primitives +C mpiTypeXFaceThread_xy_r8 for communicating edges of all blocks +C mpiTypeYFaceThread_xy_r4 owned by a thread. +C mpiTypeYFaceThread_xy_r8 +C mpiTypeXFaceThread_xyz_r4 +C mpiTypeXFaceThread_xyz_r8 +C mpiTypeYFaceThread_xyz_r4 +C mpiTypeYFaceBlock_xyz_r8 +C mpiTagE - Tags are needed to mark requests when MPI is running +C mpiTagW between multithreaded processes or when the same process. +C mpiTagS is a neighbor in more than one direction. The tags ensure that +C mpiTagN a thread will get the message it is looking for. +C mpiTagSW The scheme adopted is to tag messages according to +C mpiTagSE the direction they are travelling. Thus a message +C mpiTagNW travelling east is tagged mpiTagE. However, in a +C mpiTagNE multi-threaded environemnt several messages could +C be travelling east from the same process at the +C same time. The tag is therefore modified to +C be mpiTag[EWS...]*nThreads+myThid. This requires that +C each thread also know the thread ids of its "neighbor" +C threads. + COMMON /EESUPP_MPI_I/ + & mpiPidW, mpiPidE, mpiPidS, mpiPidN, + & mpiPidSE, mpiPidSW, mpiPidNE, mpiPidNW, + & mpiPidIo, mpiMyId, mpiNProcs, mpiComm, + & mpiPx, mpiPy, mpiXGlobalLo, mpiYGlobalLo, + & mpiTypeXFaceBlock_xy_r4, mpiTypeXFaceBlock_xy_r8, + & mpiTypeYFaceBlock_xy_r4, mpiTypeYFaceBlock_xy_r8, + & mpiTypeXFaceBlock_xyz_r4, mpiTypeXFaceBlock_xyz_r8, + & mpiTypeYFaceBlock_xyz_r4, mpiTypeYFaceBlock_xyz_r8, + & mpiTypeXFaceThread_xy_r4, mpiTypeXFaceThread_xy_r8, + & mpiTypeYFaceThread_xy_r4, mpiTypeYFaceThread_xy_r8, + & mpiTypeXFaceThread_xyz_r4, mpiTypeXFaceThread_xyz_r8, + & mpiTypeYFaceThread_xyz_r4, mpiTypeYFaceThread_xyz_r8, + & mpiTagE, mpiTagW, mpiTagN, mpiTagS, + & mpiTagSE, mpiTagSW, mpiTagNW, mpiTagNE + + INTEGER mpiPidW + INTEGER mpiPidE + INTEGER mpiPidS + INTEGER mpiPidN + INTEGER mpiPidSW + INTEGER mpiPidSE + INTEGER mpiPidNW + INTEGER mpiPidNE + INTEGER mpiPidIO + INTEGER mpiMyId + INTEGER mpiNProcs + INTEGER mpiComm + INTEGER mpiPx + INTEGER mpiPy + INTEGER mpiXGlobalLo + INTEGER mpiYGlobalLo + INTEGER mpiTypeXFaceBlock_xy_r4 + INTEGER mpiTypeXFaceBlock_xy_r8 + INTEGER mpiTypeYFaceBlock_xy_r4 + INTEGER mpiTypeYFaceBlock_xy_r8 + INTEGER mpiTypeXFaceBlock_xyz_r4 + INTEGER mpiTypeXFaceBlock_xyz_r8 + INTEGER mpiTypeYFaceBlock_xyz_r4 + INTEGER mpiTypeYFaceBlock_xyz_r8 + INTEGER mpiTypeXFaceThread_xy_r4(MAX_NO_THREADS) + INTEGER mpiTypeXFaceThread_xy_r8(MAX_NO_THREADS) + INTEGER mpiTypeYFaceThread_xy_r4(MAX_NO_THREADS) + INTEGER mpiTypeYFaceThread_xy_r8(MAX_NO_THREADS) + INTEGER mpiTypeXFaceThread_xyz_r4(MAX_NO_THREADS) + INTEGER mpiTypeXFaceThread_xyz_r8(MAX_NO_THREADS) + INTEGER mpiTypeYFaceThread_xyz_r4(MAX_NO_THREADS) + INTEGER mpiTypeYFaceThread_xyz_r8(MAX_NO_THREADS) + INTEGER mpiTagNW + INTEGER mpiTagNE + INTEGER mpiTagSW + INTEGER mpiTagSE + INTEGER mpiTagW + INTEGER mpiTagE + INTEGER mpiTagN + INTEGER mpiTagS + +C-- COMMON /MPI_FULLMAP_I/ holds integer arrays of the full list of MPI process +C mpi_myXGlobalLo :: List of all processors bottom-left X-index in global domain +C mpi_myYGlobalLo :: List of all processors bottom-left Y-index in global domain +C Note: needed for mpi gather/scatter routines & singleCpuIO. + COMMON /MPI_FULLMAP_I/ + & mpi_myXGlobalLo, mpi_myYGlobalLo + INTEGER mpi_myXGlobalLo(nPx*nPy) + INTEGER mpi_myYGlobalLo(nPx*nPy) + +C MPI communicator describing this model realization + COMMON /MPI_COMMS/ + & MPI_COMM_MODEL + INTEGER MPI_COMM_MODEL + +#endif /* ALLOW_USE_MPI */ Index: /issm/trunk-jpl/test/MITgcm/code/cpl_issm.F =================================================================== --- /issm/trunk-jpl/test/MITgcm/code/cpl_issm.F (revision 22532) +++ /issm/trunk-jpl/test/MITgcm/code/cpl_issm.F (revision 22532) @@ -0,0 +1,85 @@ +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" + +#define ALLOW_CPL_ISSM + +CBOP +C !ROUTINE: CPL_ISSM +C !INTERFACE: + SUBROUTINE CPL_ISSM( myTime, myIter, myThid ) + +C !DESCRIPTION: \bv +C *================================================================== +C | SUBROUTINE cpl_issm +C | o Couple MITgcm ocean model with ISSM ice sheet model +C *================================================================== +C \ev + +C !USES: + IMPLICIT NONE +C == Global variables == +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "DYNVARS.h" +#include "GRID.h" +#include "FFIELDS.h" +#include "SHELFICE_OPTIONS.h" +#include "SHELFICE.h" +#ifdef ALLOW_EXF +# include "EXF_OPTIONS.h" +# include "EXF_FIELDS.h" +#endif + + LOGICAL DIFFERENT_MULTIPLE + EXTERNAL DIFFERENT_MULTIPLE + +C !LOCAL VARIABLES: +C mytime - time counter for this thread (seconds) +C myiter - iteration counter for this thread +C mythid - thread number for this instance of the routine. + _RL mytime + INTEGER myiter, mythid +CEOP + +#ifdef ALLOW_CPL_ISSM +#include "EESUPPORT.h" + COMMON /CPL_MPI_ID/ + & mpiMyWid, toissmcomm + INTEGER mpiMyWid, toissmcomm + integer mpistatus(MPI_STATUS_SIZE), mpiRC + real*8 CouplingTime, IceModelTime + + print*,'choubi 1' + + IF( myTime .EQ. startTime ) THEN + print*,'choubi 2' + +C Send deltatimestep + _BEGIN_MASTER( myThid ) + + print*,'choubi 3' + + call MPI_Recv(CouplingTime,1,MPI_DOUBLE,0,10001000,toissmcomm, + & mpistatus,mpiRC) + print*, 'Ocean received CouplingTime: ', CouplingTime + call MPI_Recv(IceModelTime,1,MPI_DOUBLE,0,10001001,toissmcomm, + & mpistatus,mpiRC) + print*, 'Ocean received IceModelTime: ', IceModelTime + call MPI_Send(myTime,1,MPI_DOUBLE,0,10001002,toissmcomm, + & mpistatus,mpiRC) + call MPI_Send(Nx,1,MPI_INT,0,10001003,toissmcomm, + & mpistatus,mpiRC) + call MPI_Send(Ny,1,MPI_INT,0,10001004,toissmcomm, + & mpistatus,mpiRC) + +cdb call MPI_Bcast(dummy1,1,MPI_INT,0,MPI_COMM_MODEL,mpiRC) +cdb print*, 'Ocean : dummy received:',dummy1 + + _END_MASTER( myThid ) + ENDIF + +#endif /* ALLOW_CPL_MPMICE */ + + RETURN + END Index: /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F =================================================================== --- /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F (revision 22532) +++ /issm/trunk-jpl/test/MITgcm/code/do_oceanic_phys.F (revision 22532) @@ -0,0 +1,1196 @@ +C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.149 2017/02/12 20:18:24 gforget Exp $ +C $Name: $ + +#include "PACKAGES_CONFIG.h" +#include "CPP_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif +#ifdef ALLOW_CTRL +# include "CTRL_OPTIONS.h" +#endif +#ifdef ALLOW_SALT_PLUME +# include "SALT_PLUME_OPTIONS.h" +#endif +#ifdef ALLOW_ECCO +# include "ECCO_OPTIONS.h" +#endif + +#ifdef ALLOW_AUTODIFF +# ifdef ALLOW_GGL90 +# include "GGL90_OPTIONS.h" +# endif +# ifdef ALLOW_GMREDI +# include "GMREDI_OPTIONS.h" +# endif +# ifdef ALLOW_KPP +# include "KPP_OPTIONS.h" +# endif +# ifdef ALLOW_SEAICE +# include "SEAICE_OPTIONS.h" +# endif +# ifdef ALLOW_EXF +# include "EXF_OPTIONS.h" +# endif +#endif /* ALLOW_AUTODIFF */ + +CBOP +C !ROUTINE: DO_OCEANIC_PHYS +C !INTERFACE: + SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid) +C !DESCRIPTION: \bv +C *==========================================================* +C | SUBROUTINE DO_OCEANIC_PHYS +C | o Controlling routine for oceanic physics and +C | parameterization +C *==========================================================* +C | o originally, part of S/R thermodynamics +C *==========================================================* +C \ev + +C !CALLING SEQUENCE: +C DO_OCEANIC_PHYS +C | +C |-- OBCS_CALC +C | +C |-- OCN_APPLY_IMPORT +C | +C |-- FRAZIL_CALC_RHS +C | +C |-- THSICE_MAIN +C | +C |-- SEAICE_FAKE +C |-- SEAICE_MODEL +C |-- SEAICE_COST_SENSI +C | +C |-- OCN_EXPORT_DATA +C | +C |-- SHELFICE_THERMODYNAMICS +C | +C |-- ICEFRONT_THERMODYNAMICS +C | +C |-- SALT_PLUME_DO_EXCH +C | +C |-- FREEZE_SURFACE +C | +C |-- EXTERNAL_FORCING_SURF +C | +C |- k loop (Nr:1): +C | - DWNSLP_CALC_RHO +C | - BBL_CALC_RHO +C | - FIND_RHO_2D @ p(k) +C | - FIND_RHO_2D @ p(k-1) +C | - GRAD_SIGMA +C | - CALC_IVDC +C | - DIAGS_RHO_L +C |- end k loop. +C | +C |-- CALC_OCE_MXLAYER +C | +C |-- SALT_PLUME_CALC_DEPTH +C |-- SALT_PLUME_VOLFRAC +C |-- SALT_PLUME_APPLY +C |-- SALT_PLUME_APPLY +C |-- SALT_PLUME_FORCING_SURF +C | +C |-- KPP_CALC +C |-- KPP_CALC_DUMMY +C | +C |-- PP81_CALC +C | +C |-- KL10_CALC +C | +C |-- MY82_CALC +C | +C |-- GGL90_CALC +C | +C |-- TIMEAVE_SURF_FLUX +C | +C |-- GMREDI_CALC_TENSOR +C |-- GMREDI_CALC_TENSOR_DUMMY +C | +C |-- DWNSLP_CALC_FLOW +C |-- DWNSLP_CALC_FLOW +C | +C |-- OFFLINE_GET_DIFFUS +C | +C |-- BBL_CALC_RHS +C | +C |-- MYPACKAGE_CALC_RHS +C | +C |-- GMREDI_DO_EXCH +C | +C |-- KPP_DO_EXCH +C | +C |-- DIAGS_RHO_G +C |-- DIAGS_OCEANIC_SURF_FLUX +C |-- SALT_PLUME_DIAGNOSTICS_FILL +C | +C |-- ECCO_PHYS + +C !USES: + IMPLICIT NONE +C == Global variables === +#include "SIZE.h" +#include "EEPARAMS.h" +#include "PARAMS.h" +#include "GRID.h" +#include "DYNVARS.h" +#ifdef ALLOW_TIMEAVE +# include "TIMEAVE_STATV.h" +#endif +#ifdef ALLOW_OFFLINE +# include "OFFLINE_SWITCH.h" +#endif + +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_MYFIELDS.h" +# include "tamc.h" +# include "tamc_keys.h" +# include "FFIELDS.h" +# include "SURFACE.h" +# include "EOS.h" +# ifdef ALLOW_GMREDI +# include "GMREDI.h" +# endif +# ifdef ALLOW_KPP +# include "KPP.h" +# endif +# ifdef ALLOW_GGL90 +# include "GGL90.h" +# endif +# ifdef ALLOW_EBM +# include "EBM.h" +# endif +# ifdef ALLOW_EXF +# include "ctrl.h" +# include "EXF_FIELDS.h" +# ifdef ALLOW_BULKFORMULAE +# include "EXF_CONSTANTS.h" +# endif +# endif +# ifdef ALLOW_SEAICE +# include "SEAICE_SIZE.h" +# include "SEAICE.h" +# include "SEAICE_PARAMS.h" +# endif +# ifdef ALLOW_THSICE +# include "THSICE_VARS.h" +# endif +# ifdef ALLOW_SALT_PLUME +# include "SALT_PLUME.h" +# endif +# ifdef ALLOW_ECCO +# ifdef ALLOW_SIGMAR_COST_CONTRIBUTION +# include "ecco_cost.h" +# endif +# endif +#endif /* ALLOW_AUTODIFF */ + +C !INPUT/OUTPUT PARAMETERS: +C == Routine arguments == +C myTime :: Current time in simulation +C myIter :: Current iteration number in simulation +C myThid :: Thread number for this instance of the routine. + _RL myTime + INTEGER myIter + INTEGER myThid + +C !LOCAL VARIABLES: +C == Local variables +C rhoK, rhoKm1 :: Density at current level, and level above +C iMin, iMax :: Ranges and sub-block indices on which calculations +C jMin, jMax are applied. +C bi, bj :: tile indices +C msgBuf :: Temp. for building output string +C i,j,k :: loop indices + _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) + _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + INTEGER iMin, iMax + INTEGER jMin, jMax + INTEGER bi, bj + INTEGER i, j, k + CHARACTER*(MAX_LEN_MBUF) msgBuf + INTEGER doDiagsRho + LOGICAL calcGMRedi, calcKPP, calcConvect +#ifdef ALLOW_DIAGNOSTICS + LOGICAL DIAGNOSTICS_IS_ON + EXTERNAL DIAGNOSTICS_IS_ON +#endif /* ALLOW_DIAGNOSTICS */ +#ifdef ALLOW_AUTODIFF + _RL thetaRef +#endif /* ALLOW_AUTODIFF */ +CEOP + +#ifdef ALLOW_AUTODIFF_TAMC +C-- dummy statement to end declaration part + itdkey = 1 +#endif /* ALLOW_AUTODIFF_TAMC */ + +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid) +#endif + + doDiagsRho = 0 +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics .AND. fluidIsWater ) THEN + IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) ) + & doDiagsRho = doDiagsRho + 1 + IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) ) + & doDiagsRho = doDiagsRho + 2 + IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) ) + & doDiagsRho = doDiagsRho + 4 + IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) ) + & doDiagsRho = doDiagsRho + 8 + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + calcGMRedi = useGMRedi + calcKPP = useKPP + calcConvect = ivdc_kappa.NE.0. +#ifdef ALLOW_OFFLINE + IF ( useOffLine ) THEN + calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi + calcKPP = useKPP .AND. .NOT.offlineLoadKPP + calcConvect=calcConvect.AND. .NOT.offlineLoadConvec + ENDIF +#endif /* ALLOW_OFFLINE */ + +#ifdef ALLOW_OBCS + IF (useOBCS) THEN +C-- Calculate future values on open boundaries +C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte +# endif +# ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid) +# endif + CALL OBCS_CALC( myTime+deltaTClock, myIter+1, + I uVel, vVel, wVel, theta, salt, myThid ) + ENDIF +#endif /* ALLOW_OBCS */ + +#ifdef ALLOW_OCN_COMPON_INTERF +C-- Apply imported data (from coupled interface) to forcing fields +C jmc: moved here before any freezing/seaice pkg adjustment of surf-fluxes + IF ( useCoupler ) THEN + CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_OCN_COMPON_INTERF */ + +#ifdef ALLOW_AUTODIFF +# ifdef ALLOW_SALT_PLUME + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + saltPlumeDepth(i,j,bi,bj) = 0. _d 0 + saltPlumeFlux(i,j,bi,bj) = 0. _d 0 + ENDDO + ENDDO + ENDDO + ENDDO +# endif +# ifdef ALLOW_ECCO +# ifdef ALLOW_SIGMAR_COST_CONTRIBUTION + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + sigmaRfield(i,j,k,bi,bj) = 0. _d 0 + ENDDO + ENDDO + ENDDO + ENDDO + ENDDO +# endif +# endif +#endif /* ALLOW_AUTODIFF */ + +#ifdef ALLOW_FRAZIL + IF ( useFRAZIL ) THEN +C-- Freeze water in the ocean interior and let it rise to the surface +CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte + CALL FRAZIL_CALC_RHS( myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_FRAZIL */ + +#ifndef OLD_THSICE_CALL_SEQUENCE +#if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D) + IF ( useThSIce .AND. fluidIsWater ) THEN +# ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE uice,vice = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE snowHeight, Tsrf = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE hocemxl = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE icflxsw = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE salt,theta = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# ifdef NONLIN_FRSURF +CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# endif /* ALLOW_AUTODIFF_TAMC */ +# ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid) +# endif +C-- Step forward Therm.Sea-Ice variables +C and modify forcing terms including effects from ice + CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid) + CALL THSICE_MAIN( myTime, myIter, myThid ) + CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid) + ENDIF +#endif /* ALLOW_THSICE */ +#endif /* ndef OLD_THSICE_CALL_SEQUENCE */ + +#ifdef ALLOW_SEAICE +# ifdef ALLOW_AUTODIFF +CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte +#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP) +CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte +#endif + IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN + CALL SEAICE_FAKE( myTime, myIter, myThid ) + ENDIF +CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE qnet = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE qsw = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte +#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP) +CADJ STORE evap = comlev1, key=ikey_dynamics, kind=isbyte +#endif +# endif /* ALLOW_AUTODIFF */ +#endif /* ALLOW_SEAICE */ + +#ifdef ALLOW_SEAICE + IF ( useSEAICE ) THEN +# ifdef ALLOW_AUTODIFF_TAMC +cph-adj-test( +CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE empmr, qnet = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte +CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte +cCADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte +cCADJ STORE salt = comlev1, key=ikey_dynamics, kind=isbyte +cph-adj-test) +c#ifdef ALLOW_EXF +CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE evap = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE uwind,vwind = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +c#endif +CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# ifdef SEAICE_CGRID +CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# ifdef SEAICE_ALLOW_DYNAMICS +CADJ STORE uice = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE vice = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE dwatn = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# ifdef SEAICE_ALLOW_EVP +CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# endif +# ifdef SEAICE_VARIABLE_SALINITY +CADJ STORE hsalt = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# ifdef ATMOSPHERIC_LOADING +CADJ STORE pload = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE siceload = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# ifdef NONLIN_FRSURF +CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# ifdef ANNUAL_BALANCE +CADJ STORE balance_itcount = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif /* ANNUAL_BALANCE */ +# ifdef ALLOW_THSICE +C-- store thSIce vars before advection (called from SEAICE_MODEL) update them: +CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE Qice1, Qice2 = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif /* ALLOW_THSICE */ +# endif /* ALLOW_AUTODIFF_TAMC */ +# ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid) +# endif + CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid) + CALL SEAICE_MODEL( myTime, myIter, myThid ) + CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid) +# ifdef ALLOW_COST + CALL SEAICE_COST_SENSI ( myTime, myIter, myThid ) +# endif + ENDIF +#endif /* ALLOW_SEAICE */ + +#if (defined ALLOW_OCN_COMPON_INTERF) && (defined ALLOW_THSICE) +C-- After seaice-dyn and advection of pkg/thsice fields, +C Export ocean coupling fields to coupled interface (only with pkg/thsice) + IF ( useCoupler ) THEN +# ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('OCN_EXPORT_DATA',myThid) +# endif + CALL TIMER_START('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid) + CALL OCN_EXPORT_DATA( myTime, myIter, myThid ) + CALL TIMER_STOP ('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid) + ENDIF +#endif /* ALLOW_OCN_COMPON_INTERF & ALLOW_THSICE */ + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE sst, sss = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE qsw = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# ifdef ALLOW_SEAICE +CADJ STORE area = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +#endif + +#ifdef OLD_THSICE_CALL_SEQUENCE +#if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D) + IF ( useThSIce .AND. fluidIsWater ) THEN +# ifdef ALLOW_AUTODIFF_TAMC +cph( +# ifdef NONLIN_FRSURF +CADJ STORE uice,vice = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE salt,theta = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE qnet,qsw, empmr = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE hFac_surfC = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +# endif +# endif +# ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid) +# endif +C-- Step forward Therm.Sea-Ice variables +C and modify forcing terms including effects from ice + CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid) + CALL THSICE_MAIN( myTime, myIter, myThid ) + CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid) + ENDIF +#endif /* ALLOW_THSICE */ +#endif /* OLD_THSICE_CALL_SEQUENCE */ + +#define ALLOW_CPL_ISSM +#ifdef ALLOW_CPL_ISSM + CALL CPL_ISSM( myTime, myIter, myThid ) +#endif + +#ifdef ALLOW_SHELFICE + IF ( useShelfIce .AND. fluidIsWater ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid) +#endif +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE salt, theta = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +#endif +C compute temperature and (virtual) salt flux at the +C shelf-ice ocean interface + CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]', + & myThid) + CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid ) + CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]', + & myThid) + ENDIF +#endif /* ALLOW_SHELFICE */ + +#ifdef ALLOW_ICEFRONT + IF ( useICEFRONT .AND. fluidIsWater ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid) +#endif +C compute temperature and (virtual) salt flux at the +C ice-front ocean interface + CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]', + & myThid) + CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid ) + CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]', + & myThid) + ENDIF +#endif /* ALLOW_ICEFRONT */ + +#ifdef ALLOW_SALT_PLUME + IF ( useSALT_PLUME ) THEN +Catn: exchanging saltPlumeFlux: + CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_SALT_PLUME */ + +C-- Freeze water at the surface + IF ( allowFreezing ) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE theta = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +#endif + CALL FREEZE_SURFACE( myTime, myIter, myThid ) + ENDIF + + iMin = 1-OLx + iMax = sNx+OLx + jMin = 1-OLy + jMax = sNy+OLy + +C--- Determines forcing terms based on external fields +C relaxation terms, etc. +#ifdef ALLOW_AUTODIFF +CADJ STORE salt, theta = comlev1, key = ikey_dynamics, +CADJ & kind = isbyte +#else /* ALLOW_AUTODIFF */ +C-- if fluid is not water, by-pass surfaceForcing, find_rho, gmredi +C and all vertical mixing schemes, but keep OBCS_CALC + IF ( fluidIsWater ) THEN +#endif /* ALLOW_AUTODIFF */ +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid) +#endif + CALL EXTERNAL_FORCING_SURF( + I iMin, iMax, jMin, jMax, + I myTime, myIter, myThid ) + +#ifdef ALLOW_AUTODIFF_TAMC +C-- HPF directive to help TAMC +CHPF$ INDEPENDENT +#endif /* ALLOW_AUTODIFF_TAMC */ + DO bj=myByLo(myThid),myByHi(myThid) +#ifdef ALLOW_AUTODIFF_TAMC +C-- HPF directive to help TAMC +CHPF$ INDEPENDENT +#endif /* ALLOW_AUTODIFF_TAMC */ + DO bi=myBxLo(myThid),myBxHi(myThid) + +#ifdef ALLOW_AUTODIFF_TAMC + act1 = bi - myBxLo(myThid) + max1 = myBxHi(myThid) - myBxLo(myThid) + 1 + act2 = bj - myByLo(myThid) + max2 = myByHi(myThid) - myByLo(myThid) + 1 + act3 = myThid - 1 + max3 = nTx*nTy + act4 = ikey_dynamics - 1 + itdkey = (act1 + 1) + act2*max1 + & + act3*max1*max2 + & + act4*max1*max2*max3 +#endif /* ALLOW_AUTODIFF_TAMC */ + +C-- Set up work arrays with valid (i.e. not NaN) values +C These inital values do not alter the numerical results. They +C just ensure that all memory references are to valid floating +C point numbers. This prevents spurious hardware signals due to +C uninitialised but inert locations. + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx +C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics + sigmaX(i,j,k) = 0. _d 0 + sigmaY(i,j,k) = 0. _d 0 + sigmaR(i,j,k) = 0. _d 0 + ENDDO + ENDDO + ENDDO + +#ifdef ALLOW_AUTODIFF + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + rhoKm1 (i,j) = 0. _d 0 + rhoKp1 (i,j) = 0. _d 0 + ENDDO + ENDDO +cph all the following init. are necessary for TAF +cph although some of these are re-initialised later. + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + rhoInSitu(i,j,k,bi,bj) = 0. +# ifdef ALLOW_GGL90 + GGL90viscArU(i,j,k,bi,bj) = 0. _d 0 + GGL90viscArV(i,j,k,bi,bj) = 0. _d 0 + GGL90diffKr(i,j,k,bi,bj) = 0. _d 0 +# endif /* ALLOW_GGL90 */ +# ifdef ALLOW_SALT_PLUME +# ifdef SALT_PLUME_VOLUME + SPforcingS(i,j,k,bi,bj) = 0. _d 0 + SPforcingT(i,j,k,bi,bj) = 0. _d 0 +# endif +# endif /* ALLOW_SALT_PLUME */ + ENDDO + ENDDO + ENDDO +#ifdef ALLOW_OFFLINE + IF ( calcConvect ) THEN +#endif + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + IVDConvCount(i,j,k,bi,bj) = 0. + ENDDO + ENDDO + ENDDO +#ifdef ALLOW_OFFLINE + ENDIF + IF ( calcGMRedi ) THEN +#endif +# ifdef ALLOW_GMREDI + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + Kwx(i,j,k,bi,bj) = 0. _d 0 + Kwy(i,j,k,bi,bj) = 0. _d 0 + Kwz(i,j,k,bi,bj) = 0. _d 0 +# ifdef GM_NON_UNITY_DIAGONAL + Kux(i,j,k,bi,bj) = 0. _d 0 + Kvy(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_EXTRA_DIAGONAL + Kuz(i,j,k,bi,bj) = 0. _d 0 + Kvz(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_BOLUS_ADVEC + GM_PsiX(i,j,k,bi,bj) = 0. _d 0 + GM_PsiY(i,j,k,bi,bj) = 0. _d 0 +# endif +# ifdef GM_VISBECK_VARIABLE_K + VisbeckK(i,j,bi,bj) = 0. _d 0 +# endif + ENDDO + ENDDO + ENDDO +# endif /* ALLOW_GMREDI */ +#ifdef ALLOW_OFFLINE + ENDIF + IF ( calcKPP ) THEN +#endif +# ifdef ALLOW_KPP + DO k=1,Nr + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0 + KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0 + ENDDO + ENDDO + ENDDO +# endif /* ALLOW_KPP */ +#ifdef ALLOW_OFFLINE + ENDIF +#endif +#endif /* ALLOW_AUTODIFF */ + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE totphihyd(:,:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +# ifdef ALLOW_KPP +CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +# endif +# ifdef ALLOW_SALT_PLUME +CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +C-- Always compute density (stored in common block) here; even when it is not +C needed here, will be used anyway in calc_phi_hyd (data flow easier this way) +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid) +#endif +#ifdef ALLOW_AUTODIFF + IF ( fluidIsWater ) THEN +#endif /* ALLOW_AUTODIFF */ +#ifdef ALLOW_DOWN_SLOPE + IF ( useDOWN_SLOPE ) THEN + DO k=1,Nr + CALL DWNSLP_CALC_RHO( + I theta, salt, + O rhoInSitu(1-OLx,1-OLy,k,bi,bj), + I k, bi, bj, myTime, myIter, myThid ) + ENDDO + ENDIF +#endif /* ALLOW_DOWN_SLOPE */ +#ifdef ALLOW_BBL + IF ( useBBL ) THEN +C pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl. +C To reduce computation and storage requirement, these densities are stored in the +C dry grid boxes of rhoInSitu. See BBL_CALC_RHO for details. + DO k=Nr,1,-1 + CALL BBL_CALC_RHO( + I theta, salt, + O rhoInSitu, + I k, bi, bj, myTime, myIter, myThid ) + + ENDDO + ENDIF +#endif /* ALLOW_BBL */ + IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN + DO k=1,Nr + CALL FIND_RHO_2D( + I iMin, iMax, jMin, jMax, k, + I theta(1-OLx,1-OLy,k,bi,bj), + I salt (1-OLx,1-OLy,k,bi,bj), + O rhoInSitu(1-OLx,1-OLy,k,bi,bj), + I k, bi, bj, myThid ) + ENDDO + ENDIF +#ifdef ALLOW_AUTODIFF + ELSE +C- fluid is not water: + DO k=1,Nr + IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN +C- isothermal (theta=const) reference state + thetaRef = thetaConst + ELSE +C- horizontally uniform (tRef) reference state + thetaRef = tRef(k) + ENDIF + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + rhoInSitu(i,j,k,bi,bj) = + & ( theta(i,j,k,bi,bj) + & *( salt(i,j,k,bi,bj)*atm_Rq + oneRL ) + & - thetaRef )*maskC(i,j,k,bi,bj) + ENDDO + ENDDO + ENDDO + ENDIF +#endif /* ALLOW_AUTODIFF */ + +#ifdef ALLOW_DEBUG + IF (debugMode) THEN + WRITE(msgBuf,'(A,2(I4,A))') + & 'ENTERING UPWARD K LOOP (bi=', bi, ', bj=', bj,')' + CALL DEBUG_MSG(msgBuf(1:43),myThid) + ENDIF +#endif + +C-- Start of diagnostic loop + DO k=Nr,1,-1 + +#ifdef ALLOW_AUTODIFF_TAMC +C? Patrick, is this formula correct now that we change the loop range? +C? Do we still need this? +cph kkey formula corrected. +cph Needed for rhoK, rhoKm1, in the case useGMREDI. + kkey = (itdkey-1)*Nr + k +#endif /* ALLOW_AUTODIFF_TAMC */ + +c#ifdef ALLOW_AUTODIFF_TAMC +cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, +cCADJ & kind = isbyte +cCADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, +cCADJ & kind = isbyte +c#endif /* ALLOW_AUTODIFF_TAMC */ + +C-- Calculate gradients of potential density for isoneutral +C slope terms (e.g. GM/Redi tensor or IVDC diffusivity) + IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect) + & .OR. usePP81 .OR. useKL10 + & .OR. useMY82 .OR. useGGL90 + & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN + IF (k.GT.1) THEN +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, +CADJ & kind = isbyte +CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, +CADJ & kind = isbyte +CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey, +CADJ & kind = isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + CALL FIND_RHO_2D( + I iMin, iMax, jMin, jMax, k, + I theta(1-OLx,1-OLy,k-1,bi,bj), + I salt (1-OLx,1-OLy,k-1,bi,bj), + O rhoKm1, + I k-1, bi, bj, myThid ) + ENDIF +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid) +#endif +cph Avoid variable aliasing for adjoint !!! + DO j=jMin,jMax + DO i=iMin,iMax + rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj) + ENDDO + ENDDO + CALL GRAD_SIGMA( + I bi, bj, iMin, iMax, jMin, jMax, k, + I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1, + O sigmaX, sigmaY, sigmaR, + I myThid ) +#ifdef ALLOW_ECCO +# ifdef ALLOW_SIGMAR_COST_CONTRIBUTION + DO j=jMin,jMax + DO i=iMin,iMax + sigmaRfield(i,j,k,bi,bj)=sigmaR(i,j,k) + ENDDO + ENDDO +# endif +#endif /* ALLOW_ECCO */ +#ifdef ALLOW_AUTODIFF +#ifdef GMREDI_WITH_STABLE_ADJOINT +cgf zero out adjoint fields to stabilize pkg/gmredi adjoint +cgf -> cuts adjoint dependency from slope to state + CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid) + CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid) + CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid) +#endif +#endif /* ALLOW_AUTODIFF */ + ENDIF + +C-- Implicit Vertical Diffusion for Convection + IF (k.GT.1 .AND. calcConvect) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid) +#endif + CALL CALC_IVDC( + I bi, bj, iMin, iMax, jMin, jMax, k, + I sigmaR, + I myTime, myIter, myThid) + ENDIF + +#ifdef ALLOW_DIAGNOSTICS + IF ( doDiagsRho.GE.4 ) THEN + CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj, + I rhoInSitu(1-OLx,1-OLy,1,bi,bj), + I rhoKm1, wVel, + I myTime, myIter, myThid ) + ENDIF +#endif + +C-- end of diagnostic k loop (Nr:1) + ENDDO + +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE IVDConvCount(:,:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +#endif + +C-- Diagnose Mixed Layer Depth: + IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN + CALL CALC_OCE_MXLAYER( + I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR, + I bi, bj, myTime, myIter, myThid ) + ENDIF + +#ifdef ALLOW_SALT_PLUME + IF ( useSALT_PLUME ) THEN + CALL SALT_PLUME_CALC_DEPTH( + I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR, + I bi, bj, myTime, myIter, myThid ) +#ifdef SALT_PLUME_VOLUME + CALL SALT_PLUME_VOLFRAC( + I bi, bj, myTime, myIter, myThid ) +C-- get forcings for kpp + CALL SALT_PLUME_APPLY( + I 1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj), + I theta, 0, + I myTime, myIter, myThid ) + CALL SALT_PLUME_APPLY( + I 2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj), + I salt, 0, + I myTime, myIter, myThid ) +C-- need to call this S/R from here to apply just before kpp + CALL SALT_PLUME_FORCING_SURF( + I bi, bj, iMin, iMax, jMin, jMax, + I myTime, myIter, myThid ) +#endif /* SALT_PLUME_VOLUME */ + ENDIF +#endif /* ALLOW_SALT_PLUME */ + +#ifdef ALLOW_DIAGNOSTICS + IF ( MOD(doDiagsRho,4).GE.2 ) THEN + CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr, + & 2, bi, bj, myThid) + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + +C-- This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called; +C now called earlier, before bi,bj loop. + +#ifdef ALLOW_AUTODIFF_TAMC +cph needed for KPP +CADJ STORE surfaceForcingU(:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE surfaceForcingV(:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE surfaceForcingS(:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE surfaceForcingT(:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE surfaceForcingTice(:,:,bi,bj) +CADJ & = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ + +#ifdef ALLOW_KPP +C-- Compute KPP mixing coefficients + IF ( calcKPP ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid) +#endif + CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid) + CALL KPP_CALC( + I bi, bj, myTime, myIter, myThid ) + CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid) +#if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE) + ELSE + CALL KPP_CALC_DUMMY( + I bi, bj, myTime, myIter, myThid ) +#endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */ + ENDIF +#endif /* ALLOW_KPP */ + +#ifdef ALLOW_PP81 +C-- Compute PP81 mixing coefficients + IF (usePP81) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid) +#endif + CALL PP81_CALC( + I bi, bj, sigmaR, myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_PP81 */ + +#ifdef ALLOW_KL10 +C-- Compute KL10 mixing coefficients + IF (useKL10) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('KL10_CALC',myThid) +#endif + CALL KL10_CALC( + I bi, bj, sigmaR, myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_KL10 */ + +#ifdef ALLOW_MY82 +C-- Compute MY82 mixing coefficients + IF (useMY82) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid) +#endif + CALL MY82_CALC( + I bi, bj, sigmaR, myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_MY82 */ + +#ifdef ALLOW_GGL90 +#ifdef ALLOW_AUTODIFF_TAMC +CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +#endif /* ALLOW_AUTODIFF_TAMC */ +C-- Compute GGL90 mixing coefficients + IF (useGGL90) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid) +#endif + CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid) + CALL GGL90_CALC( + I bi, bj, sigmaR, myTime, myIter, myThid ) + CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid) + ENDIF +#endif /* ALLOW_GGL90 */ + +#ifdef ALLOW_TIMEAVE + IF ( taveFreq.GT. 0. _d 0 ) THEN + CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid) + ENDIF + IF ( taveFreq.GT.0. .AND. calcConvect ) THEN + CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount, + I Nr, deltaTClock, bi, bj, myThid) + ENDIF +#endif /* ALLOW_TIMEAVE */ + +#ifdef ALLOW_GMREDI +#ifdef ALLOW_AUTODIFF_TAMC +# ifndef GM_EXCLUDE_CLIPPING +cph storing here is needed only for one GMREDI_OPTIONS: +cph define GM_BOLUS_ADVEC +cph keep it although TAF says you dont need to. +cph but I have avoided the #ifdef for now, in case more things change +CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, +CADJ & kind = isbyte +# endif +#endif /* ALLOW_AUTODIFF_TAMC */ + +C-- Calculate iso-neutral slopes for the GM/Redi parameterisation + IF ( calcGMRedi ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid) +#endif + CALL GMREDI_CALC_TENSOR( + I iMin, iMax, jMin, jMax, + I sigmaX, sigmaY, sigmaR, + I bi, bj, myTime, myIter, myThid ) +#if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE) + ELSE + CALL GMREDI_CALC_TENSOR_DUMMY( + I iMin, iMax, jMin, jMax, + I sigmaX, sigmaY, sigmaR, + I bi, bj, myTime, myIter, myThid ) +#endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */ + ENDIF +#endif /* ALLOW_GMREDI */ + +#ifdef ALLOW_DOWN_SLOPE + IF ( useDOWN_SLOPE ) THEN +C-- Calculate Downsloping Flow for Down_Slope parameterization + IF ( usingPCoords ) THEN + CALL DWNSLP_CALC_FLOW( + I bi, bj, kSurfC, rhoInSitu, + I myTime, myIter, myThid ) + ELSE + CALL DWNSLP_CALC_FLOW( + I bi, bj, kLowC, rhoInSitu, + I myTime, myIter, myThid ) + ENDIF + ENDIF +#endif /* ALLOW_DOWN_SLOPE */ + +C-- end bi,bj loops. + ENDDO + ENDDO + +#ifndef ALLOW_AUTODIFF +C--- if fluid Is Water: end + ENDIF +#endif + +#ifdef ALLOW_OFFLINE + IF ( useOffLine ) THEN +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_CALL('OFFLINE_GET_DIFFUS',myThid) +#endif /* ALLOW_DEBUG */ + CALL OFFLINE_GET_DIFFUS( myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_OFFLINE */ + +#ifdef ALLOW_BBL + IF ( useBBL ) THEN + CALL BBL_CALC_RHS( + I myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_BBL */ + +#ifdef ALLOW_MYPACKAGE + IF ( useMYPACKAGE ) THEN + CALL MYPACKAGE_CALC_RHS( + I myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_MYPACKAGE */ + +#ifdef ALLOW_GMREDI + IF ( calcGMRedi ) THEN + CALL GMREDI_DO_EXCH( myTime, myIter, myThid ) + ENDIF +#endif /* ALLOW_GMREDI */ + +#ifdef ALLOW_KPP + IF ( calcKPP ) THEN + CALL KPP_DO_EXCH( myThid ) + ENDIF +#endif /* ALLOW_KPP */ + +#ifdef ALLOW_DIAGNOSTICS + IF ( fluidIsWater .AND. useDiagnostics ) THEN + CALL DIAGS_RHO_G( + I rhoInSitu, uVel, vVel, wVel, + I myTime, myIter, myThid ) + ENDIF + IF ( useDiagnostics ) THEN + CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid ) + ENDIF + IF ( calcConvect .AND. useDiagnostics ) THEN + CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ', + & 0, Nr, 0, 1, 1, myThid ) + ENDIF +#ifdef ALLOW_SALT_PLUME + IF ( useDiagnostics ) + & CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid) +#endif +#endif + +#ifdef ALLOW_DEBUG + IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid) +#endif + + RETURN + END Index: /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F =================================================================== --- /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F (revision 22531) +++ /issm/trunk-jpl/test/MITgcm/code/eeboot_minimal.F (revision 22532) @@ -67,10 +67,9 @@ #ifdef ALLOW_CPL_ISSM COMMON /CPL_MPI_ID/ - & mpiMyWid - integer :: mpiMyWid, numprocsworld - integer my_local_rank,my_local_size - integer toissmcomm - integer dummy1(1),dummy2(1) - integer, dimension(mpi_status_size) :: status + & mpiMyWid, toissmcomm + INTEGER mpiMyWid, toissmcomm + integer my_local_rank,my_local_size, numprocsworld + integer dummy1(1),dummy2(1) + integer status(MPI_STATUS_SIZE) #endif /* ALLOW_CPL_ISSM */ #if defined(ALLOW_NEST_PARENT) || defined(ALLOW_NEST_CHILD) @@ -210,11 +209,10 @@ if (my_local_rank .eq. 0) then - call MPI_Send(dummy2,1,MPI_INT,0,1,toissmcomm,mpiRC) - call MPI_Recv(dummy1,1,MPI_INT,0,1,toissmcomm,status,mpiRC) +cdb call MPI_Send(dummy2,1,MPI_INT,0,1,toissmcomm,mpiRC) +cdb call MPI_Recv(dummy1,1,MPI_INT,0,1,toissmcomm,status,mpiRC) end if - call MPI_Bcast(dummy1,1,MPI_INT,0,MPI_COMM_MODEL,mpiRC) - print*, 'Ocean : dummy received:',dummy1 - +cdb call MPI_Bcast(dummy1,1,MPI_INT,0,MPI_COMM_MODEL,mpiRC) +cdb print*, 'Ocean : dummy received:',dummy1 #endif /* ALLOW_CPL_ISSM */