C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_fluxcalc.F,v 1.17 2015/06/15 21:34:07 jmc Exp $ C $Name: $ #include "LAYERS_OPTIONS.h" #ifdef ALLOW_GMREDI #include "GMREDI_OPTIONS.h" #endif C-- File layers_fluxcalc.F: C-- Contents C-- o LAYERS_FLUXCALC C-- o LAYERS_DIAPYCNAL C-- o LAYERS_LOCATE CBOP 0 C !ROUTINE: LAYERS_FLUXCALC C !INTERFACE: SUBROUTINE LAYERS_FLUXCALC( I uVel,vVel,tracer,iLa, O UH,VH,Hw,Hs,PIw,PIs,Uw,Vs, I myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE LAYERS_FLUXCALC C | Calculate the transport in isotracer layers, for a chosen C | tracer. This is the meat of the LAYERS package. C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "LAYERS_SIZE.h" #include "LAYERS.h" #ifdef ALLOW_GMREDI # include "GMREDI.h" #endif C !INPUT PARAMETERS: C myThid :: my Thread Id number C uVel :: zonal velocity (m/s, i=1 held at western face) C vVel :: meridional velocity (m/s, j=1 held at southern face) C tracer :: potential temperature, salt or potential density prho C UH :: U integrated over layer (m^2/s) C VH :: V integrated over layer (m^2/s) C Hw :: Layer thickness at the U point (m) C Hs :: Layer thickness at the V point (m) C PIw :: 1 if layer exists, 0 otherwise (at U point) C PIs :: 1 if layer exists, 0 otherwise (at V point) C Uw :: average U over layer (m/s) C Vs :: average V over layer (m/s) C iLa :: layer coordinate index INTEGER myThid _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy) _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy) _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy) #ifdef LAYERS_UFLUX _RL UH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) # ifdef LAYERS_THICKNESS _RL Hw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL PIw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL Uw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) else _RL Hw(1), PIw(1), Uw(1) # endif #else _RL UH(1), Hw(1), PIw(1), Uw(1) #endif #ifdef LAYERS_VFLUX _RL VH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) # ifdef LAYERS_THICKNESS _RL Hs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL PIs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL Vs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) else _RL Hs(1), PIs(1), Vs(1) # endif #else _RL VH(1), Hs(1), PIs(1), Vs(1) #endif INTEGER iLa CEOP #ifdef ALLOW_LAYERS C !LOCAL VARIABLES: C bi, bj :: tile indices C i,j :: horizontal indices C k :: vertical index for model grid C kci :: index from CellIndex C kg :: index for looping though layers_bounds C kk :: vertical index for ZZ (fine) grid C kgu,kgv :: vertical index for isopycnal grid C kloc :: local copy of kgu/v to reduce accesses to index arrays C mSteps :: maximum number of steps for bisection method C prho :: pot. density (less 1000) referenced to layers_krho pressure C TatU :: temperature at U point C TatV :: temperature at V point C dzfac :: temporary sublayer thickness C Tloc,Tp1 :: horizontally interpolated tracer values INTEGER bi, bj INTEGER i,j,k,kk,kg,kci,kp1,kloc INTEGER mSteps INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1) _RL TatU(sNx+1,sNy+1), TatV(sNx+1,sNy+1) _RL dzfac #ifdef ALLOW_GMREDI INTEGER kcip1 _RL delPsi, maskp1 #endif LOGICAL errorFlag CHARACTER*(MAX_LEN_MBUF) msgBuf C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C compute maximum number of steps for bisection method (approx. C log2(Nlayers)) as log2(Nlayers) + 1 for safety mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1 C --- The tile loops DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C Initialize the search indices DO j = 1,sNy+1 DO i = 1,sNx+1 C The temperature index (layer_G) goes from cold to warm. C The water column goes from warm (k=1) to cold (k=Nr). C So initialize the search with the warmest value. kgu(i,j) = Nlayers kgv(i,j) = Nlayers ENDDO ENDDO C Reset the arrays DO kg=1,Nlayers DO j = 1-OLy,sNy+OLy DO i = 1-OLx,sNx+OLx #ifdef LAYERS_UFLUX UH (i,j,kg,bi,bj) = 0. _d 0 #ifdef LAYERS_THICKNESS Hw(i,j,kg,bi,bj) = 0. _d 0 PIw(i,j,kg,bi,bj) = 0. _d 0 Uw(i,j,kg,bi,bj) = 0. _d 0 #endif /* LAYERS_THICKNESS */ #endif /* UH */ #ifdef LAYERS_VFLUX VH (i,j,kg,bi,bj) = 0. _d 0 #ifdef LAYERS_THICKNESS Hs(i,j,kg,bi,bj) = 0. _d 0 PIs(i,j,kg,bi,bj) = 0. _d 0 Vs(i,j,kg,bi,bj) = 0. _d 0 #endif /* LAYERS_THICKNESS */ #endif /* VH */ ENDDO ENDDO ENDDO DO kk=1,NZZ k = MapIndex(kk) kci = CellIndex(kk) #ifdef ALLOW_GMREDI kcip1 = MIN(kci+1,Nr) maskp1 = 1. IF (kci.GE.Nr) maskp1 = 0. #endif /* ALLOW_GMREDI */ #ifdef LAYERS_UFLUX DO j = 1,sNy+1 DO i = 1,sNx+1 C ------ Find theta at the U point (west) on the fine Z grid kp1=k+1 IF (maskW(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k TatU(i,j) = MapFact(kk) * & 0.5 _d 0 * (tracer(i-1,j,k,bi,bj)+tracer(i,j,k,bi,bj)) + & (1. _d 0 -MapFact(kk)) * & 0.5 _d 0 * (tracer(i-1,j,kp1,bi,bj)+tracer(i,j,kp1,bi,bj)) ENDDO ENDDO C ------ Now that we know T everywhere, determine the binning. C find the layer indices kgu CALL LAYERS_LOCATE( I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatU, O kgu, I myThid ) #ifndef TARGET_NEC_SX C check for failures IF ( debugLevel .GE. debLevC ) THEN errorFlag = .FALSE. DO j = 1,sNy+1 DO i = 1,sNx+1 IF ( kgu(i,j) .LE. 0 ) THEN WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)') & 'S/R LAYERS_LOCATE: Could not find a bin in ', & 'layers_bounds for TatU(',i,',',j,',)=',TatU(i,j) CALL PRINT_ERROR( msgBuf, myThid ) errorFlag = .TRUE. ENDIF ENDDO ENDDO IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC' ENDIF #endif /* ndef TARGET_NEC_SX */ C DO j = 1,sNy+1 DO i = 1,sNx+1 kloc = kgu(i,j) dzfac = dZZf(kk) * hFacW(i,j,kci,bi,bj) C ------ Augment the bin values UH(i,j,kloc,bi,bj) = & UH(i,j,kloc,bi,bj) + & dzfac * uVel(i,j,kci,bi,bj) #ifdef ALLOW_GMREDI IF ( layers_bolus(iLa) ) THEN IF ( .NOT.GM_AdvForm ) THEN delPsi = 0.25 _d 0 *( & ( rA(i-1,j,bi,bj)*Kwx(i-1,j,kcip1,bi,bj) & +rA( i ,j,bi,bj)*Kwx( i ,j,kcip1,bi,bj) & ) * maskW(i,j,kcip1,bi,bj) * maskp1 & - ( rA(i-1,j,bi,bj)*Kwx(i-1,j, kci ,bi,bj) & +rA( i ,j,bi,bj)*Kwx( i ,j, kci ,bi,bj) & ) * maskW(i,j, kci ,bi,bj) & ) * recip_rAw(i,j,bi,bj) #ifdef GM_BOLUS_ADVEC ELSE delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1 & - GM_PsiX(i,j, kci, bi,bj) #endif ENDIF UH(i,j,kloc,bi,bj) = UH(i,j,kloc,bi,bj) & + delPsi*recip_drF(kci)*_recip_hFacW(i,j,kci,bi,bj) & * dzfac ENDIF #endif /* ALLOW_GMREDI */ #ifdef LAYERS_THICKNESS Hw(i,j,kloc,bi,bj) = Hw(i,j,kloc,bi,bj) + dzfac #endif /* LAYERS_THICKNESS */ ENDDO ENDDO #endif /* LAYERS_UFLUX */ #ifdef LAYERS_VFLUX DO j = 1,sNy+1 DO i = 1,sNx+1 C ------ Find theta at the V point (south) on the fine Z grid kp1=k+1 IF (maskS(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k TatV(i,j) = MapFact(kk) * & 0.5 _d 0 * (tracer(i,j-1,k,bi,bj)+tracer(i,j,k,bi,bj)) + & (1. _d 0 -MapFact(kk)) * & 0.5 _d 0 * (tracer(i,j-1,kp1,bi,bj)+tracer(i,j,kp1,bi,bj)) ENDDO ENDDO C ------ Now that we know T everywhere, determine the binning. C find the layer indices kgv CALL LAYERS_LOCATE( I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatV, O kgv, I myThid ) #ifndef TARGET_NEC_SX IF ( debugLevel .GE. debLevC ) THEN C check for failures errorFlag = .FALSE. DO j = 1,sNy+1 DO i = 1,sNx+1 IF ( kgv(i,j) .LE. 0 ) THEN WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)') & 'S/R LAYERS_LOCATE: Could not find a bin in ', & 'layers_bounds for TatV(',i,',',j,',)=',TatV(i,j) CALL PRINT_ERROR( msgBuf, myThid ) errorFlag = .TRUE. ENDIF ENDDO ENDDO IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC' ENDIF #endif /* ndef TARGET_NEC_SX */ C DO j = 1,sNy+1 DO i = 1,sNx+1 kloc = kgv(i,j) dzfac = dZZf(kk) * hFacS(i,j,kci,bi,bj) C ------ debugging stuff C IF (i.EQ.10 .AND. j.EQ.10) THEN C WRITE(msgBuf,'(A,I3,A,F10.2,A,F6.2)') C & ' kloc=', kloc, C & ', TatV=',TatV(i,j), C & ', dzfac=',dzfac C CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, C & SQUEEZE_RIGHT, myThid ) C ENDIF C ------ Augment the bin values VH(i,j,kloc,bi,bj) = & VH(i,j,kloc,bi,bj) + dzfac * vVel(i,j,kci,bi,bj) #ifdef ALLOW_GMREDI IF ( layers_bolus(iLa) ) THEN IF ( .NOT.GM_AdvForm ) THEN delPsi = 0.25 _d 0 *( & ( rA(i,j-1,bi,bj)*Kwy(i,j-1,kcip1,bi,bj) & +rA(i, j ,bi,bj)*Kwy(i, j ,kcip1,bi,bj) & ) * maskS(i,j,kcip1,bi,bj) * maskp1 & - ( rA(i,j-1,bi,bj)*Kwy(i,j-1, kci ,bi,bj) & +rA(i, j ,bi,bj)*Kwy(i, j , kci ,bi,bj) & ) * maskS(i,j, kci ,bi,bj) & ) * recip_rAs(i,j,bi,bj) #ifdef GM_BOLUS_ADVEC ELSE delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1 & - GM_PsiY(i,j, kci, bi,bj) #endif ENDIF VH(i,j,kloc,bi,bj) = VH(i,j,kloc,bi,bj) & + delPsi*recip_drF(kci)*_recip_hFacS(i,j,kci,bi,bj) & * dzfac ENDIF #endif /* ALLOW_GMREDI */ #ifdef LAYERS_THICKNESS Hs(i,j,kloc,bi,bj) = Hs(i,j,kloc,bi,bj) + dzfac #endif /* LAYERS_THICKNESS */ ENDDO ENDDO #endif /* LAYERS_VFLUX */ ENDDO C-- Now that we know the thicknesses, compute the heaviside function C-- (Needs another loop through Ng) #ifdef LAYERS_THICKNESS DO kg=1,Nlayers DO j = 1,sNy+1 DO i = 1,sNx+1 #ifdef LAYERS_UFLUX IF (Hw(i,j,kg,bi,bj) .GT. 0.) THEN PIw(i,j,kg,bi,bj) = 1. _d 0 Uw(i,j,kg,bi,bj) = & UH(i,j,kg,bi,bj) / Hw(i,j,kg,bi,bj) ENDIF #endif /* LAYERS_UFLUX */ #ifdef LAYERS_VFLUX IF (Hs(i,j,kg,bi,bj) .GT. 0.) THEN PIs(i,j,kg,bi,bj) = 1. _d 0 Vs(i,j,kg,bi,bj) = & VH(i,j,kg,bi,bj) / Hs(i,j,kg,bi,bj) ENDIF #endif /* LAYERS_VFLUX */ ENDDO ENDDO ENDDO #endif /* LAYERS_THICKNESS */ C --- End bi,bj loop ENDDO ENDDO RETURN END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP 0 C !ROUTINE: LAYERS_DIAPYCNAL C !INTERFACE: SUBROUTINE LAYERS_DIAPYCNAL( I tracer,iLa, O TtendSurf, TtendDiffh, TtendDiffr, O TtendAdvh, TtendAdvr, Ttendtot, O StendSurf, StendDiffh, StendDiffr, O StendAdvh, StendAdvr, Stendtot, O Hc, PIc, I myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE LAYERS_DIAPYCNAL C | Calculate the diapycnal velocity in isotracer layers, for a chosen C | tracer. C *==========================================================* C \ev IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "LAYERS_SIZE.h" #include "LAYERS.h" C !INPUT PARAMETERS: C myThid :: my Thread Id number C tracer :: potential temperature, salt or potential density prho C iLa :: layer coordinate index C TtendSurf :: temperature tendency due to surface forcing times thickness C TtendDiffh:: temperature tendency due to horiz. diffusion times thickness C TtendDiffr:: temperature tendency due to vert. diffusion times thickness C TtendAdvh:: salinity tendency due to horiz. advection times thickness C TtendAdvr:: salinity tendency due to vert. advection times thickness C StendSurf :: salinity tendency due to surface forcing times thickness C StendDiffh:: salinity tendency due to horiz. diffusion times thickness C StendDiffr:: salinity tendency due to vert. diffusion times thickness C StendAdvh :: salinity tendency due to horiz. advection times thickness C StendAdvr :: salinity tendency due to vert. advection times thickness C Hc :: Layer thickness at the tracer point (m) C PIw :: 1 if layer exists, 0 otherwise (at tracer point) INTEGER iLa, myThid _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy) _RL Hc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL PIc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy) _RL TtendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL TtendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL TtendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL TtendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL TtendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL Ttendtot (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL StendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL StendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL StendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL StendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL StendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) _RL Stendtot (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) CEOP #ifdef LAYERS_THERMODYNAMICS C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C !LOCAL VARIABLES: C bi, bj :: tile indices C i,j :: horizontal indices C k :: vertical index for model grid C kp1 :: vertical index for model grid next cell C kci :: index from CellIndex C kg :: index for looping though layers_bounds C kk :: vertical index for ZZ (fine) grid C kloc :: local copy of kgu/v to reduce accesses to index arrays C mSteps :: maximum number of steps for bisection method C TatC :: temperature at C point _RL Hcw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy) INTEGER bi, bj INTEGER i,j,k,kk,kg,kci,kloc INTEGER mSteps INTEGER kgc(sNx+1,sNy+1) INTEGER kgcw(sNx+1,sNy+1) _RL TatC(sNx+1,sNy+1), dzfac, Tfac, Sfac LOGICAL errorFlag CHARACTER*(MAX_LEN_MBUF) msgBuf #ifdef LAYERS_FINEGRID_DIAPYCNAL INTEGER kp1 #endif C -- constants for T and S forcing, gets reset later for rho Tfac = 1. _d 0 Sfac = 1. _d 0 C compute maximum number of steps for bisection method (approx. C log2(Nlayers)) as log2(Nlayers) + 1 for safety mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1 C STOP 'DEBUG END: S/R LAYERS_DIAPYCNAL' C --- The tile loops DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C Initialize the search indices DO j = 1,sNy+1 DO i = 1,sNx+1 C The temperature index (layer_G) goes from cold to warm. C The water column goes from warm (k=1) to cold (k=Nr). C So initialize the search with the warmest value. kgc(i,j) = Nlayers kgcw(i,j) = Nlayers-1 ENDDO ENDDO C Reset the arrays C --- These are at the w point DO kg=1,Nlayers-1 DO j = 1-OLy,sNy+OLy DO i = 1-OLx,sNx+OLx TtendSurf (i,j,kg,bi,bj) = 0. _d 0 TtendDiffh(i,j,kg,bi,bj) = 0. _d 0 TtendDiffr(i,j,kg,bi,bj) = 0. _d 0 TtendAdvh(i,j,kg,bi,bj) = 0. _d 0 TtendAdvr(i,j,kg,bi,bj) = 0. _d 0 Ttendtot(i,j,kg,bi,bj) = 0. _d 0 StendSurf (i,j,kg,bi,bj) = 0. _d 0 StendDiffh(i,j,kg,bi,bj) = 0. _d 0 StendDiffr(i,j,kg,bi,bj) = 0. _d 0 StendAdvh(i,j,kg,bi,bj) = 0. _d 0 StendAdvr(i,j,kg,bi,bj) = 0. _d 0 Stendtot(i,j,kg,bi,bj) = 0. _d 0 Hcw(i,j,kg,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO C --- These are at the c point DO kg=1,Nlayers DO j = 1-OLy,sNy+OLy DO i = 1-OLx,sNx+OLx Hc(i,j,kg,bi,bj) = 0. _d 0 PIc(i,j,kg,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO #ifdef LAYERS_FINEGRID_DIAPYCNAL DO kk=1,NZZ k = MapIndex(kk) kci = CellIndex(kk) DO j = 1,sNy+1 DO i = 1,sNx+1 C ------ Find theta at the V point (south) on the fine Z grid kp1=k+1 IF (maskC(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k TatC(i,j) = MapFact(kk) * tracer(i,j,k,bi,bj) + & (1. _d 0 -MapFact(kk)) * tracer(i,j,kp1,bi,bj) ENDDO ENDDO #else DO kk=1,Nr k = kk kci = kk DO j = 1,sNy+1 DO i = 1,sNx+1 TatC(i,j) = tracer(i,j,k,bi,bj) ENDDO ENDDO #endif /* LAYERS_FINEGRID_DIAPYCNAL */ C ------ debugging stuff c IF (i.EQ.38 .AND. j.EQ.4 .AND. bi.EQ.1 .AND. bj.EQ.1) THEN c i=38 c j=4 c WRITE(msgBuf, c & '(A,I3,A,I3,A,I3,A,F7.2,A,F7.2,A,F7.2,A,F7.2,A,F3.1)') c & 'LAYERS_DEBUG: iLa=', iLa, c & ', kk=', kk, c & ', k=', k, c & ', tracer=', tracer(i,j,k,bi,bj), c & ', TatC=',TatC(i,j), c & ', hFacC=',hFacC(i,j,k,bi,bj) c CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, c & SQUEEZE_RIGHT, myThid ) c ENDIF C ------ Now that we know T everywhere, determine the binning. C find the layer indices kgc for the center point CALL LAYERS_LOCATE( I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatC, O kgc, I myThid ) #ifndef TARGET_NEC_SX C check for failures IF ( debugLevel .GE. debLevC ) THEN errorFlag = .FALSE. DO j = 1,sNy+1 DO i = 1,sNx+1 IF ( kgc(i,j) .LE. 0 ) THEN WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)') & 'S/R LAYERS_LOCATE: Could not find a bin in ', & 'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j) CALL PRINT_ERROR( msgBuf, myThid ) errorFlag = .TRUE. ENDIF ENDDO ENDDO IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL' ENDIF #endif /* ndef TARGET_NEC_SX */ C find the layer indices kgcw for the w point CALL LAYERS_LOCATE( I layers_bounds_w(1,iLa),Nlayers-1,mSteps,sNx,sNy,TatC, O kgcw, I myThid ) #ifndef TARGET_NEC_SX C check for failures IF ( debugLevel .GE. debLevC ) THEN errorFlag = .FALSE. DO j = 1,sNy+1 DO i = 1,sNx+1 IF ( kgcw(i,j) .LE. 0 ) THEN WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)') & 'S/R LAYERS_LOCATE: Could not find a bin in ', & 'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j) CALL PRINT_ERROR( msgBuf, myThid ) errorFlag = .TRUE. ENDIF ENDDO ENDDO IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL' ENDIF #endif /* ndef TARGET_NEC_SX */ C ------ Augment the bin values DO j = 1,sNy+1 DO i = 1,sNx+1 #ifdef LAYERS_FINEGRID_DIAPYCNAL dzfac = dZZf(kk) * hFacC(i,j,kci,bi,bj) #else dzfac = dRf(kci) * hFacC(i,j,kci,bi,bj) #endif /* LAYERS_FINEGRID_DIAPYCNAL */ kloc = kgcw(i,j) C ------- Thickness at w point Hcw(i,j,kloc,bi,bj) = Hcw(i,j,kloc,bi,bj) & + dzfac C ------- Thickness at c point Hc(i,j,kgc(i,j),bi,bj) = Hc(i,j,kgc(i,j),bi,bj) & + dzfac C ------- Now rescale dzfac to include the layer coordinate spacing dzfac = dzfac * layers_recip_delta(kloc,iLa) IF ( layers_num(iLa) .EQ. 3 ) THEN Tfac = layers_alpha(i,j,kci,bi,bj) Sfac = layers_beta(i,j,kci,bi,bj) ENDIF IF (kci.EQ.1) THEN C ------- We are in the surface layer TtendSurf(i,j,kloc,bi,bj) = & TtendSurf(i,j,kloc,bi,bj) + & Tfac * dzfac * layers_surfflux(i,j,1,1,bi,bj) StendSurf(i,j,kloc,bi,bj) = & StendSurf(i,j,kloc,bi,bj) + & Sfac * dzfac * layers_surfflux(i,j,1,2,bi,bj) ENDIF #ifdef SHORTWAVE_HEATING TtendSurf(i,j,kloc,bi,bj) = & TtendSurf(i,j,kloc,bi,bj) + & Tfac * dzfac * layers_sw(i,j,kci,1,bi,bj) #endif /* SHORTWAVE_HEATING */ C ------- Diffusion TtendDiffh(i,j,kloc,bi,bj) = & TtendDiffh(i,j,kloc,bi,bj) + dzfac * Tfac * & (layers_dfx(i,j,kci,1,bi,bj)+ & layers_dfy(i,j,kci,1,bi,bj)) TtendDiffr(i,j,kloc,bi,bj) = & TtendDiffr(i,j,kloc,bi,bj) + & dzfac * Tfac * layers_dfr(i,j,kci,1,bi,bj) StendDiffh(i,j,kloc,bi,bj) = & StendDiffh(i,j,kloc,bi,bj) + dzfac * Sfac * & (layers_dfx(i,j,kci,2,bi,bj)+ & layers_dfy(i,j,kci,2,bi,bj)) StendDiffr(i,j,kloc,bi,bj) = & StendDiffr(i,j,kloc,bi,bj) + & dzfac * Sfac * layers_dfr(i,j,kci,2,bi,bj) C ------- Advection TtendAdvh(i,j,kloc,bi,bj) = & TtendAdvh(i,j,kloc,bi,bj) + dzfac * Tfac * & (layers_afx(i,j,kci,1,bi,bj)+ & layers_afy(i,j,kci,1,bi,bj)) TtendAdvr(i,j,kloc,bi,bj) = & TtendAdvr(i,j,kloc,bi,bj) + & dzfac * Tfac * layers_afr(i,j,kci,1,bi,bj) StendAdvh(i,j,kloc,bi,bj) = & StendAdvh(i,j,kloc,bi,bj) + dzfac * Sfac * & (layers_afx(i,j,kci,2,bi,bj)+ & layers_afy(i,j,kci,2,bi,bj)) StendAdvr(i,j,kloc,bi,bj) = & StendAdvr(i,j,kloc,bi,bj) + & dzfac * Sfac * layers_afr(i,j,kci,2,bi,bj) C -------- Total Tendency Ttendtot(i,j,kloc,bi,bj) = & Ttendtot(i,j,kloc,bi,bj) + & dzfac * Tfac * layers_tottend(i,j,kci,1,bi,bj) Stendtot(i,j,kloc,bi,bj) = & Stendtot(i,j,kloc,bi,bj) + & dzfac * Sfac * layers_tottend(i,j,kci,2,bi,bj) ENDDO ENDDO ENDDO C-- Now that we know the thicknesses, compute the heaviside function C-- (Needs another loop through Ng) DO kg=1,Nlayers DO j = 1,sNy+1 DO i = 1,sNx+1 IF (Hc(i,j,kg,bi,bj) .GT. 0.) THEN PIc(i,j,kg,bi,bj) = 1. _d 0 ENDIF ENDDO ENDDO ENDDO C --- End bi,bj loop ENDDO ENDDO #endif /* LAYERS_THERMODYNAMICS */ RETURN END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE LAYERS_LOCATE( I xx,n,m,sNx,sNy,x, O k, I myThid ) C !DESCRIPTION: \bv C *==========================================================* C | Find the index(-array) k such that x is between xx(k) C | and xx(k+1) by bisection, following Press et al., C | Numerical Recipes in Fortran. xx must be monotonic. C *==========================================================* C \ev C !USES: IMPLICIT NONE C !INPUT PARAMETERS: C xx :: array of bin-boundaries (layers_boundaries) C n :: length of xx C m :: int(log2(n)) + 1 = length of bisection loop C sNx,sNy :: size of index array and input x C x :: input array of values C k :: index array (output) C myThid :: my Thread Id number INTEGER n,m,sNx,sNy _RL xx(1:n+1) _RL x(snx+1,sny+1) INTEGER k(snx+1,sny+1) INTEGER myThid C !LOCAL VARIABLES: C i,j :: horizontal indices C l :: bisection loop index C kl,ku,km :: work arrays and variables INTEGER i,j CEOP #ifdef TARGET_NEC_SX INTEGER l, km INTEGER kl(sNx+1,sNy+1), ku(sNx+1,sNy+1) C bisection, following Press et al., Numerical Recipes in Fortran, C mostly, because it can be vectorized DO j = 1,sNy+1 DO i = 1,sNx+1 kl(i,j)=1 ku(i,j)=n+1 END
DO END
DO DO l = 1,m DO j = 1,sNy+1 DO i = 1,sNx+1 IF (ku(i,j)-kl(i,j).GT.1) THEN km=(ku(i,j)+kl(i,j))/2 CML IF ((xx(n).GE.xx(1)).EQV.(x(i,j).GE.xx(km))) THEN IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR. & ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN kl(i,j)=km ELSE ku(i,j)=km END
IF END
IF END
DO END
DO END
DO DO j = 1,sNy+1 DO i = 1,sNx+1 IF ( x(i,j).LT.xx(2) ) THEN k(i,j)=1 ELSE IF ( x(i,j).GE.xx(n) ) THEN k(i,j)=n ELSE k(i,j)=kl(i,j) END
IF END
DO END
DO #else C the old way DO j = 1,sNy+1 DO i = 1,sNx+1 IF (x(i,j) .GE. xx(n)) THEN C the point is in the hottest bin or hotter k(i,j) = n ELSE IF (x(i,j) .LT. xx(2)) THEN C the point is in the coldest bin or colder k(i,j) = 1 ELSE IF ( (x(i,j) .GE. xx(k(i,j))) & .AND. (x(i,j) .LT. xx(k(i,j)+1)) ) THEN C already on the right bin -- do nothing ELSE IF (x(i,j) .GE. xx(k(i,j))) THEN C have to hunt for the right bin by getting hotter DO WHILE (x(i,j) .GE. xx(k(i,j)+1)) k(i,j) = k(i,j) + 1 ENDDO C now xx(k) < x <= xx(k+1) ELSE IF (x(i,j) .LT. xx(k(i,j)+1)) THEN C have to hunt for the right bin by getting colder DO WHILE (x(i,j) .LT. xx(k(i,j))) k(i,j) = k(i,j) - 1 ENDDO C now xx(k) <= x < xx(k+1) ELSE C that should have covered all the options k(i,j) = -1 ENDIF ENDDO ENDDO #endif /* TARGET_NEC_SX */ #endif /* ALLOW_LAYERS */ RETURN END