/[MITgcm]/MITgcm/verification/rotating_tank/code/apply_forcing.F
ViewVC logotype

Contents of /MITgcm/verification/rotating_tank/code/apply_forcing.F

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


Revision 1.1 - (show annotations) (download)
Fri Jul 11 18:46:36 2014 UTC (7 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65, HEAD
 new file "apply_forcing.F" containing all the code previously in
  external_forcing.F

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 "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6
7 C-- File apply_forcing.F:
8 C-- Contents
9 C-- o APPLY_FORCING_U
10 C-- o APPLY_FORCING_V
11 C-- o APPLY_FORCING_T
12 C-- o APPLY_FORCING_S
13
14 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
15 CBOP
16 C !ROUTINE: APPLY_FORCING_U
17 C !INTERFACE:
18 SUBROUTINE APPLY_FORCING_U(
19 U gU_arr,
20 I iMin,iMax,jMin,jMax, k, bi, bj,
21 I myTime, myIter, myThid )
22 C !DESCRIPTION: \bv
23 C *==========================================================*
24 C | S/R APPLY_FORCING_U
25 C | o Contains problem specific forcing for zonal velocity.
26 C *==========================================================*
27 C | Adds terms to gU for forcing by external sources
28 C | e.g. wind stress, bottom friction etc ...
29 C *==========================================================*
30 C \ev
31
32 C !USES:
33 IMPLICIT NONE
34 C == Global data ==
35 #include "SIZE.h"
36 #include "EEPARAMS.h"
37 #include "PARAMS.h"
38 #include "GRID.h"
39 #include "DYNVARS.h"
40 #include "FFIELDS.h"
41
42 C !INPUT/OUTPUT PARAMETERS:
43 C gU_arr :: the tendency array
44 C iMin,iMax :: Working range of x-index for applying forcing.
45 C jMin,jMax :: Working range of y-index for applying forcing.
46 C k :: Current vertical level index
47 C bi,bj :: Current tile indices
48 C myTime :: Current time in simulation
49 C myIter :: Current iteration number
50 C myThid :: my Thread Id number
51 _RL gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 INTEGER iMin, iMax, jMin, jMax
53 INTEGER k, bi, bj
54 _RL myTime
55 INTEGER myIter
56 INTEGER myThid
57
58 C !LOCAL VARIABLES:
59 C i,j :: Loop counters
60 C kSurface :: index of surface level
61 INTEGER i, j
62 #ifdef USE_OLD_EXTERNAL_FORCING
63 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RL tmpVar
65 #else
66 INTEGER kSurface
67 #endif
68 CEOP
69
70 #ifdef USE_OLD_EXTERNAL_FORCING
71
72 DO j=1-OLy,sNy+OLy
73 DO i=1-OLx,sNx+OLx
74 locVar(i,j) = gU(i,j,k,bi,bj)
75 ENDDO
76 ENDDO
77 CALL EXTERNAL_FORCING_U(
78 I iMin, iMax, jMin, jMax, bi, bj, k,
79 I myTime, myThid )
80 DO j=1-OLy,sNy+OLy
81 DO i=1-OLx,sNx+OLx
82 tmpVar = gU(i,j,k,bi,bj) - locVar(i,j)
83 gU(i,j,k,bi,bj) = locVar(i,j)
84 gU_arr(i,j) = gU_arr(i,j) + tmpVar
85 ENDDO
86 ENDDO
87
88 #else /* USE_OLD_EXTERNAL_FORCING */
89
90 IF ( fluidIsAir ) THEN
91 kSurface = 0
92 ELSEIF ( usingPCoords ) THEN
93 kSurface = Nr
94 ELSE
95 kSurface = 1
96 ENDIF
97
98 C-- Forcing term
99 #ifdef ALLOW_AIM
100 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
101 U gU_arr,
102 I iMin,iMax,jMin,jMax, k, bi,bj,
103 I myTime, 0, myThid )
104 #endif /* ALLOW_AIM */
105
106 #ifdef ALLOW_ATM_PHYS
107 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
108 U gU_arr,
109 I iMin,iMax,jMin,jMax, k, bi,bj,
110 I myTime, 0, myThid )
111 #endif /* ALLOW_ATM_PHYS */
112
113 #ifdef ALLOW_FIZHI
114 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
115 U gU_arr,
116 I iMin,iMax,jMin,jMax, k, bi,bj,
117 I myTime, 0, myThid )
118 #endif /* ALLOW_FIZHI */
119
120 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
121 IF ( k .EQ. kSurface ) THEN
122 c DO j=1,sNy
123 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
124 DO j=0,sNy+1
125 DO i=1,sNx+1
126 gU_arr(i,j) = gU_arr(i,j)
127 & +foFacMom*surfaceForcingU(i,j,bi,bj)
128 & *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
129 ENDDO
130 ENDDO
131 ELSEIF ( kSurface.EQ.-1 ) THEN
132 DO j=0,sNy+1
133 DO i=1,sNx+1
134 IF ( kSurfW(i,j,bi,bj).EQ.k ) THEN
135 gU_arr(i,j) = gU_arr(i,j)
136 & +foFacMom*surfaceForcingU(i,j,bi,bj)
137 & *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
138 ENDIF
139 ENDDO
140 ENDDO
141 ENDIF
142
143 #ifdef ALLOW_EDDYPSI
144 CALL TAUEDDY_TENDENCY_APPLY_U(
145 U gU_arr,
146 I iMin,iMax,jMin,jMax, k, bi,bj,
147 I myTime, 0, myThid )
148 #endif
149
150 #ifdef ALLOW_RBCS
151 IF (useRBCS) THEN
152 CALL RBCS_ADD_TENDENCY(
153 U gU_arr,
154 I k, bi, bj, -1,
155 I myTime, 0, myThid )
156
157 ENDIF
158 #endif /* ALLOW_RBCS */
159
160 #ifdef ALLOW_OBCS
161 IF (useOBCS) THEN
162 CALL OBCS_SPONGE_U(
163 U gU_arr,
164 I iMin,iMax,jMin,jMax, k, bi,bj,
165 I myTime, 0, myThid )
166 ENDIF
167 #endif /* ALLOW_OBCS */
168
169 #ifdef ALLOW_MYPACKAGE
170 IF ( useMYPACKAGE ) THEN
171 CALL MYPACKAGE_TENDENCY_APPLY_U(
172 U gU_arr,
173 I iMin,iMax,jMin,jMax, k, bi,bj,
174 I myTime, 0, myThid )
175 ENDIF
176 #endif /* ALLOW_MYPACKAGE */
177
178 #endif /* USE_OLD_EXTERNAL_FORCING */
179
180 RETURN
181 END
182
183 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
184 CBOP
185 C !ROUTINE: APPLY_FORCING_V
186 C !INTERFACE:
187 SUBROUTINE APPLY_FORCING_V(
188 U gV_arr,
189 I iMin,iMax,jMin,jMax, k, bi, bj,
190 I myTime, myIter, myThid )
191 C !DESCRIPTION: \bv
192 C *==========================================================*
193 C | S/R APPLY_FORCING_V
194 C | o Contains problem specific forcing for merid velocity.
195 C *==========================================================*
196 C | Adds terms to gV for forcing by external sources
197 C | e.g. wind stress, bottom friction etc ...
198 C *==========================================================*
199 C \ev
200
201 C !USES:
202 IMPLICIT NONE
203 C == Global data ==
204 #include "SIZE.h"
205 #include "EEPARAMS.h"
206 #include "PARAMS.h"
207 #include "GRID.h"
208 #include "DYNVARS.h"
209 #include "FFIELDS.h"
210
211 C !INPUT/OUTPUT PARAMETERS:
212 C gV_arr :: the tendency array
213 C iMin,iMax :: Working range of x-index for applying forcing.
214 C jMin,jMax :: Working range of y-index for applying forcing.
215 C k :: Current vertical level index
216 C bi,bj :: Current tile indices
217 C myTime :: Current time in simulation
218 C myIter :: Current iteration number
219 C myThid :: my Thread Id number
220 _RL gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
221 INTEGER iMin, iMax, jMin, jMax
222 INTEGER k, bi, bj
223 _RL myTime
224 INTEGER myIter
225 INTEGER myThid
226
227 C !LOCAL VARIABLES:
228 C i,j :: Loop counters
229 C kSurface :: index of surface level
230 INTEGER i, j
231 #ifdef USE_OLD_EXTERNAL_FORCING
232 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
233 _RL tmpVar
234 #else
235 INTEGER kSurface
236 #endif
237 CEOP
238
239 #ifdef USE_OLD_EXTERNAL_FORCING
240
241 DO j=1-OLy,sNy+OLy
242 DO i=1-OLx,sNx+OLx
243 locVar(i,j) = gV(i,j,k,bi,bj)
244 ENDDO
245 ENDDO
246 CALL EXTERNAL_FORCING_V(
247 I iMin, iMax, jMin, jMax, bi, bj, k,
248 I myTime, myThid )
249 DO j=1-OLy,sNy+OLy
250 DO i=1-OLx,sNx+OLx
251 tmpVar = gV(i,j,k,bi,bj) - locVar(i,j)
252 gV(i,j,k,bi,bj) = locVar(i,j)
253 gV_arr(i,j) = gV_arr(i,j) + tmpVar
254 ENDDO
255 ENDDO
256
257 #else /* USE_OLD_EXTERNAL_FORCING */
258
259 IF ( fluidIsAir ) THEN
260 kSurface = 0
261 ELSEIF ( usingPCoords ) THEN
262 kSurface = Nr
263 ELSE
264 kSurface = 1
265 ENDIF
266
267 C-- Forcing term
268 #ifdef ALLOW_AIM
269 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
270 U gV_arr,
271 I iMin,iMax,jMin,jMax, k, bi,bj,
272 I myTime, 0, myThid )
273 #endif /* ALLOW_AIM */
274
275 #ifdef ALLOW_ATM_PHYS
276 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
277 U gV_arr,
278 I iMin,iMax,jMin,jMax, k, bi,bj,
279 I myTime, 0, myThid )
280 #endif /* ALLOW_ATM_PHYS */
281
282 #ifdef ALLOW_FIZHI
283 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
284 U gV_arr,
285 I iMin,iMax,jMin,jMax, k, bi,bj,
286 I myTime, 0, myThid )
287 #endif /* ALLOW_FIZHI */
288
289 C Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
290 IF ( k .EQ. kSurface ) THEN
291 DO j=1,sNy+1
292 c DO i=1,sNx
293 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
294 DO i=0,sNx+1
295 gV_arr(i,j) = gV_arr(i,j)
296 & +foFacMom*surfaceForcingV(i,j,bi,bj)
297 & *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
298 ENDDO
299 ENDDO
300 ELSEIF ( kSurface.EQ.-1 ) THEN
301 DO j=1,sNy+1
302 DO i=0,sNx+1
303 IF ( kSurfS(i,j,bi,bj).EQ.k ) THEN
304 gV_arr(i,j) = gV_arr(i,j)
305 & +foFacMom*surfaceForcingV(i,j,bi,bj)
306 & *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
307 ENDIF
308 ENDDO
309 ENDDO
310 ENDIF
311
312 #ifdef ALLOW_EDDYPSI
313 CALL TAUEDDY_TENDENCY_APPLY_V(
314 U gV_arr,
315 I iMin,iMax,jMin,jMax, k, bi,bj,
316 I myTime, 0, myThid )
317 #endif
318
319 #ifdef ALLOW_RBCS
320 IF (useRBCS) THEN
321 CALL RBCS_ADD_TENDENCY(
322 U gV_arr,
323 I k, bi, bj, -2,
324 I myTime, 0, myThid )
325 ENDIF
326 #endif /* ALLOW_RBCS */
327
328 #ifdef ALLOW_OBCS
329 IF (useOBCS) THEN
330 CALL OBCS_SPONGE_V(
331 U gV_arr,
332 I iMin,iMax,jMin,jMax, k, bi,bj,
333 I myTime, 0, myThid )
334 ENDIF
335 #endif /* ALLOW_OBCS */
336
337 #ifdef ALLOW_MYPACKAGE
338 IF ( useMYPACKAGE ) THEN
339 CALL MYPACKAGE_TENDENCY_APPLY_V(
340 U gV_arr,
341 I iMin,iMax,jMin,jMax, k, bi,bj,
342 I myTime, 0, myThid )
343 ENDIF
344 #endif /* ALLOW_MYPACKAGE */
345
346 #endif /* USE_OLD_EXTERNAL_FORCING */
347
348 RETURN
349 END
350
351 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
352 CBOP
353 C !ROUTINE: APPLY_FORCING_T
354 C !INTERFACE:
355 SUBROUTINE APPLY_FORCING_T(
356 U gT_arr,
357 I iMin,iMax,jMin,jMax, k, bi, bj,
358 I myTime, myIter, myThid )
359 C !DESCRIPTION: \bv
360 C *==========================================================*
361 C | S/R APPLY_FORCING_T
362 C | o Contains problem specific forcing for temperature.
363 C *==========================================================*
364 C | Adds terms to gT for forcing by external sources
365 C | e.g. heat flux, climatalogical relaxation, etc ...
366 C *==========================================================*
367 C \ev
368
369 C !USES:
370 IMPLICIT NONE
371 C == Global data ==
372 #include "SIZE.h"
373 #include "EEPARAMS.h"
374 #include "PARAMS.h"
375 #include "GRID.h"
376 #include "DYNVARS.h"
377 #include "FFIELDS.h"
378 #include "SURFACE.h"
379
380 C !INPUT/OUTPUT PARAMETERS:
381 C gT_arr :: the tendency array
382 C iMin,iMax :: Working range of x-index for applying forcing.
383 C jMin,jMax :: Working range of y-index for applying forcing.
384 C k :: Current vertical level index
385 C bi,bj :: Current tile indices
386 C myTime :: Current time in simulation
387 C myIter :: Current iteration number
388 C myThid :: my Thread Id number
389 _RL gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
390 INTEGER iMin, iMax, jMin, jMax
391 INTEGER k, bi, bj
392 _RL myTime
393 INTEGER myIter
394 INTEGER myThid
395
396 C !LOCAL VARIABLES:
397 C i,j :: Loop counters
398 C kSurface :: index of surface level
399 INTEGER i, j
400 #ifdef USE_OLD_EXTERNAL_FORCING
401 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
402 _RL tmpVar
403 #else
404 INTEGER kSurface
405 INTEGER km, kc, kp
406 _RL tmpVar(1:sNx,1:sNy)
407 _RL tmpFac, delPI
408 _RL recip_Cp
409 CEOP
410 #ifdef SHORTWAVE_HEATING
411 _RL minusone
412 PARAMETER (minusOne=-1.)
413 _RL swfracb(2)
414 INTEGER kp1
415 #endif
416 C---- Experiment specific declaration: starts here -------------------
417 C iG, jG :: Global index temps.
418 C hC, hW, hE, hN, hS :: Fractional vertical distance open to fluid temps.
419 C dFlux[WENS] :: Diffusive flux normal to each cell face.
420 C faceArea :: Temp. for holding area normal to tempurature gradient.
421 INTEGER iG, jG
422 _RL hC, hW, hE, hN, hS
423 _RL dFluxW, dFluxE, dFluxN, dFluxS
424 _RL faceArea
425
426 C-- Forcing term
427 C Add term which represents the diffusive flux from a circular cylinder of temperature tCylIm in the
428 C interior of the simulation domain. Result is a tendency which is determined from the finite
429 C volume formulated divergence of the diffusive heat flux due to the local cylinder
430 C temperature, fluid temperature difference.
431 C kDiffCyl :: Diffusion coefficient
432 C tCylIn :: Temperature of the inner boundary of the cylinder
433 C tCylOut :: Temperature of the outer boundary cylinder
434 C iGSource :: Index space I (global) coordinate for source center.
435 C jGSource :: Index space J (global) coordinate for source center.
436 C rSource :: Extent of the source term region. Loop will skip checking points outside
437 C :: this region. Within this region the source heating will be added
438 C :: to any points that are at a land - fluid boundary. rSource is in grid
439 C :: points, so that points checked are ophi(iGlobal,jGlobal) such that
440 C :: iGlobal=iGsource +/- rSource, jGlobal = jGsource +/- rSource.
441 C :: rSource, iGSource and jGSource only need to define a quadrilateral that
442 C :: includes the cylinder and no other land, they do not need to be exact.
443 C afe:
444 C the below appears to be an attempt to make the heat flux somewhat general in regards
445 C to bathymetry, but jmc pointed out some ways it could be better. this is not
446 C an issue at this point (July 04) since all experiments are being done with straight-
447 C sided tank and cyclinder walls.
448 C some todo items:
449 C * add terms to top and bottom -- probably critical!!!
450 C * make tCyl more flexible -- maybe have as separate inner and outer variables
451 C or, eventually, a forcing field
452 C * think about possible problems with differing heat diffusion rates in wall materials
453 C (plexiglass, air, water in the compound tank case)
454 _RL kDiffCyl
455 _RL tCyl
456 INTEGER rSource
457 INTEGER iGSource
458 INTEGER jGSource
459 C---- Experiment specific declaration: ends here -------------------
460 #endif
461
462 #ifdef USE_OLD_EXTERNAL_FORCING
463
464 DO j=1-OLy,sNy+OLy
465 DO i=1-OLx,sNx+OLx
466 locVar(i,j) = gT(i,j,k,bi,bj)
467 ENDDO
468 ENDDO
469 CALL EXTERNAL_FORCING_T(
470 I iMin, iMax, jMin, jMax, bi, bj, k,
471 I myTime, myThid )
472 DO j=1-OLy,sNy+OLy
473 DO i=1-OLx,sNx+OLx
474 tmpVar = gT(i,j,k,bi,bj) - locVar(i,j)
475 gT(i,j,k,bi,bj) = locVar(i,j)
476 gT_arr(i,j) = gT_arr(i,j) + tmpVar
477 ENDDO
478 ENDDO
479
480 #else /* USE_OLD_EXTERNAL_FORCING */
481
482 IF ( fluidIsAir ) THEN
483 kSurface = 0
484 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
485 kSurface = -1
486 ELSEIF ( usingPCoords ) THEN
487 kSurface = Nr
488 ELSE
489 kSurface = 1
490 ENDIF
491 recip_Cp = 1. _d 0 / HeatCapacity_Cp
492
493 C-- Forcing term
494 #ifdef ALLOW_AIM
495 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
496 U gT_arr,
497 I iMin,iMax,jMin,jMax, k, bi,bj,
498 I myTime, 0, myThid )
499 #endif /* ALLOW_AIM */
500
501 #ifdef ALLOW_ATM_PHYS
502 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
503 U gT_arr,
504 I iMin,iMax,jMin,jMax, k, bi,bj,
505 I myTime, 0, myThid )
506 #endif /* ALLOW_ATM_PHYS */
507
508 #ifdef ALLOW_FIZHI
509 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
510 U gT_arr,
511 I iMin,iMax,jMin,jMax, k, bi,bj,
512 I myTime, 0, myThid )
513 #endif /* ALLOW_FIZHI */
514
515 #ifdef ALLOW_ADDFLUID
516 IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
517 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
518 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
519 DO j=1,sNy
520 DO i=1,sNx
521 gT_arr(i,j) = gT_arr(i,j)
522 & + addMass(i,j,k,bi,bj)*mass2rUnit
523 & *( temp_addMass - theta(i,j,k,bi,bj) )
524 & *recip_rA(i,j,bi,bj)
525 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
526 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
527 ENDDO
528 ENDDO
529 ELSE
530 DO j=1,sNy
531 DO i=1,sNx
532 gT_arr(i,j) = gT_arr(i,j)
533 & + addMass(i,j,k,bi,bj)*mass2rUnit
534 & *( temp_addMass - tRef(k) )
535 & *recip_rA(i,j,bi,bj)
536 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
537 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
538 ENDDO
539 ENDDO
540 ENDIF
541 ENDIF
542 #endif /* ALLOW_ADDFLUID */
543
544 #ifdef ALLOW_FRICTION_HEATING
545 IF ( addFrictionHeating ) THEN
546 IF ( fluidIsAir ) THEN
547 C conversion from in-situ Temp to Pot.Temp
548 tmpFac = (atm_Po/rC(k))**atm_kappa
549 C conversion from W/m^2/r_unit to K/s
550 tmpFac = (tmpFac/atm_Cp) * mass2rUnit
551 ELSE
552 C conversion from W/m^2/r_unit to K/s
553 tmpFac = recip_Cp * mass2rUnit
554 ENDIF
555 DO j=1,sNy
556 DO i=1,sNx
557 gT_arr(i,j) = gT_arr(i,j)
558 & + frictionHeating(i,j,k,bi,bj)
559 & *tmpFac*recip_rA(i,j,bi,bj)
560 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
561 ENDDO
562 ENDDO
563 ENDIF
564 #endif /* ALLOW_FRICTION_HEATING */
565
566 IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
567 C-- Compressible fluid: account for difference between moist and dry air
568 C specific volume in Enthalpy equation (+ V.dP term), since only the
569 C dry air part is accounted for in the (dry) Pot.Temp formulation.
570 C Used centered averaging from interface to center (consistent with
571 C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
572 C as for Theta_v in CALC_PHI_HYD
573
574 C conversion from in-situ Temp to Pot.Temp
575 tmpFac = (atm_Po/rC(k))**atm_kappa
576 C conversion from W/kg to K/s
577 tmpFac = tmpFac/atm_Cp
578 km = k-1
579 kc = k
580 kp = k+1
581 IF ( k.EQ.1 ) THEN
582 DO j=1,sNy
583 DO i=1,sNx
584 tmpVar(i,j) = 0.
585 ENDDO
586 ENDDO
587 ELSE
588 delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
589 & - (rC(kc)/atm_Po)**atm_kappa )
590 DO j=1,sNy
591 DO i=1,sNx
592 tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
593 & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
594 & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
595 & )*maskC(i,j,km,bi,bj)*0.25 _d 0
596 ENDDO
597 ENDDO
598 ENDIF
599 IF ( k.LT.Nr ) THEN
600 delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
601 & - (rC(kp)/atm_Po)**atm_kappa )
602 DO j=1,sNy
603 DO i=1,sNx
604 tmpVar(i,j) = tmpVar(i,j)
605 & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
606 & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
607 & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
608 & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
609 ENDDO
610 ENDDO
611 ENDIF
612 DO j=1,sNy
613 DO i=1,sNx
614 gT_arr(i,j) = gT_arr(i,j)
615 & + tmpVar(i,j)*tmpFac
616 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
617 ENDDO
618 ENDDO
619 #ifdef ALLOW_DIAGNOSTICS
620 IF ( useDiagnostics ) THEN
621 C conversion to W/m^2
622 tmpFac = rUnit2mass
623 CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
624 & 'MoistCor', kc, 1, 3, bi,bj,myThid )
625 ENDIF
626 #endif /* ALLOW_DIAGNOSTICS */
627 ENDIF
628
629 C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
630 IF ( k .EQ. kSurface ) THEN
631 DO j=1,sNy
632 DO i=1,sNx
633 gT_arr(i,j) = gT_arr(i,j)
634 & +surfaceForcingT(i,j,bi,bj)
635 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
636 ENDDO
637 ENDDO
638 ELSEIF ( kSurface.EQ.-1 ) THEN
639 DO j=1,sNy
640 DO i=1,sNx
641 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
642 gT_arr(i,j) = gT_arr(i,j)
643 & +surfaceForcingT(i,j,bi,bj)
644 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
645 ENDIF
646 ENDDO
647 ENDDO
648 ENDIF
649
650 IF (linFSConserveTr) THEN
651 DO j=1,sNy
652 DO i=1,sNx
653 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
654 gT_arr(i,j) = gT_arr(i,j)
655 & +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
656 ENDIF
657 ENDDO
658 ENDDO
659 ENDIF
660
661 #ifdef SHORTWAVE_HEATING
662 C Penetrating SW radiation
663 c IF ( usePenetratingSW ) THEN
664 swfracb(1)=abs(rF(k))
665 swfracb(2)=abs(rF(k+1))
666 CALL SWFRAC(
667 I 2, minusOne,
668 U swfracb,
669 I myTime, 1, myThid )
670 kp1 = k+1
671 IF (k.EQ.Nr) THEN
672 kp1 = k
673 swfracb(2)=0. _d 0
674 ENDIF
675 DO j=1,sNy
676 DO i=1,sNx
677 gT_arr(i,j) = gT_arr(i,j)
678 & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
679 & -swfracb(2)*maskC(i,j,kp1, bi,bj))
680 & *recip_Cp*mass2rUnit
681 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
682 ENDDO
683 ENDDO
684 c ENDIF
685 #endif
686
687 C---- Experiment specific code: starts here --------------------------
688 kDiffCyl = 3. _d -7
689 rSource = 3
690 iGSource = 30
691 jGSource = 8
692 DO j=1,sNy
693 DO i=1,sNx
694 dFluxW = 0.
695 dFluxE = 0.
696 dFluxN = 0.
697 dFluxS = 0.
698 jG = myYGlobalLo-1+(bj-1)*sNy+J
699 iG = myXGlobalLo-1+(bi-1)*sNx+I
700 c IF(jG .GE. jGSource-rSource .AND. jG .LE. jGSource+rSource) THEN
701 c the following bites the big one
702 IF(jG .LE. 10) THEN
703 tCyl = tCylIn
704 ELSE
705 tCyl = tCylOut
706 ENDIF
707 c IF(iG .GE. iGSource-rSource .AND. iG .LE. iGSource+rSource) THEN
708 hC = hFacC(i ,j ,k,bi,bj)
709 hW = hFacW(i ,j ,k,bi,bj)
710 hE = hFacW(i+1,j ,k,bi,bj)
711 hN = hFacS(i ,j+1,k,bi,bj)
712 hS = hFacS(i ,j ,k,bi,bj)
713 IF ( hC .NE. 0. .AND .hW .EQ. 0. ) THEN
714 C Cylinder to west
715 faceArea = drF(k)*dyG(i,j,bi,bj)
716 dFluxW =
717 & -faceArea*kDiffCyl*(theta(i,j,k,bi,bj) - tCyl)
718 & *recip_dxC(i,j,bi,bj)
719 ENDIF
720 IF ( hC .NE. 0. .AND .hE .EQ. 0. ) THEN
721 C Cylinder to east
722 faceArea = drF(k)*dyG(i+1,j,bi,bj)
723 dFluxE =
724 & -faceArea*kDiffCyl*(tCyl - theta(i,j,k,bi,bj))
725 & *recip_dxC(i,j,bi,bj)
726 ENDIF
727 IF ( hC .NE. 0. .AND .hN .EQ. 0. ) THEN
728 C Cylinder to north
729 faceArea = drF(k)*dxG(i,j+1,bi,bj)
730 dFluxN =
731 & -faceArea*kDiffCyl*(tCyl-theta(i,j,k,bi,bj))
732 & *recip_dyC(i,j,bi,bj)
733 ENDIF
734 IF ( hC .NE. 0. .AND .hS .EQ. 0. ) THEN
735 C Cylinder to south
736 faceArea = drF(k)*dxG(i,j,bi,bj)
737 dFluxS =
738 & -faceArea*kDiffCyl*(theta(i,j,k,bi,bj) - tCyl)
739 & *recip_dyC(i,j,bi,bj)
740 ENDIF
741 c ENDIF
742 c ENDIF
743 C Net tendency term is minus flux divergence divided by the volume.
744 gT_arr(i,j) = gT_arr(i,j)
745 & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
746 & *recip_rA(i,j,bi,bj)
747 & *(
748 & dFluxE-dFluxW
749 & +dFluxN-dFluxS
750 & )
751
752 ENDDO
753 ENDDO
754 C---- Experiment specific code: ends here --------------------------
755
756 #ifdef ALLOW_FRAZIL
757 IF ( useFRAZIL )
758 & CALL FRAZIL_TENDENCY_APPLY_T(
759 U gT_arr,
760 I iMin,iMax,jMin,jMax, k, bi,bj,
761 I myTime, 0, myThid )
762 #endif /* ALLOW_FRAZIL */
763
764 #ifdef ALLOW_SHELFICE
765 IF ( useShelfIce )
766 & CALL SHELFICE_FORCING_T(
767 U gT_arr,
768 I iMin,iMax,jMin,jMax, k, bi,bj,
769 I myTime, 0, myThid )
770 #endif /* ALLOW_SHELFICE */
771
772 #ifdef ALLOW_ICEFRONT
773 IF ( useICEFRONT )
774 & CALL ICEFRONT_TENDENCY_APPLY_T(
775 U gT_arr,
776 I k, bi, bj, myTime, 0, myThid )
777 #endif /* ALLOW_ICEFRONT */
778
779 #ifdef ALLOW_SALT_PLUME
780 IF ( useSALT_PLUME )
781 & CALL SALT_PLUME_TENDENCY_APPLY_T(
782 U gT_arr,
783 I iMin,iMax,jMin,jMax, k, bi,bj,
784 I myTime, 0, myThid )
785 #endif /* ALLOW_SALT_PLUME */
786
787 #ifdef ALLOW_RBCS
788 IF (useRBCS) THEN
789 CALL RBCS_ADD_TENDENCY(
790 U gT_arr,
791 I k, bi, bj, 1,
792 I myTime, 0, myThid )
793 ENDIF
794 #endif /* ALLOW_RBCS */
795
796 #ifdef ALLOW_OBCS
797 IF (useOBCS) THEN
798 CALL OBCS_SPONGE_T(
799 U gT_arr,
800 I iMin,iMax,jMin,jMax, k, bi,bj,
801 I myTime, 0, myThid )
802 ENDIF
803 #endif /* ALLOW_OBCS */
804
805 #ifdef ALLOW_BBL
806 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
807 U gT_arr,
808 I iMin,iMax,jMin,jMax, k, bi,bj,
809 I myTime, 0, myThid )
810 #endif /* ALLOW_BBL */
811
812 #ifdef ALLOW_MYPACKAGE
813 IF ( useMYPACKAGE ) THEN
814 CALL MYPACKAGE_TENDENCY_APPLY_T(
815 U gT_arr,
816 I iMin,iMax,jMin,jMax, k, bi,bj,
817 I myTime, 0, myThid )
818 ENDIF
819 #endif /* ALLOW_MYPACKAGE */
820
821 #endif /* USE_OLD_EXTERNAL_FORCING */
822
823 RETURN
824 END
825
826 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
827 CBOP
828 C !ROUTINE: APPLY_FORCING_S
829 C !INTERFACE:
830 SUBROUTINE APPLY_FORCING_S(
831 U gS_arr,
832 I iMin,iMax,jMin,jMax, k, bi, bj,
833 I myTime, myIter, myThid )
834 C !DESCRIPTION: \bv
835 C *==========================================================*
836 C | S/R APPLY_FORCING_S
837 C | o Contains problem specific forcing for merid velocity.
838 C *==========================================================*
839 C | Adds terms to gS for forcing by external sources
840 C | e.g. fresh-water flux, climatalogical relaxation, etc ...
841 C *==========================================================*
842 C \ev
843
844 C !USES:
845 IMPLICIT NONE
846 C == Global data ==
847 #include "SIZE.h"
848 #include "EEPARAMS.h"
849 #include "PARAMS.h"
850 #include "GRID.h"
851 #include "DYNVARS.h"
852 #include "FFIELDS.h"
853 #include "SURFACE.h"
854
855 C !INPUT/OUTPUT PARAMETERS:
856 C gS_arr :: the tendency array
857 C iMin,iMax :: Working range of x-index for applying forcing.
858 C jMin,jMax :: Working range of y-index for applying forcing.
859 C k :: Current vertical level index
860 C bi,bj :: Current tile indices
861 C myTime :: Current time in simulation
862 C myIter :: Current iteration number
863 C myThid :: my Thread Id number
864 _RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
865 INTEGER iMin, iMax, jMin, jMax
866 INTEGER k, bi, bj
867 _RL myTime
868 INTEGER myIter
869 INTEGER myThid
870
871 C !LOCAL VARIABLES:
872 C i,j :: Loop counters
873 C kSurface :: index of surface level
874 INTEGER i, j
875 #ifdef USE_OLD_EXTERNAL_FORCING
876 _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
877 _RL tmpVar
878 #else
879 INTEGER kSurface
880 #endif
881 CEOP
882
883 #ifdef USE_OLD_EXTERNAL_FORCING
884
885 DO j=1-OLy,sNy+OLy
886 DO i=1-OLx,sNx+OLx
887 locVar(i,j) = gS(i,j,k,bi,bj)
888 ENDDO
889 ENDDO
890 CALL EXTERNAL_FORCING_S(
891 I iMin, iMax, jMin, jMax, bi, bj, k,
892 I myTime, myThid )
893 DO j=1-OLy,sNy+OLy
894 DO i=1-OLx,sNx+OLx
895 tmpVar = gS(i,j,k,bi,bj) - locVar(i,j)
896 gS(i,j,k,bi,bj) = locVar(i,j)
897 gS_arr(i,j) = gS_arr(i,j) + tmpVar
898 ENDDO
899 ENDDO
900
901 #else /* USE_OLD_EXTERNAL_FORCING */
902
903 IF ( fluidIsAir ) THEN
904 kSurface = 0
905 ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
906 kSurface = -1
907 ELSEIF ( usingPCoords ) THEN
908 kSurface = Nr
909 ELSE
910 kSurface = 1
911 ENDIF
912
913 C-- Forcing term
914 #ifdef ALLOW_AIM
915 IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
916 U gS_arr,
917 I iMin,iMax,jMin,jMax, k, bi,bj,
918 I myTime, 0, myThid )
919 #endif /* ALLOW_AIM */
920
921 #ifdef ALLOW_ATM_PHYS
922 IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
923 U gS_arr,
924 I iMin,iMax,jMin,jMax, k, bi,bj,
925 I myTime, 0, myThid )
926 #endif /* ALLOW_ATM_PHYS */
927
928 #ifdef ALLOW_FIZHI
929 IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
930 U gS_arr,
931 I iMin,iMax,jMin,jMax, k, bi,bj,
932 I myTime, 0, myThid )
933 #endif /* ALLOW_FIZHI */
934
935 #ifdef ALLOW_ADDFLUID
936 IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
937 IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
938 & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
939 DO j=1,sNy
940 DO i=1,sNx
941 gS_arr(i,j) = gS_arr(i,j)
942 & + addMass(i,j,k,bi,bj)*mass2rUnit
943 & *( salt_addMass - salt(i,j,k,bi,bj) )
944 & *recip_rA(i,j,bi,bj)
945 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
946 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
947 ENDDO
948 ENDDO
949 ELSE
950 DO j=1,sNy
951 DO i=1,sNx
952 gS_arr(i,j) = gS_arr(i,j)
953 & + addMass(i,j,k,bi,bj)*mass2rUnit
954 & *( salt_addMass - sRef(k) )
955 & *recip_rA(i,j,bi,bj)
956 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
957 C & *recip_deepFac2C(k)*recip_rhoFacC(k)
958 ENDDO
959 ENDDO
960 ENDIF
961 ENDIF
962 #endif /* ALLOW_ADDFLUID */
963
964 C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
965 IF ( k .EQ. kSurface ) THEN
966 DO j=1,sNy
967 DO i=1,sNx
968 gS_arr(i,j) = gS_arr(i,j)
969 & +surfaceForcingS(i,j,bi,bj)
970 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
971 ENDDO
972 ENDDO
973 ELSEIF ( kSurface.EQ.-1 ) THEN
974 DO j=1,sNy
975 DO i=1,sNx
976 IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
977 gS_arr(i,j) = gS_arr(i,j)
978 & +surfaceForcingS(i,j,bi,bj)
979 & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
980 ENDIF
981 ENDDO
982 ENDDO
983 ENDIF
984
985 IF (linFSConserveTr) THEN
986 DO j=1,sNy
987 DO i=1,sNx
988 IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
989 gS_arr(i,j) = gS_arr(i,j)
990 & +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
991 ENDIF
992 ENDDO
993 ENDDO
994 ENDIF
995
996 #ifdef ALLOW_SHELFICE
997 IF ( useShelfIce )
998 & CALL SHELFICE_FORCING_S(
999 U gS_arr,
1000 I iMin,iMax,jMin,jMax, k, bi,bj,
1001 I myTime, 0, myThid )
1002 #endif /* ALLOW_SHELFICE */
1003
1004 #ifdef ALLOW_ICEFRONT
1005 IF ( useICEFRONT )
1006 & CALL ICEFRONT_TENDENCY_APPLY_S(
1007 U gS_arr,
1008 I k, bi, bj, myTime, 0, myThid )
1009 #endif /* ALLOW_ICEFRONT */
1010
1011 #ifdef ALLOW_SALT_PLUME
1012 IF ( useSALT_PLUME )
1013 & CALL SALT_PLUME_TENDENCY_APPLY_S(
1014 U gS_arr,
1015 I iMin,iMax,jMin,jMax, k, bi,bj,
1016 I myTime, 0, myThid )
1017 #endif /* ALLOW_SALT_PLUME */
1018
1019 #ifdef ALLOW_RBCS
1020 IF (useRBCS) THEN
1021 CALL RBCS_ADD_TENDENCY(
1022 U gS_arr,
1023 I k, bi, bj, 2,
1024 I myTime, 0, myThid )
1025 ENDIF
1026 #endif /* ALLOW_RBCS */
1027
1028 #ifdef ALLOW_OBCS
1029 IF (useOBCS) THEN
1030 CALL OBCS_SPONGE_S(
1031 U gS_arr,
1032 I iMin,iMax,jMin,jMax, k, bi,bj,
1033 I myTime, 0, myThid )
1034 ENDIF
1035 #endif /* ALLOW_OBCS */
1036
1037 #ifdef ALLOW_BBL
1038 IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
1039 U gS_arr,
1040 I iMin,iMax,jMin,jMax, k, bi,bj,
1041 I myTime, 0, myThid )
1042 #endif /* ALLOW_BBL */
1043
1044 #ifdef ALLOW_MYPACKAGE
1045 IF ( useMYPACKAGE ) THEN
1046 CALL MYPACKAGE_TENDENCY_APPLY_S(
1047 U gS_arr,
1048 I iMin,iMax,jMin,jMax, k, bi,bj,
1049 I myTime, 0, myThid )
1050 ENDIF
1051 #endif /* ALLOW_MYPACKAGE */
1052
1053 #endif /* USE_OLD_EXTERNAL_FORCING */
1054
1055 RETURN
1056 END

  ViewVC Help
Powered by ViewVC 1.1.22