1 |
C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.68 2014/05/22 22:00:36 atn Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "CPP_OPTIONS.h" |
5 |
|
6 |
C-- File apply_forcing.F: |
7 |
C-- Contents |
8 |
C-- o APPLY_FORCING_U |
9 |
C-- o APPLY_FORCING_V |
10 |
C-- o APPLY_FORCING_T |
11 |
C-- o APPLY_FORCING_S |
12 |
|
13 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
14 |
CBOP |
15 |
C !ROUTINE: APPLY_FORCING_U |
16 |
C !INTERFACE: |
17 |
SUBROUTINE APPLY_FORCING_U( |
18 |
U gU_arr, |
19 |
I iMin,iMax,jMin,jMax, k, bi, bj, |
20 |
I myTime, myIter, myThid ) |
21 |
C !DESCRIPTION: \bv |
22 |
C *==========================================================* |
23 |
C | S/R APPLY_FORCING_U |
24 |
C | o Contains problem specific forcing for zonal velocity. |
25 |
C *==========================================================* |
26 |
C | Adds terms to gU for forcing by external sources |
27 |
C | e.g. wind stress, bottom friction etc ... |
28 |
C *==========================================================* |
29 |
C \ev |
30 |
|
31 |
C !USES: |
32 |
IMPLICIT NONE |
33 |
C == Global data == |
34 |
#include "SIZE.h" |
35 |
#include "EEPARAMS.h" |
36 |
#include "PARAMS.h" |
37 |
#include "GRID.h" |
38 |
#include "DYNVARS.h" |
39 |
#include "FFIELDS.h" |
40 |
|
41 |
C !INPUT/OUTPUT PARAMETERS: |
42 |
C gU_arr :: the tendency array |
43 |
C iMin,iMax :: Working range of x-index for applying forcing. |
44 |
C jMin,jMax :: Working range of y-index for applying forcing. |
45 |
C k :: Current vertical level index |
46 |
C bi,bj :: Current tile indices |
47 |
C myTime :: Current time in simulation |
48 |
C myIter :: Current iteration number |
49 |
C myThid :: my Thread Id number |
50 |
_RL gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
51 |
INTEGER iMin, iMax, jMin, jMax |
52 |
INTEGER k, bi, bj |
53 |
_RL myTime |
54 |
INTEGER myIter |
55 |
INTEGER myThid |
56 |
|
57 |
C !LOCAL VARIABLES: |
58 |
C i,j :: Loop counters |
59 |
INTEGER i, j |
60 |
CEOP |
61 |
_RL recip_P0g, termP, kV, kF, sigma_b |
62 |
|
63 |
C-- Forcing term |
64 |
kF = 1. _d 0/86400. _d 0 |
65 |
sigma_b = 0.7 _d 0 |
66 |
c DO j=1,sNy |
67 |
C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1] |
68 |
DO j=0,sNy+1 |
69 |
DO i=1,sNx+1 |
70 |
IF ( maskW(i,j,k,bi,bj).EQ.oneRS ) THEN |
71 |
recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i-1,j,bi,bj)) |
72 |
termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL ) |
73 |
& +rF(k+1)*recip_P0g ) |
74 |
c termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g |
75 |
kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) ) |
76 |
gU_arr(i,j) = gU_arr(i,j) |
77 |
& - kV*uVel(i,j,k,bi,bj) |
78 |
ENDIF |
79 |
ENDDO |
80 |
ENDDO |
81 |
|
82 |
RETURN |
83 |
END |
84 |
|
85 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
86 |
CBOP |
87 |
C !ROUTINE: APPLY_FORCING_V |
88 |
C !INTERFACE: |
89 |
SUBROUTINE APPLY_FORCING_V( |
90 |
U gV_arr, |
91 |
I iMin,iMax,jMin,jMax, k, bi, bj, |
92 |
I myTime, myIter, myThid ) |
93 |
C !DESCRIPTION: \bv |
94 |
C *==========================================================* |
95 |
C | S/R APPLY_FORCING_V |
96 |
C | o Contains problem specific forcing for merid velocity. |
97 |
C *==========================================================* |
98 |
C | Adds terms to gV for forcing by external sources |
99 |
C | e.g. wind stress, bottom friction etc ... |
100 |
C *==========================================================* |
101 |
C \ev |
102 |
|
103 |
C !USES: |
104 |
IMPLICIT NONE |
105 |
C == Global data == |
106 |
#include "SIZE.h" |
107 |
#include "EEPARAMS.h" |
108 |
#include "PARAMS.h" |
109 |
#include "GRID.h" |
110 |
#include "DYNVARS.h" |
111 |
#include "FFIELDS.h" |
112 |
|
113 |
C !INPUT/OUTPUT PARAMETERS: |
114 |
C gV_arr :: the tendency array |
115 |
C iMin,iMax :: Working range of x-index for applying forcing. |
116 |
C jMin,jMax :: Working range of y-index for applying forcing. |
117 |
C k :: Current vertical level index |
118 |
C bi,bj :: Current tile indices |
119 |
C myTime :: Current time in simulation |
120 |
C myIter :: Current iteration number |
121 |
C myThid :: my Thread Id number |
122 |
_RL gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
123 |
INTEGER iMin, iMax, jMin, jMax |
124 |
INTEGER k, bi, bj |
125 |
_RL myTime |
126 |
INTEGER myIter |
127 |
INTEGER myThid |
128 |
|
129 |
C !LOCAL VARIABLES: |
130 |
C i,j :: Loop counters |
131 |
INTEGER i, j |
132 |
CEOP |
133 |
_RL recip_P0g, termP, kV, kF, sigma_b |
134 |
|
135 |
C-- Forcing term |
136 |
kF = 1. _d 0/86400. _d 0 |
137 |
sigma_b = 0.7 _d 0 |
138 |
DO j=1,sNy+1 |
139 |
c DO i=1,sNx |
140 |
C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1] |
141 |
DO i=0,sNx+1 |
142 |
IF ( maskS(i,j,k,bi,bj).EQ.oneRS ) THEN |
143 |
recip_P0g = MAX(recip_Rcol(i,j,bi,bj),recip_Rcol(i,j-1,bi,bj)) |
144 |
termP = 0.5 _d 0*( MIN( rF(k)*recip_P0g, oneRL ) |
145 |
& +rF(k+1)*recip_P0g ) |
146 |
c termP = 0.5 _d 0*( rF(k) + rF(k+1) )*recip_P0g |
147 |
kV = kF*MAX( zeroRL, (termP-sigma_b)/(1. _d 0-sigma_b) ) |
148 |
gV_arr(i,j) = gV_arr(i,j) |
149 |
& - kV*vVel(i,j,k,bi,bj) |
150 |
ENDIF |
151 |
ENDDO |
152 |
ENDDO |
153 |
|
154 |
RETURN |
155 |
END |
156 |
|
157 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
158 |
CBOP |
159 |
C !ROUTINE: APPLY_FORCING_T |
160 |
C !INTERFACE: |
161 |
SUBROUTINE APPLY_FORCING_T( |
162 |
U gT_arr, |
163 |
I iMin,iMax,jMin,jMax, k, bi, bj, |
164 |
I myTime, myIter, myThid ) |
165 |
C !DESCRIPTION: \bv |
166 |
C *==========================================================* |
167 |
C | S/R APPLY_FORCING_T |
168 |
C | o Contains problem specific forcing for temperature. |
169 |
C *==========================================================* |
170 |
C | Adds terms to gT for forcing by external sources |
171 |
C | e.g. heat flux, climatalogical relaxation, etc ... |
172 |
C *==========================================================* |
173 |
C \ev |
174 |
|
175 |
C !USES: |
176 |
IMPLICIT NONE |
177 |
C == Global data == |
178 |
#include "SIZE.h" |
179 |
#include "EEPARAMS.h" |
180 |
#include "PARAMS.h" |
181 |
#include "GRID.h" |
182 |
#include "DYNVARS.h" |
183 |
#include "FFIELDS.h" |
184 |
#include "SURFACE.h" |
185 |
|
186 |
C !INPUT/OUTPUT PARAMETERS: |
187 |
C gT_arr :: the tendency array |
188 |
C iMin,iMax :: Working range of x-index for applying forcing. |
189 |
C jMin,jMax :: Working range of y-index for applying forcing. |
190 |
C k :: Current vertical level index |
191 |
C bi,bj :: Current tile indices |
192 |
C myTime :: Current time in simulation |
193 |
C myIter :: Current iteration number |
194 |
C myThid :: my Thread Id number |
195 |
_RL gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
196 |
INTEGER iMin, iMax, jMin, jMax |
197 |
INTEGER k, bi, bj |
198 |
_RL myTime |
199 |
INTEGER myIter |
200 |
INTEGER myThid |
201 |
|
202 |
C !LOCAL VARIABLES: |
203 |
C i,j :: Loop counters |
204 |
C kSurface :: index of surface level |
205 |
INTEGER i, j |
206 |
CEOP |
207 |
_RL thetaLim,kT,ka,ks,sigma_b,term1,term2,thetaEq,termP |
208 |
|
209 |
C-- Forcing term |
210 |
ka = 1. _d 0/(40. _d 0*86400. _d 0) |
211 |
ks = 1. _d 0/(4. _d 0 *86400. _d 0) |
212 |
sigma_b = 0.7 _d 0 |
213 |
DO j=1,sNy |
214 |
DO i=1,sNx |
215 |
term1 = 60. _d 0*(SIN(yC(i,j,bi,bj)*deg2rad)**2) |
216 |
termP = 0.5 _d 0*( rF(k) + rF(k+1) ) |
217 |
term2 = 10. _d 0*LOG(termP/atm_po) |
218 |
& *(COS(yC(i,j,bi,bj)*deg2rad)**2) |
219 |
thetaLim = 200. _d 0/ ((termP/atm_po)**atm_kappa) |
220 |
thetaEq = 315. _d 0-term1-term2 |
221 |
thetaEq = MAX(thetaLim,thetaEq) |
222 |
termP = 0.5 _d 0*( MIN(rF(k),Ro_surf(i,j,bi,bj)) |
223 |
& + rF(k+1) ) |
224 |
kT = ka+(ks-ka) |
225 |
& *MAX( zeroRL, |
226 |
& (termP*recip_Rcol(i,j,bi,bj)-sigma_b)/(1. _d 0-sigma_b) ) |
227 |
& *COS((yC(i,j,bi,bj)*deg2rad))**4 |
228 |
gT_arr(i,j) = gT_arr(i,j) |
229 |
& - kT*( theta(i,j,k,bi,bj)-thetaEq ) |
230 |
& *maskC(i,j,k,bi,bj) |
231 |
ENDDO |
232 |
ENDDO |
233 |
|
234 |
RETURN |
235 |
END |
236 |
|
237 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
238 |
CBOP |
239 |
C !ROUTINE: APPLY_FORCING_S |
240 |
C !INTERFACE: |
241 |
SUBROUTINE APPLY_FORCING_S( |
242 |
U gS_arr, |
243 |
I iMin,iMax,jMin,jMax, k, bi, bj, |
244 |
I myTime, myIter, myThid ) |
245 |
C !DESCRIPTION: \bv |
246 |
C *==========================================================* |
247 |
C | S/R APPLY_FORCING_S |
248 |
C | o Contains problem specific forcing for merid velocity. |
249 |
C *==========================================================* |
250 |
C | Adds terms to gS for forcing by external sources |
251 |
C | e.g. fresh-water flux, climatalogical relaxation, etc ... |
252 |
C *==========================================================* |
253 |
C \ev |
254 |
|
255 |
C !USES: |
256 |
IMPLICIT NONE |
257 |
C == Global data == |
258 |
#include "SIZE.h" |
259 |
#include "EEPARAMS.h" |
260 |
#include "PARAMS.h" |
261 |
#include "GRID.h" |
262 |
#include "DYNVARS.h" |
263 |
#include "FFIELDS.h" |
264 |
#include "SURFACE.h" |
265 |
|
266 |
C !INPUT/OUTPUT PARAMETERS: |
267 |
C gS_arr :: the tendency array |
268 |
C iMin,iMax :: Working range of x-index for applying forcing. |
269 |
C jMin,jMax :: Working range of y-index for applying forcing. |
270 |
C k :: Current vertical level index |
271 |
C bi,bj :: Current tile indices |
272 |
C myTime :: Current time in simulation |
273 |
C myIter :: Current iteration number |
274 |
C myThid :: my Thread Id number |
275 |
_RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
276 |
INTEGER iMin, iMax, jMin, jMax |
277 |
INTEGER k, bi, bj |
278 |
_RL myTime |
279 |
INTEGER myIter |
280 |
INTEGER myThid |
281 |
|
282 |
C !LOCAL VARIABLES: |
283 |
C i,j :: Loop counters |
284 |
c INTEGER i, j |
285 |
CEOP |
286 |
|
287 |
C-- Forcing term |
288 |
|
289 |
RETURN |
290 |
END |