C $Header: /u/gcmpack/MITgcm/pkg/frazil/frazil_calc_rhs.F,v 1.4 2012/03/26 01:47:08 dimitri Exp $ C $Name: $ #include "FRAZIL_OPTIONS.h" CBOP C !ROUTINE: FRAZIL_CALC_RHS C !INTERFACE: ========================================================== SUBROUTINE FRAZIL_CALC_RHS( I myTime, myIter, myThid ) C !DESCRIPTION: C Check water temperature and if colder than freezing C point bring excess negative heat to the surface. C !USES: =============================================================== IMPLICIT NONE #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" #include "FFIELDS.h" #ifdef ALLOW_FRAZIL # include "FRAZIL.h" #endif C !INPUT/OUTPUT PARAMETERS: C myTime :: current time in simulation C myIter :: current iteration number in simulation C myThid :: my Thread Id number _RL myTime INTEGER myIter INTEGER myThid CEOP #ifdef ALLOW_FRAZIL C !LOCAL VARIABLES: C Tfreezing :: freezing threshold temperature INTEGER bi,bj,i,j,k,kTop _RL Tfreezing, Tresid, pLoc, sLoc, tLoc _RL a0, a1, a2, b PARAMETER( a0 = -0.0575 _d 0 ) PARAMETER( a1 = 1.710523 _d -3 ) PARAMETER( a2 = -2.154996 _d -4 ) PARAMETER( b = -7.53 _d -4 ) _RL SW_TEMP EXTERNAL DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C Initialize FrazilForcingT to zero DO k=1,Nr DO j=1-Oly,sNy+OLy DO i=1-Olx,sNx+Olx FrazilForcingT(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO C Check for water below freezing point. DO k = 2, Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx IF ( maskC(i,j,k-1,bi,bj) .NE. 0. _d 0 .AND. & maskC(i,j,k, bi,bj) .NE. 0. _d 0 ) THEN pLoc = ABS(RC(k)) sLoc = MAX(salt(i,j,k,bi,bj), 0. _d 0) tLoc = SW_TEMP(sLoc,theta(i,j,k,bi,bj),pLoc,0. _d 0) C Freezing point of seawater C REFERENCE: UNESCO TECH. PAPERS IN THE MARINE SCIENCE NO. 28. 1978 C EIGHTH REPORT JPOTS C ANNEX 6 FREEZING POINT OF SEAWATER F.J. MILLERO PP.29-35. C C UNITS: C PRESSURE P DECIBARS C SALINITY S PSS-78 C TEMPERATURE TF DEGREES CELSIUS C FREEZING PT. C************************************************************ C CHECKVALUE: TF= -2.588567 DEG. C FOR S=40.0, P=500. DECIBARS Tfreezing = (a0 + a1*sqrt(sLoc) + a2*sLoc) * sLoc + b*pLoc IF (tLoc .LT. Tfreezing) THEN C Move the negative heat to surface level. kTop = kSurfC(i,j,bi,bj) Tresid = ( Tfreezing - tloc ) & * HeatCapacity_Cp * rUnit2mass & * drF(k) * _hFacC(i,j,k,bi,bj) FrazilForcingT(i,j,k,bi,bj) = Tresid / dTtracerLev(k) FrazilForcingT(i,j,kTop,bi,bj) = & FrazilForcingT(i,j,kTop,bi,bj) & - Tresid / dTtracerLev(kTop) ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO # ifdef ALLOW_DIAGNOSTICS IF ( useDiagnostics ) THEN CALL DIAGNOSTICS_FILL( FrazilForcingT, 'FrzForcT', & 0,Nr, 0, 1, 1, myThid ) ENDIF # endif /* ALLOW_DIAGNOSTICS */ #endif /* ALLOW_FRAZIL */ RETURN END