-
Notifications
You must be signed in to change notification settings - Fork 237
/
frazil_calc_rhs.F
115 lines (98 loc) · 3.25 KB
/
frazil_calc_rhs.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#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 SW_TEMP
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