--- MITgcm/model/src/external_forcing.F 2004/10/19 02:39:58 1.28 +++ MITgcm/model/src/external_forcing.F 2005/12/15 17:47:54 1.34 @@ -1,25 +1,25 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/external_forcing.F,v 1.28 2004/10/19 02:39:58 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/external_forcing.F,v 1.34 2005/12/15 17:47:54 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" -#ifdef ALLOW_OBCS -# include "OBCS_OPTIONS.h" +#ifdef ALLOW_EXF +# include "EXF_OPTIONS.h" #endif CBOP C !ROUTINE: EXTERNAL_FORCING_U C !INTERFACE: SUBROUTINE EXTERNAL_FORCING_U( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | S/R EXTERNAL_FORCING_U -C | o Contains problem specific forcing for zonal velocity. +C | S/R EXTERNAL_FORCING_U +C | o Contains problem specific forcing for zonal velocity. C *==========================================================* -C | Adds terms to gU for forcing by external sources -C | e.g. wind stress, bottom friction etc.................. +C | Adds terms to gU for forcing by external sources +C | e.g. wind stress, bottom friction etc ... C *==========================================================* C \ev @@ -35,20 +35,21 @@ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == -C iMin - Working range of tile for applying forcing. -C iMax -C jMin -C jMax -C kLev +C iMin,iMax :: Working range of x-index for applying forcing. +C jMin,jMax :: Working range of y-index for applying forcing. +C bi,bj :: Current tile indices +C kLev :: Current vertical level index +C myTime :: Current time in simulation +C myThid :: Thread Id number INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj - _RL myCurrentTime + _RL myTime INTEGER myThid C !LOCAL VARIABLES: C == Local variables == -C Loop counters - INTEGER I, J -C number of surface interface layer +C i,j :: Loop counters +C kSurface :: index of surface layer + INTEGER i, j INTEGER kSurface CEOP @@ -64,20 +65,21 @@ #ifdef ALLOW_AIM IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_AIM */ -C AMM + #ifdef ALLOW_FIZHI IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_FIZHI */ -C AMM C Add windstress momentum impulse into the top-layer IF ( kLev .EQ. kSurface ) THEN - DO j=jMin,jMax - DO i=iMin,iMax +c DO j=1,sNy +C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1] + DO j=0,sNy+1 + DO i=1,sNx+1 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj) & +foFacMom*surfaceForcingU(i,j,bi,bj) & *recip_drF(kLev)*recip_hFacW(i,j,kLev,bi,bj) @@ -85,29 +87,37 @@ ENDDO ENDIF -#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) +#if (defined (ALLOW_TAU_EDDY)) + CALL TAUEDDY_EXTERNAL_FORCING_U( + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) +#endif + +#ifdef ALLOW_OBCS IF (useOBCS) THEN CALL OBCS_SPONGE_U( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) ENDIF #endif RETURN END + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: EXTERNAL_FORCING_V C !INTERFACE: SUBROUTINE EXTERNAL_FORCING_V( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | S/R EXTERNAL_FORCING_V -C | o Contains problem specific forcing for merid velocity. +C | S/R EXTERNAL_FORCING_V +C | o Contains problem specific forcing for merid velocity. C *==========================================================* -C | Adds terms to gV for forcing by external sources -C | e.g. wind stress, bottom friction etc.................. +C | Adds terms to gV for forcing by external sources +C | e.g. wind stress, bottom friction etc ... C *==========================================================* C \ev @@ -123,20 +133,21 @@ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == -C iMin - Working range of tile for applying forcing. -C iMax -C jMin -C jMax -C kLev +C iMin,iMax :: Working range of x-index for applying forcing. +C jMin,jMax :: Working range of y-index for applying forcing. +C bi,bj :: Current tile indices +C kLev :: Current vertical level index +C myTime :: Current time in simulation +C myThid :: Thread Id number INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj - _RL myCurrentTime + _RL myTime INTEGER myThid C !LOCAL VARIABLES: C == Local variables == -C Loop counters - INTEGER I, J -C number of surface interface layer +C i,j :: Loop counters +C kSurface :: index of surface layer + INTEGER i, j INTEGER kSurface CEOP @@ -152,20 +163,21 @@ #ifdef ALLOW_AIM IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_AIM */ -C AMM #ifdef ALLOW_FIZHI IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_FIZHI */ -C AMM + C Add windstress momentum impulse into the top-layer IF ( kLev .EQ. kSurface ) THEN - DO j=jMin,jMax - DO i=iMin,iMax + DO j=1,sNy+1 +c DO i=1,sNx +C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1] + DO i=0,sNx+1 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj) & +foFacMom*surfaceForcingV(i,j,bi,bj) & *recip_drF(kLev)*recip_hFacS(i,j,kLev,bi,bj) @@ -173,29 +185,37 @@ ENDDO ENDIF -#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) +#if (defined (ALLOW_TAU_EDDY)) + CALL TAUEDDY_EXTERNAL_FORCING_V( + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) +#endif + +#ifdef ALLOW_OBCS IF (useOBCS) THEN CALL OBCS_SPONGE_V( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) ENDIF #endif RETURN END + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: EXTERNAL_FORCING_T C !INTERFACE: SUBROUTINE EXTERNAL_FORCING_T( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | S/R EXTERNAL_FORCING_T -C | o Contains problem specific forcing for temperature. +C | S/R EXTERNAL_FORCING_T +C | o Contains problem specific forcing for temperature. C *==========================================================* -C | Adds terms to gT for forcing by external sources -C | e.g. heat flux, climatalogical relaxation.............. +C | Adds terms to gT for forcing by external sources +C | e.g. heat flux, climatalogical relaxation, etc ... C *==========================================================* C \ev @@ -211,22 +231,23 @@ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == -C iMin - Working range of tile for applying forcing. -C iMax -C jMin -C jMax -C kLev +C iMin,iMax :: Working range of x-index for applying forcing. +C jMin,jMax :: Working range of y-index for applying forcing. +C bi,bj :: Current tile indices +C kLev :: Current vertical level index +C myTime :: Current time in simulation +C myThid :: Thread Id number INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj - _RL myCurrentTime + _RL myTime INTEGER myThid -CEndOfInterface C !LOCAL VARIABLES: C == Local variables == -C Loop counters - INTEGER I, J -C number of surface interface layer +C i,j :: Loop counters +C kSurface :: index of surface layer + INTEGER i, j INTEGER kSurface +CEOP #ifdef SHORTWAVE_HEATING integer two _RL minusone @@ -234,7 +255,6 @@ _RL swfracb(two) INTEGER kp1 #endif -CEOP IF ( fluidIsAir ) THEN kSurface = 0 @@ -248,21 +268,19 @@ #ifdef ALLOW_AIM IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_AIM */ -C AMM #ifdef ALLOW_FIZHI IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_FIZHI */ -C AMM C Add heat in top-layer IF ( kLev .EQ. kSurface ) THEN - DO j=jMin,jMax - DO i=iMin,iMax + DO j=1,sNy + DO i=1,sNx gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj) & +surfaceForcingT(i,j,bi,bj) & *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj) @@ -272,52 +290,69 @@ #ifdef SHORTWAVE_HEATING C Penetrating SW radiation - kp1 = klev+1 - swfracb(1)=abs(rF(klev)) - swfracb(2)=abs(rF(klev+1)) - CALL SWFRAC( +c IF ( usePenetratingSW ) THEN + swfracb(1)=abs(rF(klev)) + swfracb(2)=abs(rF(klev+1)) + CALL SWFRAC( I two,minusone, - I myCurrentTime,myThid, + I myTime,myThid, U swfracb) - IF (klev.EQ.Nr) THEN + kp1 = klev+1 + IF (klev.EQ.Nr) THEN kp1 = klev swfracb(2)=0. _d 0 - ENDIF - DO j=jMin,jMax - DO i=iMin,iMax - gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj) + ENDIF + DO j=1,sNy + DO i=1,sNx + gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj) & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,klev,bi,bj) & -swfracb(2)*maskC(i,j,kp1, bi,bj)) & *recip_Cp*recip_rhoConst & *recip_drF(klev)*recip_hFacC(i,j,kLev,bi,bj) + ENDDO ENDDO - ENDDO +c ENDIF +#endif + +#ifdef ALLOW_CLIMTEMP_RELAXATION + IF ( tauThetaClimRelax3Dim .NE. 0. ) THEN + DO j=1,sNy + DO i=1,sNx + gT(i,j,klev,bi,bj) = gT(i,j,klev,bi,bj) + & -1./tauThetaClimRelax3Dim + & *(theta(i,j,klev,bi,bj)-thetaStar(i,j,klev,bi,bj)) + & *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj) + ENDDO + ENDDO + ENDIF #endif -#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) +#ifdef ALLOW_OBCS IF (useOBCS) THEN CALL OBCS_SPONGE_T( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) ENDIF #endif RETURN END + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: EXTERNAL_FORCING_S C !INTERFACE: SUBROUTINE EXTERNAL_FORCING_S( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | S/R EXTERNAL_FORCING_S -C | o Contains problem specific forcing for merid velocity. +C | S/R EXTERNAL_FORCING_S +C | o Contains problem specific forcing for merid velocity. C *==========================================================* -C | Adds terms to gS for forcing by external sources -C | e.g. fresh-water flux, climatalogical relaxation....... +C | Adds terms to gS for forcing by external sources +C | e.g. fresh-water flux, climatalogical relaxation, etc ... C *==========================================================* C \ev @@ -333,20 +368,21 @@ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == -C iMin - Working range of tile for applying forcing. -C iMax -C jMin -C jMax -C kLev +C iMin,iMax :: Working range of x-index for applying forcing. +C jMin,jMax :: Working range of y-index for applying forcing. +C bi,bj :: Current tile indices +C kLev :: Current vertical level index +C myTime :: Current time in simulation +C myThid :: Thread Id number INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj - _RL myCurrentTime + _RL myTime INTEGER myThid C !LOCAL VARIABLES: C == Local variables == -C Loop counters - INTEGER I, J -C number of surface interface layer +C i,j :: Loop counters +C kSurface :: index of surface layer + INTEGER i, j INTEGER kSurface CEOP @@ -362,21 +398,19 @@ #ifdef ALLOW_AIM IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_AIM */ -C AMM #ifdef ALLOW_FIZHI IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S( & iMin,iMax, jMin,jMax, bi,bj, kLev, - & myCurrentTime, myThid ) + & myTime, myThid ) #endif /* ALLOW_FIZHI */ -C AMM C Add fresh-water in top-layer IF ( kLev .EQ. kSurface ) THEN - DO j=jMin,jMax - DO i=iMin,iMax + DO j=1,sNy + DO i=1,sNx gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj) & +surfaceForcingS(i,j,bi,bj) & *recip_drF(kLev)*recip_hFacC(i,j,kLev,bi,bj) @@ -384,11 +418,24 @@ ENDDO ENDIF -#if (defined (ALLOW_OBCS) && defined (ALLOW_OBCS_SPONGE)) +#ifdef ALLOW_CLIMSALT_RELAXATION + IF ( tauSaltClimRelax3Dim .NE. 0. ) THEN + DO j=1,sNy + DO i=1,sNx + gS(i,j,klev,bi,bj) = gS(i,j,klev,bi,bj) + & -1./tauSaltClimRelax3Dim + & *(salt(i,j,klev,bi,bj)-saltStar(i,j,klev,bi,bj)) + & *hFacC(i,j,klev,bi,bj)*recip_hFacC(i,j,kLev,bi,bj) + ENDDO + ENDDO + ENDIF +#endif + +#ifdef ALLOW_OBCS IF (useOBCS) THEN CALL OBCS_SPONGE_S( - I iMin, iMax, jMin, jMax,bi,bj,kLev, - I myCurrentTime,myThid) + I iMin,iMax, jMin,jMax, bi,bj, kLev, + I myTime, myThid ) ENDIF #endif