--- MITgcm/model/src/external_forcing_surf.F 2003/10/24 05:29:35 1.16 +++ MITgcm/model/src/external_forcing_surf.F 2009/12/21 00:24:58 1.49 @@ -1,20 +1,20 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/external_forcing_surf.F,v 1.16 2003/10/24 05:29:35 edhill Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/external_forcing_surf.F,v 1.49 2009/12/21 00:24:58 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" - + CBOP C !ROUTINE: EXTERNAL_FORCING_SURF C !INTERFACE: - SUBROUTINE EXTERNAL_FORCING_SURF( + SUBROUTINE EXTERNAL_FORCING_SURF( I bi, bj, iMin, iMax, jMin, jMax, I myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* -C | SUBROUTINE EXTERNAL_FORCING_SURF -C | o Determines forcing terms based on external fields -C | relaxation terms etc. +C | SUBROUTINE EXTERNAL_FORCING_SURF +C | o Determines forcing terms based on external fields +C | relaxation terms etc. C *==========================================================* C \ev @@ -28,102 +28,225 @@ #include "DYNVARS.h" #include "GRID.h" #include "SURFACE.h" - +#ifdef ALLOW_SEAICE +#include "SEAICE_PARAMS.h" +#include "SEAICE.h" +#endif /* ALLOW_SEAICE */ +#ifdef ALLOW_SHELFICE +#include "SHELFICE.h" +#endif /* ALLOW_SHELFICE */ + C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === -C myTime - Current time in simulation -C myIter - Current iteration number in simulation +C bi,bj :: tile indices +C iMin,iMax, jMin,jMax :: Range of points for calculation +C myTime :: Current time in simulation +C myIter :: Current iteration number in simulation C myThid :: Thread no. that called this routine. - _RL myTime - INTEGER myIter - INTEGER myThid INTEGER bi,bj INTEGER iMin, iMax INTEGER jMin, jMax + _RL myTime + INTEGER myIter + INTEGER myThid C !LOCAL VARIABLES: C === Local variables === +C i,j :: loop indices +C ks :: index of surface interface layer INTEGER i,j -C number of surface interface layer - INTEGER kSurface + INTEGER ks CEOP +#ifdef ALLOW_DIAGNOSTICS + _RL tmpFac +#endif /* ALLOW_DIAGNOSTICS */ - if ( buoyancyRelation .eq. 'OCEANICP' ) then - kSurface = Nr - else - kSurface = 1 - endif + IF ( usingPCoords ) THEN + ks = Nr + ELSE + ks = 1 + ENDIF -C-- Surface Fluxes : +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| - DO j = jMin, jMax - DO i = iMin, iMax + IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN +C-- Start with surface restoring term : -c Zonal wind stress fu: - surfaceTendencyU(i,j,bi,bj) = - & fu(i,j,bi,bj)*horiVertRatio*recip_rhoConst - & *recip_drF(kSurface)*recip_hFacW(i,j,kSurface,bi,bj) -c Meridional wind stress fv: - surfaceTendencyV(i,j,bi,bj) = - & fv(i,j,bi,bj)*horiVertRatio*recip_rhoConst - & *recip_drF(kSurface)*recip_hFacS(i,j,kSurface,bi,bj) -c Net heat flux Qnet: - surfaceTendencyT(i,j,bi,bj) = - & -Qnet(i,j,bi,bj)*recip_Cp*horiVertRatio*recip_rhoConst - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) -C Net Salt Flux : - surfaceTendencyS(i,j,bi,bj) = - & -saltFlux(i,j,bi,bj)*horiVertRatio*recip_rhoConst - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) - -#ifdef ALLOW_PASSIVE_TRACER -c *** define the tracer surface tendency here *** -#endif /* ALLOW_PASSIVE_TRACER */ +#ifdef ALLOW_SEAICE + IF ( useSEAICE .AND. (.NOT. SEAICErestoreUnderIce) ) THEN +C Do not restore under sea-ice + DO j = jMin, jMax + DO i = iMin, iMax +C Heat Flux (restoring term) : + surfaceForcingT(i,j,bi,bj) = + & -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,bi,bj)) + & *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj)) + & *drF(ks)*_hFacC(i,j,ks,bi,bj) +C Salt Flux (restoring term) : + surfaceForcingS(i,j,bi,bj) = + & -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,bi,bj)) + & *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj)) + & *drF(ks)*_hFacC(i,j,ks,bi,bj) + ENDDO + ENDDO + ELSE +#endif /* ALLOW_SEAICE */ + DO j = jMin, jMax + DO i = iMin, iMax +C Heat Flux (restoring term) : + surfaceForcingT(i,j,bi,bj) = + & -lambdaThetaClimRelax(i,j,bi,bj) + & *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj)) + & *drF(ks)*_hFacC(i,j,ks,bi,bj) +C Salt Flux (restoring term) : + surfaceForcingS(i,j,bi,bj) = + & -lambdaSaltClimRelax(i,j,bi,bj) + & *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj)) + & *drF(ks)*_hFacC(i,j,ks,bi,bj) + ENDDO + ENDDO +#ifdef ALLOW_SEAICE + ENDIF +#endif /* ALLOW_SEAICE */ +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +#ifdef NONLIN_FRSURF +C- T,S surface forcing will be applied (thermodynamics) after the update +C of surf.thickness (hFac): account for change in surf.thickness + IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN + IF (select_rStar.GT.0) THEN +# ifndef DISABLE_RSTAR_CODE + DO j=jMin,jMax + DO i=iMin,iMax + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) + & * rStarExpC(i,j,bi,bj) + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) + & * rStarExpC(i,j,bi,bj) + ENDDO ENDDO - ENDDO - -C-- Surface restoring term : +# endif /* DISABLE_RSTAR_CODE */ + ELSE + DO j=jMin,jMax + DO i=iMin,iMax + IF (ks.EQ.ksurfC(i,j,bi,bj)) THEN + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) + & *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj) + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) + & *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj) + ENDIF + ENDDO + ENDDO + ENDIF + ENDIF +#endif /* NONLIN_FRSURF */ + +#ifdef ALLOW_DIAGNOSTICS + IF ( useDiagnostics ) THEN + +C tRelax (temperature relaxation [W/m2], positive <-> increasing Theta) + tmpFac = HeatCapacity_Cp*rUnit2mass + CALL DIAGNOSTICS_SCALE_FILL( + & surfaceForcingT(1-OLx,1-OLy,bi,bj),tmpFac,1, + & 'TRELAX ',0, 1,2,bi,bj,myThid) + +C sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt) + tmpFac = rUnit2mass + CALL DIAGNOSTICS_SCALE_FILL( + & surfaceForcingS(1-OLx,1-OLy,bi,bj),tmpFac,1, + & 'SRELAX ',0, 1,2,bi,bj,myThid) + + ENDIF +#endif /* ALLOW_DIAGNOSTICS */ + + ELSE +C-- No restoring for T & S : set surfaceForcingT,S to zero : - IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN DO j = jMin, jMax DO i = iMin, iMax -C Heat Flux (restoring term) : - IF ( abs(yC(i,j,bi,bj)).LE.latBandClimRelax ) THEN - surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj) - & -lambdaThetaClimRelax - & *(theta(i,j,kSurface,bi,bj)-SST(i,j,bi,bj)) -C Salt Flux (restoring term) : - surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj) - & -lambdaSaltClimRelax - & *(salt(i,j,kSurface,bi,bj)-SSS(i,j,bi,bj)) - ENDIF + surfaceForcingT(i,j,bi,bj) = 0. _d 0 + surfaceForcingS(i,j,bi,bj) = 0. _d 0 ENDDO ENDDO + +C-- end restoring / no restoring block. ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| -C-- Fresh-water flux -#ifdef NONLIN_FRSURF - IF ( (nonlinFreeSurf.GT.0 .OR. buoyancyRelation.EQ.'OCEANICP') +C-- Surface Fluxes : + + DO j = jMin, jMax + DO i = iMin, iMax + +C Zonal wind stress fu: + surfaceForcingU(i,j,bi,bj) = + & fu(i,j,bi,bj)*mass2rUnit +C Meridional wind stress fv: + surfaceForcingV(i,j,bi,bj) = + & fv(i,j,bi,bj)*mass2rUnit +C Net heat flux Qnet: + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) + & - ( Qnet(i,j,bi,bj) +#ifdef SHORTWAVE_HEATING + & -Qsw(i,j,bi,bj) +#endif + & ) *recip_Cp*mass2rUnit +C Net Salt Flux : + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) + & -saltFlux(i,j,bi,bj)*mass2rUnit + + ENDDO + ENDDO + +#ifdef ALLOW_SALT_PLUME +C saltPlume is the amount of salt rejected by ice while freezing; +C it is here subtracted from surfaceForcingS and will be redistributed +C to multiple vertical levels later on as per Duffy et al. (GRL 1999) + IF ( useSALT_PLUME ) THEN + CALL SALT_PLUME_FORCING_SURF( + I bi, bj, iMin, iMax, jMin, jMax, + I myTime,myIter,myThid ) + ENDIF +#endif /* ALLOW_SALT_PLUME */ + +C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| +C-- Fresh-water flux + +C- Apply mask on Fresh-Water flux +C (needed for SSH forcing, whether or not exactConserv is used) + IF ( useRealFreshWaterFlux ) THEN + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + EmPmR(i,j,bi,bj) = EmPmR(i,j,bi,bj)*maskInC(i,j,bi,bj) + ENDDO + ENDDO + ENDIF + +#ifdef EXACT_CONSERV +C NB: synchronous time step: PmEpR lag 1 time step behind EmPmR +C to stay consitent with volume change (=d/dt etaH). + IF ( staggerTimeStep ) THEN + DO j=1-OLy,sNy+OLy + DO i=1-OLx,sNx+OLx + PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj) + ENDDO + ENDDO + ENDIF + + IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords) & .AND. useRealFreshWaterFlux ) THEN -c- NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes -c the water column height ; temp., salt, (tracer) flux associated -c with this input/output of water is added here to the surface tendency. -c -c NB: PmEpR lag 1 time step behind EmPmR ( PmEpR_n = - EmPmR_n-1 ) to stay -c consitent with volume change (=d/dt etaN). +C-- NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes +C the water column height ; temp., salt, (tracer) flux associated +C with this input/output of water is added here to the surface tendency. IF (temp_EvPrRn.NE.UNSET_RL) THEN DO j = jMin, jMax DO i = iMin, iMax - surfaceTendencyT(i,j,bi,bj) = surfaceTendencyT(i,j,bi,bj) + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) & + PmEpR(i,j,bi,bj) - & *( temp_EvPrRn - theta(i,j,kSurface,bi,bj) ) - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) - & *convertEmP2rUnit + & *( temp_EvPrRn - theta(i,j,ks,bi,bj) ) + & *mass2rUnit ENDDO ENDDO ENDIF @@ -131,87 +254,152 @@ IF (salt_EvPrRn.NE.UNSET_RL) THEN DO j = jMin, jMax DO i = iMin, iMax - surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj) + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) & + PmEpR(i,j,bi,bj) - & *( salt_EvPrRn - salt(i,j,kSurface,bi,bj) ) - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) - & *convertEmP2rUnit + & *( salt_EvPrRn - salt(i,j,ks,bi,bj) ) + & *mass2rUnit ENDDO ENDDO ENDIF -#ifdef ALLOW_PASSIVE_TRACER -c *** add the tracer flux associated with P-E+R here *** -c IF (trac_EvPrRn.NE.UNSET_RL) THEN -c & + PmEpR(i,j,bi,bj)*( trac_EvPrRn - tr1(i,j,kSurface,bi,bj) ) -c & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) -c ENDIF -#endif /* ALLOW_PASSIVE_TRACER */ - C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| ELSE -#else /* NONLIN_FRSURF */ +#else /* EXACT_CONSERV */ IF (.TRUE.) THEN -#endif /* NONLIN_FRSURF */ +#endif /* EXACT_CONSERV */ -c- EmPmR does not really affect the water column height (for tracer budget) -c and is converted to a salt tendency. +C-- EmPmR does not really affect the water column height (for tracer budget) +C and is converted to a salt tendency. IF (convertFW2Salt .EQ. -1.) THEN -c- converts EmPmR to salinity tendency using surface local salinity - DO j = jMin, jMax - DO i = iMin, iMax - surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj) - & + EmPmR(i,j,bi,bj)*salt(i,j,kSurface,bi,bj) - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) - & *convertEmP2rUnit +C- use local surface tracer field to calculate forcing term: + + IF (temp_EvPrRn.NE.UNSET_RL) THEN +C account for Rain/Evap heat content (temp_EvPrRn) using local SST + DO j = jMin, jMax + DO i = iMin, iMax + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) + & + EmPmR(i,j,bi,bj) + & *( theta(i,j,ks,bi,bj) - temp_EvPrRn ) + & *mass2rUnit + ENDDO ENDDO - ENDDO - ELSE -c- converts EmPmR to virtual salt flux using uniform salinity (default=35) - DO j = jMin, jMax - DO i = iMin, iMax - surfaceTendencyS(i,j,bi,bj) = surfaceTendencyS(i,j,bi,bj) - & + EmPmR(i,j,bi,bj)*convertFW2Salt - & *recip_drF(kSurface)*recip_hFacC(i,j,kSurface,bi,bj) - & *convertEmP2rUnit + ENDIF + IF (salt_EvPrRn.NE.UNSET_RL) THEN +C converts EmPmR to salinity tendency using surface local salinity + DO j = jMin, jMax + DO i = iMin, iMax + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) + & + EmPmR(i,j,bi,bj) + & *( salt(i,j,ks,bi,bj) - salt_EvPrRn ) + & *mass2rUnit + ENDDO ENDDO - ENDDO + ENDIF + + ELSE +C- use uniform tracer value to calculate forcing term: + + IF (temp_EvPrRn.NE.UNSET_RL) THEN +C account for Rain/Evap heat content (temp_EvPrRn) assuming uniform SST (=tRef) + DO j = jMin, jMax + DO i = iMin, iMax + surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj) + & + EmPmR(i,j,bi,bj) + & *( tRef(ks) - temp_EvPrRn ) + & *mass2rUnit + ENDDO + ENDDO + ENDIF + IF (salt_EvPrRn.NE.UNSET_RL) THEN +C converts EmPmR to virtual salt flux using uniform salinity (default=35) + DO j = jMin, jMax + DO i = iMin, iMax + surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj) + & + EmPmR(i,j,bi,bj) + & *( convertFW2Salt - salt_EvPrRn ) + & *mass2rUnit + ENDDO + ENDDO + ENDIF + +C- end local-surface-tracer / uniform-value distinction ENDIF ENDIF + C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| + #ifdef ALLOW_PTRACERS IF ( usePTRACERS ) THEN CALL PTRACERS_FORCING_SURF( - I bi, bj, iMin, iMax, jMin, jMax, - I myTime,myIter,myThid ) + I bi, bj, iMin, iMax, jMin, jMax, + I myTime,myIter,myThid ) ENDIF #endif /* ALLOW_PTRACERS */ #ifdef ATMOSPHERIC_LOADING +C-- Atmospheric surface Pressure loading : added to phi0surf when using Z-coord; +C Not yet implemented for Ocean in P: would need to be applied to the other end +C of the column, as a vertical velocity (omega); (meaningless for Atmos in P). +C- Note: +C Using P-coord., a hack (now directly applied from S/R INI_FORCING) +C is sometime used to read phi0surf from a file (pLoadFile) instead +C of computing it from bathymetry & density ref. profile. -C-- Atmospheric surface Pressure loading : - - IF (buoyancyRelation .eq. 'OCEANIC' ) THEN + IF ( usingZCoords ) THEN +C The true atmospheric P-loading is not yet implemented for P-coord +C (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS). + IF ( useRealFreshWaterFlux ) THEN + DO j = jMin, jMax + DO i = iMin, iMax + phi0surf(i,j,bi,bj) = ( pLoad(i,j,bi,bj) + & +sIceLoad(i,j,bi,bj)*gravity + & )*recip_rhoConst + ENDDO + ENDDO + ELSE DO j = jMin, jMax DO i = iMin, iMax - phi0surf(i,j,bi,bj) = pload(i,j,bi,bj)*recip_rhoConst + phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)*recip_rhoConst ENDDO ENDDO - ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN -C-- This is a hack used to read phi0surf from a file (ploadFile) + ENDIF +c ELSEIF ( usingPCoords ) THEN +C-- This is a hack used to read phi0surf from a file (pLoadFile) C instead of computing it from bathymetry & density ref. profile. -C The true atmospheric P-loading is not yet implemented for P-coord -C (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS). +C ==> now done only once, in S/R INI_FORCING +c DO j = jMin, jMax +c DO i = iMin, iMax +c phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj) +c ENDDO +c ENDDO + ENDIF +#endif /* ATMOSPHERIC_LOADING */ + +#ifdef ALLOW_SHELFICE + IF ( usingZCoords ) THEN + IF ( useSHELFICE) THEN DO j = jMin, jMax DO i = iMin, iMax - phi0surf(i,j,bi,bj) = pload(i,j,bi,bj) + phi0surf(i,j,bi,bj) = phi0surf(i,j,bi,bj) + & + shelficeLoadAnomaly(i,j,bi,bj)*recip_rhoConst ENDDO ENDDO + ENDIF ENDIF +#endif /* ALLOW_SHELFICE */ -#endif /* ATMOSPHERIC_LOADING */ +#ifdef ALLOW_EBM +c-- Values for surfaceForcingT, surfaceForcingS +c are overwritten by those produced by EBM +cph AD recomputation problems if these IF useEBM are used +cph IF ( useEBM ) THEN + CALL EBM_FORCING_SURF( + I bi, bj, iMin, iMax, jMin, jMax, + I myTime,myIter,myThid ) +cph ENDIF +#endif RETURN END