/[MITgcm]/MITgcm/pkg/layers/layers_thermodynamics.F
ViewVC logotype

Annotation of /MITgcm/pkg/layers/layers_thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (hide annotations) (download)
Tue Jul 8 19:03:39 2014 UTC (9 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64z, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65l, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65
Changes since 1.2: +20 -21 lines
- add CVS Header and Name
- remove unused variables

1 jmc 1.3 C $Header: $
2     C $Name: $
3    
4 rpa 1.1 #include "LAYERS_OPTIONS.h"
5     C-- File layers_thermodynamics.F:
6     C-- Contents
7     C-- o LAYERS_CALC_RHS
8    
9     CBOP 0
10     C !ROUTINE: LAYERS_CALC_RHS
11     C !INTERFACE:
12     SUBROUTINE LAYERS_CALC_RHS(
13     I myThid )
14    
15     C !DESCRIPTION: \bv
16     C *==========================================================*
17     C | SUBROUTINE LAYERS_CALC_RHS
18     C | Recalculate the divergence of the RHS terms in T and S eqns.
19     C | Replaces the values of layers_surfflux, layers_df? IN PLACE
20     C | with the corresponding tendencies (same units as GT and GS)
21     C *==========================================================*
22     C \ev
23    
24     C !USES:
25     IMPLICIT NONE
26     C == Global variables ===
27     #include "SIZE.h"
28     #include "EEPARAMS.h"
29     #include "PARAMS.h"
30     #include "GRID.h"
31     #include "LAYERS_SIZE.h"
32     #include "LAYERS.h"
33    
34     C !INPUT PARAMETERS:
35     C myThid :: my Thread Id number
36     INTEGER myThid
37     CEOP
38    
39     #ifdef ALLOW_LAYERS
40 rpa 1.2 #ifdef LAYERS_THERMODYNAMICS
41 rpa 1.1 C !LOCAL VARIABLES:
42     C bi, bj :: tile indices
43     C i,j :: horizontal indices
44 jmc 1.3 C k :: vertical index for model grid
45 rpa 1.1 C kdown :: temporary placeholder
46     C fluxfac :: scaling factor for converting surface flux to tendency
47 jmc 1.3 C fluxfac :: scaling factor for converting diffusive flux to tendency
48 rpa 1.1 C downfac :: mask for lower point
49    
50     INTEGER bi, bj
51 jmc 1.3 INTEGER i,j,k,kdown,iTracer
52 rpa 1.1 _RL fluxfac(2), downfac, tmpfac
53    
54     C -- These factors convert the units of TFLUX and SFLUX diagnostics
55     C -- back to surfaceForcingT and surfaceForcingS units
56     fluxfac(1) = 1.0/(HeatCapacity_Cp*rUnit2mass)
57     fluxfac(2) = 1.0/rUnit2mass
58    
59     DO bj = myByLo(myThid), myByHi(myThid)
60 jmc 1.3 DO bi = myBxLo(myThid), myBxHi(myThid)
61     DO iTracer = 1,2
62 rpa 1.1 k = 1
63 jmc 1.3 C -- Loop for surface fluxes
64 rpa 1.1 DO j=1-OLy,sNy+OLy
65     DO i=1-OLx,sNx+OLx
66     layers_surfflux(i,j,k,iTracer,bi,bj) =
67     & layers_surfflux(i,j,k,iTracer,bi,bj)
68 jmc 1.3 & *recip_drF(1)*_recip_hFacC(i,j,1,bi,bj)
69 rpa 1.1 & *fluxfac(iTracer)
70     ENDDO
71 jmc 1.3 ENDDO
72 rpa 1.1 C -- Loop for diffusive fluxes
73     C -- If done correctly, we can overwrite the flux array in place
74     C -- with its own divergence
75     DO k=1,Nr
76     kdown= MIN(k+1,Nr)
77     IF (k.EQ.Nr) THEN
78     downfac = 0. _d 0
79     ELSE
80     downfac = 1. _d 0
81 jmc 1.3 ENDIF
82 rpa 1.1 DO j=1-OLy,sNy+OLy
83     DO i=1-OLx,sNx+OLx
84     tmpfac = -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
85     & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
86     layers_dfx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
87 jmc 1.3 & tmpfac * ( layers_dfx(i+1,j,k,iTracer,bi,bj) -
88 rpa 1.1 & layers_dfx(i,j,k,iTracer,bi,bj) )
89     layers_dfy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
90 jmc 1.3 & tmpfac * ( layers_dfy(i,j+1,k,iTracer,bi,bj) -
91     & layers_dfy(i,j,k,iTracer,bi,bj) )
92 rpa 1.1 layers_dfr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
93     & ( layers_dfr(i,j,kdown,iTracer,bi,bj)*downfac -
94     & layers_dfr(i,j,k,iTracer,bi,bj) )
95     ENDDO
96 jmc 1.3 ENDDO
97 rpa 1.1 ENDDO
98     ENDDO
99     ENDDO
100     ENDDO
101    
102     C- TFLUX (=total heat flux, match heat-content variations, [W/m2])
103     C IF ( fluidIsWater .AND.
104     C & DIAGNOSTICS_IS_ON('TFLUX ',myThid) ) THEN
105     C DO bj = myByLo(myThid), myByHi(myThid)
106     C DO bi = myBxLo(myThid), myBxHi(myThid)
107     C DO j = 1,sNy
108     C DO i = 1,sNx
109     C tmp1k(i,j,bi,bj) =
110     C #ifdef SHORTWAVE_HEATING
111     C & -Qsw(i,j,bi,bj)+
112     C #endif
113     C & (surfaceForcingT(i,j,bi,bj)+surfaceForcingTice(i,j,bi,bj))
114     C & *HeatCapacity_Cp*rUnit2mass
115     C ENDDO
116     C ENDDO
117     C #ifdef NONLIN_FRSURF
118     C IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
119     C & .AND. useRealFreshWaterFlux ) THEN
120     C DO j=1,sNy
121     C DO i=1,sNx
122     C tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
123     C & + PmEpR(i,j,bi,bj)*theta(i,j,ks,bi,bj)*HeatCapacity_Cp
124     C ENDDO
125     C ENDDO
126     C ENDIF
127     C #endif /* NONLIN_FRSURF */
128     C ENDDO
129     C ENDDO
130     C CALL DIAGNOSTICS_FILL( tmp1k,'TFLUX ',0,1,0,1,1,myThid )
131     C ENDIF
132 jmc 1.3 C
133 rpa 1.1 C C- SFLUX (=total salt flux, match salt-content variations [g/m2/s])
134     C IF ( fluidIsWater .AND.
135     C & DIAGNOSTICS_IS_ON('SFLUX ',myThid) ) THEN
136     C DO bj = myByLo(myThid), myByHi(myThid)
137     C DO bi = myBxLo(myThid), myBxHi(myThid)
138     C DO j = 1,sNy
139     C DO i = 1,sNx
140     C tmp1k(i,j,bi,bj) =
141     C & surfaceForcingS(i,j,bi,bj)*rUnit2mass
142     C ENDDO
143     C ENDDO
144 jmc 1.3 C
145 rpa 1.1 C #ifdef NONLIN_FRSURF
146     C IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
147     C & .AND. useRealFreshWaterFlux ) THEN
148     C DO j=1,sNy
149     C DO i=1,sNx
150     C tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
151     C & + PmEpR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
152     C ENDDO
153     C ENDDO
154     C ENDIF
155     C #endif /* NONLIN_FRSURF */
156 jmc 1.3 C
157 rpa 1.1 C ENDDO
158     C ENDDO
159     C CALL DIAGNOSTICS_FILL( tmp1k,'SFLUX ',0,1,0,1,1,myThid )
160     C ENDIF
161    
162     C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
163     C IF ( kLev .EQ. kSurface ) THEN
164     C DO j=1,sNy
165     C DO i=1,sNx
166     C gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
167     C & +surfaceForcingT(i,j,bi,bj)
168     C & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
169     C ENDDO
170     C ENDDO
171     C ELSEIF ( kSurface.EQ.-1 ) THEN
172     C DO j=1,sNy
173     C DO i=1,sNx
174     C IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
175     C gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
176     C & +surfaceForcingT(i,j,bi,bj)
177     C & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
178     C ENDIF
179     C ENDDO
180     C ENDDO
181     C ENDIF
182    
183     C-- Divergence of fluxes
184     C Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
185     C for Stevens OBC: keep only vertical diffusive contribution on boundaries
186     C DO j=1-OLy,sNy+OLy-1
187     C DO i=1-OLx,sNx+OLx-1
188     C gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
189     C & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
190     C & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
191     C & *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
192     C & +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
193     C & +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
194     C & -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
195     C & +(vTrans(i,j+1)-vTrans(i,j))*advFac
196     C & +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
197     C & )*maskInC(i,j,bi,bj)
198     C & )
199     C ENDDO
200 jmc 1.3 C ENDDO
201 rpa 1.1
202     #endif /* LAYERS_THERMODYNAMICS */
203 rpa 1.2 #endif /* USE_LAYERS */
204 rpa 1.1 RETURN
205     END
206     C -- end of S/R LAYERS_CALC_RHS

  ViewVC Help
Powered by ViewVC 1.1.22