/[MITgcm]/MITgcm/model/src/external_forcing.F
ViewVC logotype

Contents of /MITgcm/model/src/external_forcing.F

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

Revision 1.72 - (show annotations) (download)
Wed Jan 21 14:35:16 2015 UTC (9 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.71: +2 -3 lines 
change units of frictionHeating field (from W to W/m^2)
1 C $Header: /u/gcmpack/MITgcm/model/src/external_forcing.F,v 1.71 2014/08/15 19:19:23 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 C-- File external_forcing.F:
8 C-- Contents
9 C-- o EXTERNAL_FORCING_U
10 C-- o EXTERNAL_FORCING_V
11 C-- o EXTERNAL_FORCING_T
12 C-- o EXTERNAL_FORCING_S
13
14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
15 CBOP
16 C !ROUTINE: EXTERNAL_FORCING_U
17 C !INTERFACE:
18 SUBROUTINE EXTERNAL_FORCING_U(
19 I iMin,iMax, jMin,jMax, bi,bj, kLev,
20 I myTime, myThid )
21 C !DESCRIPTION: \bv
22 C *==========================================================*
23 C | S/R EXTERNAL_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 == Routine arguments ==
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 bi,bj :: Current tile indices
46 C kLev :: Current vertical level index
47 C myTime :: Current time in simulation
48 C myThid :: Thread Id number
49 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
50 _RL myTime
51 INTEGER myThid
52
53 #ifdef USE_OLD_EXTERNAL_FORCING
54 C !LOCAL VARIABLES:
55 C == Local variables ==
56 C i,j :: Loop counters
57 C kSurface :: index of surface level
58 INTEGER i, j
59 INTEGER kSurface
60 CEOP
61
62 IF ( fluidIsAir ) THEN
63 kSurface = 0
64 ELSEIF ( usingPCoords ) THEN
65 kSurface = Nr
66 ELSE
67 kSurface = 1
68 ENDIF
69
70 C-- Forcing term
71 #ifdef ALLOW_AIM
72 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
73 U gU(1-OLx,1-OLy,kLev,bi,bj),
74 I iMin,iMax,jMin,jMax, kLev, bi,bj,
75 I myTime, 0, myThid )
76 #endif /* ALLOW_AIM */
77
78 #ifdef ALLOW_ATM_PHYS
79 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
80 U gU(1-OLx,1-OLy,kLev,bi,bj),
81 I iMin,iMax,jMin,jMax, kLev, bi,bj,
82 I myTime, 0, myThid )
83 #endif /* ALLOW_ATM_PHYS */
84
85 #ifdef ALLOW_FIZHI
86 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
87 U gU(1-OLx,1-OLy,kLev,bi,bj),
88 I iMin,iMax,jMin,jMax, kLev, bi,bj,
89 I myTime, 0, myThid )
90 #endif /* ALLOW_FIZHI */
91
92 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
93 IF ( kLev .EQ. kSurface ) THEN
94 c DO j=1,sNy
95 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
96 DO j=0,sNy+1
97 DO i=1,sNx+1
98 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
99 & +foFacMom*surfaceForcingU(i,j,bi,bj)
100 & *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
101 ENDDO
102 ENDDO
103 ELSEIF ( kSurface.EQ.-1 ) THEN
104 DO j=0,sNy+1
105 DO i=1,sNx+1
106 IF ( kSurfW(i,j,bi,bj).EQ.kLev ) THEN
107 gU(i,j,kLev,bi,bj) = gU(i,j,kLev,bi,bj)
108 & +foFacMom*surfaceForcingU(i,j,bi,bj)
109 & *recip_drF(kLev)*_recip_hFacW(i,j,kLev,bi,bj)
110 ENDIF
111 ENDDO
112 ENDDO
113 ENDIF
114
115 #ifdef ALLOW_EDDYPSI
116 CALL TAUEDDY_TENDENCY_APPLY_U(
117 U gU(1-OLx,1-OLy,kLev,bi,bj),
118 I iMin,iMax,jMin,jMax, kLev, bi,bj,
119 I myTime, 0, myThid )
120 #endif
121
122 #ifdef ALLOW_RBCS
123 IF (useRBCS) THEN
124 CALL RBCS_ADD_TENDENCY(
125 U gU(1-OLx,1-OLy,kLev,bi,bj),
126 I kLev, bi, bj, -1,
127 I myTime, 0, myThid )
128
129 ENDIF
130 #endif /* ALLOW_RBCS */
131
132 #ifdef ALLOW_OBCS
133 IF (useOBCS) THEN
134 CALL OBCS_SPONGE_U(
135 U gU(1-OLx,1-OLy,kLev,bi,bj),
136 I iMin,iMax,jMin,jMax, kLev, bi,bj,
137 I myTime, 0, myThid )
138 ENDIF
139 #endif /* ALLOW_OBCS */
140
141 #ifdef ALLOW_MYPACKAGE
142 IF ( useMYPACKAGE ) THEN
143 CALL MYPACKAGE_TENDENCY_APPLY_U(
144 U gU(1-OLx,1-OLy,kLev,bi,bj),
145 I iMin,iMax,jMin,jMax, kLev, bi,bj,
146 I myTime, 0, myThid )
147 ENDIF
148 #endif /* ALLOW_MYPACKAGE */
149
150 #endif /* USE_OLD_EXTERNAL_FORCING */
151 RETURN
152 END
153
154 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
155 CBOP
156 C !ROUTINE: EXTERNAL_FORCING_V
157 C !INTERFACE:
158 SUBROUTINE EXTERNAL_FORCING_V(
159 I iMin,iMax, jMin,jMax, bi,bj, kLev,
160 I myTime, myThid )
161 C !DESCRIPTION: \bv
162 C *==========================================================*
163 C | S/R EXTERNAL_FORCING_V
164 C | o Contains problem specific forcing for merid velocity.
165 C *==========================================================*
166 C | Adds terms to gV for forcing by external sources
167 C | e.g. wind stress, bottom friction etc ...
168 C *==========================================================*
169 C \ev
170
171 C !USES:
172 IMPLICIT NONE
173 C == Global data ==
174 #include "SIZE.h"
175 #include "EEPARAMS.h"
176 #include "PARAMS.h"
177 #include "GRID.h"
178 #include "DYNVARS.h"
179 #include "FFIELDS.h"
180
181 C !INPUT/OUTPUT PARAMETERS:
182 C == Routine arguments ==
183 C iMin,iMax :: Working range of x-index for applying forcing.
184 C jMin,jMax :: Working range of y-index for applying forcing.
185 C bi,bj :: Current tile indices
186 C kLev :: Current vertical level index
187 C myTime :: Current time in simulation
188 C myThid :: Thread Id number
189 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
190 _RL myTime
191 INTEGER myThid
192
193 #ifdef USE_OLD_EXTERNAL_FORCING
194 C !LOCAL VARIABLES:
195 C == Local variables ==
196 C i,j :: Loop counters
197 C kSurface :: index of surface level
198 INTEGER i, j
199 INTEGER kSurface
200 CEOP
201
202 IF ( fluidIsAir ) THEN
203 kSurface = 0
204 ELSEIF ( usingPCoords ) THEN
205 kSurface = Nr
206 ELSE
207 kSurface = 1
208 ENDIF
209
210 C-- Forcing term
211 #ifdef ALLOW_AIM
212 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
213 U gV(1-OLx,1-OLy,kLev,bi,bj),
214 I iMin,iMax,jMin,jMax, kLev, bi,bj,
215 I myTime, 0, myThid )
216 #endif /* ALLOW_AIM */
217
218 #ifdef ALLOW_ATM_PHYS
219 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
220 U gV(1-OLx,1-OLy,kLev,bi,bj),
221 I iMin,iMax,jMin,jMax, kLev, bi,bj,
222 I myTime, 0, myThid )
223 #endif /* ALLOW_ATM_PHYS */
224
225 #ifdef ALLOW_FIZHI
226 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
227 U gV(1-OLx,1-OLy,kLev,bi,bj),
228 I iMin,iMax,jMin,jMax, kLev, bi,bj,
229 I myTime, 0, myThid )
230 #endif /* ALLOW_FIZHI */
231
232 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
233 IF ( kLev .EQ. kSurface ) THEN
234 DO j=1,sNy+1
235 c DO i=1,sNx
236 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
237 DO i=0,sNx+1
238 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
239 & +foFacMom*surfaceForcingV(i,j,bi,bj)
240 & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
241 ENDDO
242 ENDDO
243 ELSEIF ( kSurface.EQ.-1 ) THEN
244 DO j=1,sNy+1
245 DO i=0,sNx+1
246 IF ( kSurfS(i,j,bi,bj).EQ.kLev ) THEN
247 gV(i,j,kLev,bi,bj) = gV(i,j,kLev,bi,bj)
248 & +foFacMom*surfaceForcingV(i,j,bi,bj)
249 & *recip_drF(kLev)*_recip_hFacS(i,j,kLev,bi,bj)
250 ENDIF
251 ENDDO
252 ENDDO
253 ENDIF
254
255 #ifdef ALLOW_EDDYPSI
256 CALL TAUEDDY_TENDENCY_APPLY_V(
257 U gV(1-OLx,1-OLy,kLev,bi,bj),
258 I iMin,iMax,jMin,jMax, kLev, bi,bj,
259 I myTime, 0, myThid )
260 #endif
261
262 #ifdef ALLOW_RBCS
263 IF (useRBCS) THEN
264 CALL RBCS_ADD_TENDENCY(
265 U gV(1-OLx,1-OLy,kLev,bi,bj),
266 I kLev, bi, bj, -2,
267 I myTime, 0, myThid )
268 ENDIF
269 #endif /* ALLOW_RBCS */
270
271 #ifdef ALLOW_OBCS
272 IF (useOBCS) THEN
273 CALL OBCS_SPONGE_V(
274 U gV(1-OLx,1-OLy,kLev,bi,bj),
275 I iMin,iMax,jMin,jMax, kLev, bi,bj,
276 I myTime, 0, myThid )
277 ENDIF
278 #endif /* ALLOW_OBCS */
279
280 #ifdef ALLOW_MYPACKAGE
281 IF ( useMYPACKAGE ) THEN
282 CALL MYPACKAGE_TENDENCY_APPLY_V(
283 U gV(1-OLx,1-OLy,kLev,bi,bj),
284 I iMin,iMax,jMin,jMax, kLev, bi,bj,
285 I myTime, 0, myThid )
286 ENDIF
287 #endif /* ALLOW_MYPACKAGE */
288
289 #endif /* USE_OLD_EXTERNAL_FORCING */
290 RETURN
291 END
292
293 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
294 CBOP
295 C !ROUTINE: EXTERNAL_FORCING_T
296 C !INTERFACE:
297 SUBROUTINE EXTERNAL_FORCING_T(
298 I iMin,iMax, jMin,jMax, bi,bj, kLev,
299 I myTime, myThid )
300 C !DESCRIPTION: \bv
301 C *==========================================================*
302 C | S/R EXTERNAL_FORCING_T
303 C | o Contains problem specific forcing for temperature.
304 C *==========================================================*
305 C | Adds terms to gT for forcing by external sources
306 C | e.g. heat flux, climatalogical relaxation, etc ...
307 C *==========================================================*
308 C \ev
309
310 C !USES:
311 IMPLICIT NONE
312 C == Global data ==
313 #include "SIZE.h"
314 #include "EEPARAMS.h"
315 #include "PARAMS.h"
316 #include "GRID.h"
317 #include "DYNVARS.h"
318 #include "FFIELDS.h"
319 #include "SURFACE.h"
320
321 C !INPUT/OUTPUT PARAMETERS:
322 C == Routine arguments ==
323 C iMin,iMax :: Working range of x-index for applying forcing.
324 C jMin,jMax :: Working range of y-index for applying forcing.
325 C bi,bj :: Current tile indices
326 C kLev :: Current vertical level index
327 C myTime :: Current time in simulation
328 C myThid :: Thread Id number
329 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
330 _RL myTime
331 INTEGER myThid
332
333 #ifdef USE_OLD_EXTERNAL_FORCING
334 C !LOCAL VARIABLES:
335 C == Local variables ==
336 C i,j :: Loop counters
337 C kSurface :: index of surface level
338 INTEGER i, j
339 INTEGER kSurface
340 INTEGER km, kc, kp
341 _RL tmpVar(1:sNx,1:sNy)
342 _RL tmpFac, delPI
343 _RL recip_Cp
344 CEOP
345 #ifdef SHORTWAVE_HEATING
346 _RL minusone
347 PARAMETER (minusOne=-1.)
348 _RL swfracb(2)
349 INTEGER kp1
350 #endif
351
352 IF ( fluidIsAir ) THEN
353 kSurface = 0
354 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
355 kSurface = -1
356 ELSEIF ( usingPCoords ) THEN
357 kSurface = Nr
358 ELSE
359 kSurface = 1
360 ENDIF
361 recip_Cp = 1. _d 0 / HeatCapacity_Cp
362
363 C-- Forcing term
364 #ifdef ALLOW_AIM
365 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
366 U gT(1-OLx,1-OLy,kLev,bi,bj),
367 I iMin,iMax,jMin,jMax, kLev, bi,bj,
368 I myTime, 0, myThid )
369 #endif /* ALLOW_AIM */
370
371 #ifdef ALLOW_ATM_PHYS
372 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
373 U gT(1-OLx,1-OLy,kLev,bi,bj),
374 I iMin,iMax,jMin,jMax, kLev, bi,bj,
375 I myTime, 0, myThid )
376 #endif /* ALLOW_ATM_PHYS */
377
378 #ifdef ALLOW_FIZHI
379 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
380 U gT(1-OLx,1-OLy,kLev,bi,bj),
381 I iMin,iMax,jMin,jMax, kLev, bi,bj,
382 I myTime, 0, myThid )
383 #endif /* ALLOW_FIZHI */
384
385 #ifdef ALLOW_ADDFLUID
386 IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
387 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
388 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
389 DO j=1,sNy
390 DO i=1,sNx
391 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
392 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
393 & *( temp_addMass - theta(i,j,kLev,bi,bj) )
394 & *recip_rA(i,j,bi,bj)
395 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
396 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
397 ENDDO
398 ENDDO
399 ELSE
400 DO j=1,sNy
401 DO i=1,sNx
402 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
403 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
404 & *( temp_addMass - tRef(kLev) )
405 & *recip_rA(i,j,bi,bj)
406 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
407 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
408 ENDDO
409 ENDDO
410 ENDIF
411 ENDIF
412 #endif /* ALLOW_ADDFLUID */
413
414 #ifdef ALLOW_FRICTION_HEATING
415 IF ( addFrictionHeating ) THEN
416 IF ( fluidIsAir ) THEN
417 C conversion from in-situ Temp to Pot.Temp
418 tmpFac = (atm_Po/rC(kLev))**atm_kappa
419 C conversion from W/m^2/r_unit to K/s
420 tmpFac = (tmpFac/atm_Cp) * mass2rUnit
421 ELSE
422 C conversion from W/m^2/r_unit to K/s
423 tmpFac = recip_Cp * mass2rUnit
424 ENDIF
425 DO j=1,sNy
426 DO i=1,sNx
427 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
428 & + frictionHeating(i,j,k,bi,bj)*tmpFac
429 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
430 ENDDO
431 ENDDO
432 ENDIF
433 #endif /* ALLOW_FRICTION_HEATING */
434
435 IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
436 C-- Compressible fluid: account for difference between moist and dry air
437 C specific volume in Enthalpy equation (+ V.dP term), since only the
438 C dry air part is accounted for in the (dry) Pot.Temp formulation.
439 C Used centered averaging from interface to center (consistent with
440 C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
441 C as for Theta_v in CALC_PHI_HYD
442
443 C conversion from in-situ Temp to Pot.Temp
444 tmpFac = (atm_Po/rC(kLev))**atm_kappa
445 C conversion from W/kg to K/s
446 tmpFac = tmpFac/atm_Cp
447 km = kLev-1
448 kc = kLev
449 kp = kLev+1
450 IF ( kLev.EQ.1 ) THEN
451 DO j=1,sNy
452 DO i=1,sNx
453 tmpVar(i,j) = 0.
454 ENDDO
455 ENDDO
456 ELSE
457 delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
458 & - (rC(kc)/atm_Po)**atm_kappa )
459 DO j=1,sNy
460 DO i=1,sNx
461 tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
462 & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
463 & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
464 & )*maskC(i,j,km,bi,bj)*0.25 _d 0
465 ENDDO
466 ENDDO
467 ENDIF
468 IF ( kLev.LT.Nr ) THEN
469 delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
470 & - (rC(kp)/atm_Po)**atm_kappa )
471 DO j=1,sNy
472 DO i=1,sNx
473 tmpVar(i,j) = tmpVar(i,j)
474 & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
475 & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
476 & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
477 & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
478 ENDDO
479 ENDDO
480 ENDIF
481 DO j=1,sNy
482 DO i=1,sNx
483 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
484 & + tmpVar(i,j)*tmpFac
485 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
486 ENDDO
487 ENDDO
488 #ifdef ALLOW_DIAGNOSTICS
489 IF ( useDiagnostics ) THEN
490 C conversion to W/m^2
491 tmpFac = rUnit2mass
492 CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
493 & 'MoistCor', kc, 1, 3, bi,bj,myThid )
494 ENDIF
495 #endif /* ALLOW_DIAGNOSTICS */
496 ENDIF
497
498 C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
499 IF ( kLev .EQ. kSurface ) THEN
500 DO j=1,sNy
501 DO i=1,sNx
502 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
503 & +surfaceForcingT(i,j,bi,bj)
504 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
505 ENDDO
506 ENDDO
507 ELSEIF ( kSurface.EQ.-1 ) THEN
508 DO j=1,sNy
509 DO i=1,sNx
510 IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
511 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
512 & +surfaceForcingT(i,j,bi,bj)
513 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
514 ENDIF
515 ENDDO
516 ENDDO
517 ENDIF
518
519 IF (linFSConserveTr) THEN
520 DO j=1,sNy
521 DO i=1,sNx
522 IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
523 gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
524 & +TsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
525 ENDIF
526 ENDDO
527 ENDDO
528 ENDIF
529
530 #ifdef SHORTWAVE_HEATING
531 C Penetrating SW radiation
532 c IF ( usePenetratingSW ) THEN
533 swfracb(1)=abs(rF(kLev))
534 swfracb(2)=abs(rF(kLev+1))
535 CALL SWFRAC(
536 I 2, minusOne,
537 U swfracb,
538 I myTime, 1, myThid )
539 kp1 = kLev+1
540 IF (kLev.EQ.Nr) THEN
541 kp1 = kLev
542 swfracb(2)=0. _d 0
543 ENDIF
544 DO j=1,sNy
545 DO i=1,sNx
546 gT(i,j,kLev,bi,bj) = gT(i,j,kLev,bi,bj)
547 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,kLev,bi,bj)
548 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
549 & *recip_Cp*mass2rUnit
550 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
551 ENDDO
552 ENDDO
553 c ENDIF
554 #endif
555
556 #ifdef ALLOW_FRAZIL
557 IF ( useFRAZIL )
558 & CALL FRAZIL_TENDENCY_APPLY_T(
559 U gT(1-OLx,1-OLy,kLev,bi,bj),
560 I iMin,iMax,jMin,jMax, kLev, bi,bj,
561 I myTime, 0, myThid )
562 #endif /* ALLOW_FRAZIL */
563
564 #ifdef ALLOW_SHELFICE
565 IF ( useShelfIce )
566 & CALL SHELFICE_FORCING_T(
567 U gT(1-OLx,1-OLy,kLev,bi,bj),
568 I iMin,iMax,jMin,jMax, kLev, bi,bj,
569 I myTime, 0, myThid )
570 #endif /* ALLOW_SHELFICE */
571
572 #ifdef ALLOW_ICEFRONT
573 IF ( useICEFRONT )
574 & CALL ICEFRONT_TENDENCY_APPLY_T(
575 U gT(1-OLx,1-OLy,kLev,bi,bj),
576 I kLev, bi, bj, myTime, 0, myThid )
577 #endif /* ALLOW_ICEFRONT */
578
579 #ifdef ALLOW_SALT_PLUME
580 IF ( useSALT_PLUME )
581 & CALL SALT_PLUME_TENDENCY_APPLY_T(
582 U gT(1-OLx,1-OLy,kLev,bi,bj),
583 I iMin,iMax,jMin,jMax, kLev, bi,bj,
584 I myTime, 0, myThid )
585 #endif /* ALLOW_SALT_PLUME */
586
587 #ifdef ALLOW_RBCS
588 IF (useRBCS) THEN
589 CALL RBCS_ADD_TENDENCY(
590 U gT(1-OLx,1-OLy,kLev,bi,bj),
591 I kLev, bi, bj, 1,
592 I myTime, 0, myThid )
593 ENDIF
594 #endif /* ALLOW_RBCS */
595
596 #ifdef ALLOW_OBCS
597 IF (useOBCS) THEN
598 CALL OBCS_SPONGE_T(
599 U gT(1-OLx,1-OLy,kLev,bi,bj),
600 I iMin,iMax,jMin,jMax, kLev, bi,bj,
601 I myTime, 0, myThid )
602 ENDIF
603 #endif /* ALLOW_OBCS */
604
605 #ifdef ALLOW_BBL
606 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
607 U gT(1-OLx,1-OLy,kLev,bi,bj),
608 I iMin,iMax,jMin,jMax, kLev, bi,bj,
609 I myTime, 0, myThid )
610 #endif /* ALLOW_BBL */
611
612 #ifdef ALLOW_MYPACKAGE
613 IF ( useMYPACKAGE ) THEN
614 CALL MYPACKAGE_TENDENCY_APPLY_T(
615 U gT(1-OLx,1-OLy,kLev,bi,bj),
616 I iMin,iMax,jMin,jMax, kLev, bi,bj,
617 I myTime, 0, myThid )
618 ENDIF
619 #endif /* ALLOW_MYPACKAGE */
620
621 #endif /* USE_OLD_EXTERNAL_FORCING */
622 RETURN
623 END
624
625 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
626 CBOP
627 C !ROUTINE: EXTERNAL_FORCING_S
628 C !INTERFACE:
629 SUBROUTINE EXTERNAL_FORCING_S(
630 I iMin,iMax, jMin,jMax, bi,bj, kLev,
631 I myTime, myThid )
632
633 C !DESCRIPTION: \bv
634 C *==========================================================*
635 C | S/R EXTERNAL_FORCING_S
636 C | o Contains problem specific forcing for merid velocity.
637 C *==========================================================*
638 C | Adds terms to gS for forcing by external sources
639 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
640 C *==========================================================*
641 C \ev
642
643 C !USES:
644 IMPLICIT NONE
645 C == Global data ==
646 #include "SIZE.h"
647 #include "EEPARAMS.h"
648 #include "PARAMS.h"
649 #include "GRID.h"
650 #include "DYNVARS.h"
651 #include "FFIELDS.h"
652 #include "SURFACE.h"
653
654 C !INPUT/OUTPUT PARAMETERS:
655 C == Routine arguments ==
656 C iMin,iMax :: Working range of x-index for applying forcing.
657 C jMin,jMax :: Working range of y-index for applying forcing.
658 C bi,bj :: Current tile indices
659 C kLev :: Current vertical level index
660 C myTime :: Current time in simulation
661 C myThid :: Thread Id number
662 INTEGER iMin, iMax, jMin, jMax, kLev, bi, bj
663 _RL myTime
664 INTEGER myThid
665
666 #ifdef USE_OLD_EXTERNAL_FORCING
667 C !LOCAL VARIABLES:
668 C == Local variables ==
669 C i,j :: Loop counters
670 C kSurface :: index of surface level
671 INTEGER i, j
672 INTEGER kSurface
673 CEOP
674
675 IF ( fluidIsAir ) THEN
676 kSurface = 0
677 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
678 kSurface = -1
679 ELSEIF ( usingPCoords ) THEN
680 kSurface = Nr
681 ELSE
682 kSurface = 1
683 ENDIF
684
685 C-- Forcing term
686 #ifdef ALLOW_AIM
687 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
688 U gS(1-OLx,1-OLy,kLev,bi,bj),
689 I iMin,iMax,jMin,jMax, kLev, bi,bj,
690 I myTime, 0, myThid )
691 #endif /* ALLOW_AIM */
692
693 #ifdef ALLOW_ATM_PHYS
694 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
695 U gS(1-OLx,1-OLy,kLev,bi,bj),
696 I iMin,iMax,jMin,jMax, kLev, bi,bj,
697 I myTime, 0, myThid )
698 #endif /* ALLOW_ATM_PHYS */
699
700 #ifdef ALLOW_FIZHI
701 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
702 U gS(1-OLx,1-OLy,kLev,bi,bj),
703 I iMin,iMax,jMin,jMax, kLev, bi,bj,
704 I myTime, 0, myThid )
705 #endif /* ALLOW_FIZHI */
706
707 #ifdef ALLOW_ADDFLUID
708 IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
709 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
710 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
711 DO j=1,sNy
712 DO i=1,sNx
713 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
714 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
715 & *( salt_addMass - salt(i,j,kLev,bi,bj) )
716 & *recip_rA(i,j,bi,bj)
717 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
718 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
719 ENDDO
720 ENDDO
721 ELSE
722 DO j=1,sNy
723 DO i=1,sNx
724 gS(i,j,kLev,bi,bj) = gS(i,j,kLev,bi,bj)
725 & + addMass(i,j,kLev,bi,bj)*mass2rUnit
726 & *( salt_addMass - sRef(kLev) )
727 & *recip_rA(i,j,bi,bj)
728 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
729 C & *recip_deepFac2C(kLev)*recip_rhoFacC(kLev)
730 ENDDO
731 ENDDO
732 ENDIF
733 ENDIF
734 #endif /* ALLOW_ADDFLUID */
735
736 C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
737 IF ( kLev .EQ. kSurface ) THEN
738 DO j=1,sNy
739 DO i=1,sNx
740 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
741 & +surfaceForcingS(i,j,bi,bj)
742 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
743 ENDDO
744 ENDDO
745 ELSEIF ( kSurface.EQ.-1 ) THEN
746 DO j=1,sNy
747 DO i=1,sNx
748 IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
749 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
750 & +surfaceForcingS(i,j,bi,bj)
751 & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
752 ENDIF
753 ENDDO
754 ENDDO
755 ENDIF
756
757 IF (linFSConserveTr) THEN
758 DO j=1,sNy
759 DO i=1,sNx
760 IF (kLev .EQ. kSurfC(i,j,bi,bj)) THEN
761 gS(i,j,kLev,bi,bj)=gS(i,j,kLev,bi,bj)
762 & +SsurfCor*recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
763 ENDIF
764 ENDDO
765 ENDDO
766 ENDIF
767
768 #ifdef ALLOW_SHELFICE
769 IF ( useShelfIce )
770 & CALL SHELFICE_FORCING_S(
771 U gS(1-OLx,1-OLy,kLev,bi,bj),
772 I iMin,iMax,jMin,jMax, kLev, bi,bj,
773 I myTime, 0, myThid )
774 #endif /* ALLOW_SHELFICE */
775
776 #ifdef ALLOW_ICEFRONT
777 IF ( useICEFRONT )
778 & CALL ICEFRONT_TENDENCY_APPLY_S(
779 U gS(1-OLx,1-OLy,kLev,bi,bj),
780 I kLev, bi, bj, myTime, 0, myThid )
781 #endif /* ALLOW_ICEFRONT */
782
783 #ifdef ALLOW_SALT_PLUME
784 IF ( useSALT_PLUME )
785 & CALL SALT_PLUME_TENDENCY_APPLY_S(
786 U gS(1-OLx,1-OLy,kLev,bi,bj),
787 I iMin,iMax,jMin,jMax, kLev, bi,bj,
788 I myTime, 0, myThid )
789 #endif /* ALLOW_SALT_PLUME */
790
791 #ifdef ALLOW_RBCS
792 IF (useRBCS) THEN
793 CALL RBCS_ADD_TENDENCY(
794 U gS(1-OLx,1-OLy,kLev,bi,bj),
795 I kLev, bi, bj, 2,
796 I myTime, 0, myThid )
797 ENDIF
798 #endif /* ALLOW_RBCS */
799
800 #ifdef ALLOW_OBCS
801 IF (useOBCS) THEN
802 CALL OBCS_SPONGE_S(
803 U gS(1-OLx,1-OLy,kLev,bi,bj),
804 I iMin,iMax,jMin,jMax, kLev, bi,bj,
805 I myTime, 0, myThid )
806 ENDIF
807 #endif /* ALLOW_OBCS */
808
809 #ifdef ALLOW_BBL
810 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
811 U gS(1-OLx,1-OLy,kLev,bi,bj),
812 I iMin,iMax,jMin,jMax, kLev, bi,bj,
813 I myTime, 0, myThid )
814 #endif /* ALLOW_BBL */
815
816 #ifdef ALLOW_MYPACKAGE
817 IF ( useMYPACKAGE ) THEN
818 CALL MYPACKAGE_TENDENCY_APPLY_S(
819 U gS(1-OLx,1-OLy,kLev,bi,bj),
820 I iMin,iMax,jMin,jMax, kLev, bi,bj,
821 I myTime, 0, myThid )
822 ENDIF
823 #endif /* ALLOW_MYPACKAGE */
824
825 #endif /* USE_OLD_EXTERNAL_FORCING */
826 RETURN
827 END
  ViewVC Help
Powered by ViewVC 1.1.22