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

Annotation of /MITgcm/model/src/apply_forcing.F

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

Revision 1.1 - (hide annotations) (download)
Fri Jul 11 18:43:41 2014 UTC (9 years, 10 months ago) by jmc
Branch: MAIN 
- new file "apply_forcing.F" containing all the code previously in
  external_forcing.F, but with new argument list: pass, as new argument,
  the current level tendency array to update (instead of a direct update
  of the common bloc array). Change the corresponding calls.
1 jmc 1.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     #endif
417    
418     #ifdef USE_OLD_EXTERNAL_FORCING
419    
420     DO j=1-OLy,sNy+OLy
421     DO i=1-OLx,sNx+OLx
422     locVar(i,j) = gT(i,j,k,bi,bj)
423     ENDDO
424     ENDDO
425     CALL EXTERNAL_FORCING_T(
426     I iMin, iMax, jMin, jMax, bi, bj, k,
427     I myTime, myThid )
428     DO j=1-OLy,sNy+OLy
429     DO i=1-OLx,sNx+OLx
430     tmpVar = gT(i,j,k,bi,bj) - locVar(i,j)
431     gT(i,j,k,bi,bj) = locVar(i,j)
432     gT_arr(i,j) = gT_arr(i,j) + tmpVar
433     ENDDO
434     ENDDO
435    
436     #else /* USE_OLD_EXTERNAL_FORCING */
437    
438     IF ( fluidIsAir ) THEN
439     kSurface = 0
440     ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
441     kSurface = -1
442     ELSEIF ( usingPCoords ) THEN
443     kSurface = Nr
444     ELSE
445     kSurface = 1
446     ENDIF
447     recip_Cp = 1. _d 0 / HeatCapacity_Cp
448    
449     C-- Forcing term
450     #ifdef ALLOW_AIM
451     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
452     U gT_arr,
453     I iMin,iMax,jMin,jMax, k, bi,bj,
454     I myTime, 0, myThid )
455     #endif /* ALLOW_AIM */
456    
457     #ifdef ALLOW_ATM_PHYS
458     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
459     U gT_arr,
460     I iMin,iMax,jMin,jMax, k, bi,bj,
461     I myTime, 0, myThid )
462     #endif /* ALLOW_ATM_PHYS */
463    
464     #ifdef ALLOW_FIZHI
465     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
466     U gT_arr,
467     I iMin,iMax,jMin,jMax, k, bi,bj,
468     I myTime, 0, myThid )
469     #endif /* ALLOW_FIZHI */
470    
471     #ifdef ALLOW_ADDFLUID
472     IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
473     IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
474     & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
475     DO j=1,sNy
476     DO i=1,sNx
477     gT_arr(i,j) = gT_arr(i,j)
478     & + addMass(i,j,k,bi,bj)*mass2rUnit
479     & *( temp_addMass - theta(i,j,k,bi,bj) )
480     & *recip_rA(i,j,bi,bj)
481     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
482     C & *recip_deepFac2C(k)*recip_rhoFacC(k)
483     ENDDO
484     ENDDO
485     ELSE
486     DO j=1,sNy
487     DO i=1,sNx
488     gT_arr(i,j) = gT_arr(i,j)
489     & + addMass(i,j,k,bi,bj)*mass2rUnit
490     & *( temp_addMass - tRef(k) )
491     & *recip_rA(i,j,bi,bj)
492     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
493     C & *recip_deepFac2C(k)*recip_rhoFacC(k)
494     ENDDO
495     ENDDO
496     ENDIF
497     ENDIF
498     #endif /* ALLOW_ADDFLUID */
499    
500     #ifdef ALLOW_FRICTION_HEATING
501     IF ( addFrictionHeating ) THEN
502     IF ( fluidIsAir ) THEN
503     C conversion from in-situ Temp to Pot.Temp
504     tmpFac = (atm_Po/rC(k))**atm_kappa
505     C conversion from W/m^2/r_unit to K/s
506     tmpFac = (tmpFac/atm_Cp) * mass2rUnit
507     ELSE
508     C conversion from W/m^2/r_unit to K/s
509     tmpFac = recip_Cp * mass2rUnit
510     ENDIF
511     DO j=1,sNy
512     DO i=1,sNx
513     gT_arr(i,j) = gT_arr(i,j)
514     & + frictionHeating(i,j,k,bi,bj)
515     & *tmpFac*recip_rA(i,j,bi,bj)
516     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
517     ENDDO
518     ENDDO
519     ENDIF
520     #endif /* ALLOW_FRICTION_HEATING */
521    
522     IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
523     C-- Compressible fluid: account for difference between moist and dry air
524     C specific volume in Enthalpy equation (+ V.dP term), since only the
525     C dry air part is accounted for in the (dry) Pot.Temp formulation.
526     C Used centered averaging from interface to center (consistent with
527     C conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
528     C as for Theta_v in CALC_PHI_HYD
529    
530     C conversion from in-situ Temp to Pot.Temp
531     tmpFac = (atm_Po/rC(k))**atm_kappa
532     C conversion from W/kg to K/s
533     tmpFac = tmpFac/atm_Cp
534     km = k-1
535     kc = k
536     kp = k+1
537     IF ( k.EQ.1 ) THEN
538     DO j=1,sNy
539     DO i=1,sNx
540     tmpVar(i,j) = 0.
541     ENDDO
542     ENDDO
543     ELSE
544     delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
545     & - (rC(kc)/atm_Po)**atm_kappa )
546     DO j=1,sNy
547     DO i=1,sNx
548     tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
549     & *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
550     & + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
551     & )*maskC(i,j,km,bi,bj)*0.25 _d 0
552     ENDDO
553     ENDDO
554     ENDIF
555     IF ( k.LT.Nr ) THEN
556     delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
557     & - (rC(kp)/atm_Po)**atm_kappa )
558     DO j=1,sNy
559     DO i=1,sNx
560     tmpVar(i,j) = tmpVar(i,j)
561     & + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
562     & *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
563     & + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
564     & )*maskC(i,j,kp,bi,bj)*0.25 _d 0
565     ENDDO
566     ENDDO
567     ENDIF
568     DO j=1,sNy
569     DO i=1,sNx
570     gT_arr(i,j) = gT_arr(i,j)
571     & + tmpVar(i,j)*tmpFac
572     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
573     ENDDO
574     ENDDO
575     #ifdef ALLOW_DIAGNOSTICS
576     IF ( useDiagnostics ) THEN
577     C conversion to W/m^2
578     tmpFac = rUnit2mass
579     CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
580     & 'MoistCor', kc, 1, 3, bi,bj,myThid )
581     ENDIF
582     #endif /* ALLOW_DIAGNOSTICS */
583     ENDIF
584    
585     C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
586     IF ( k .EQ. kSurface ) THEN
587     DO j=1,sNy
588     DO i=1,sNx
589     gT_arr(i,j) = gT_arr(i,j)
590     & +surfaceForcingT(i,j,bi,bj)
591     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
592     ENDDO
593     ENDDO
594     ELSEIF ( kSurface.EQ.-1 ) THEN
595     DO j=1,sNy
596     DO i=1,sNx
597     IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
598     gT_arr(i,j) = gT_arr(i,j)
599     & +surfaceForcingT(i,j,bi,bj)
600     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
601     ENDIF
602     ENDDO
603     ENDDO
604     ENDIF
605    
606     IF (linFSConserveTr) THEN
607     DO j=1,sNy
608     DO i=1,sNx
609     IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
610     gT_arr(i,j) = gT_arr(i,j)
611     & +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
612     ENDIF
613     ENDDO
614     ENDDO
615     ENDIF
616    
617     #ifdef SHORTWAVE_HEATING
618     C Penetrating SW radiation
619     c IF ( usePenetratingSW ) THEN
620     swfracb(1)=abs(rF(k))
621     swfracb(2)=abs(rF(k+1))
622     CALL SWFRAC(
623     I 2, minusOne,
624     U swfracb,
625     I myTime, 1, myThid )
626     kp1 = k+1
627     IF (k.EQ.Nr) THEN
628     kp1 = k
629     swfracb(2)=0. _d 0
630     ENDIF
631     DO j=1,sNy
632     DO i=1,sNx
633     gT_arr(i,j) = gT_arr(i,j)
634     & -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
635     & -swfracb(2)*maskC(i,j,kp1, bi,bj))
636     & *recip_Cp*mass2rUnit
637     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
638     ENDDO
639     ENDDO
640     c ENDIF
641     #endif
642    
643     #ifdef ALLOW_FRAZIL
644     IF ( useFRAZIL )
645     & CALL FRAZIL_TENDENCY_APPLY_T(
646     U gT_arr,
647     I iMin,iMax,jMin,jMax, k, bi,bj,
648     I myTime, 0, myThid )
649     #endif /* ALLOW_FRAZIL */
650    
651     #ifdef ALLOW_SHELFICE
652     IF ( useShelfIce )
653     & CALL SHELFICE_FORCING_T(
654     U gT_arr,
655     I iMin,iMax,jMin,jMax, k, bi,bj,
656     I myTime, 0, myThid )
657     #endif /* ALLOW_SHELFICE */
658    
659     #ifdef ALLOW_ICEFRONT
660     IF ( useICEFRONT )
661     & CALL ICEFRONT_TENDENCY_APPLY_T(
662     U gT_arr,
663     I k, bi, bj, myTime, 0, myThid )
664     #endif /* ALLOW_ICEFRONT */
665    
666     #ifdef ALLOW_SALT_PLUME
667     IF ( useSALT_PLUME )
668     & CALL SALT_PLUME_TENDENCY_APPLY_T(
669     U gT_arr,
670     I iMin,iMax,jMin,jMax, k, bi,bj,
671     I myTime, 0, myThid )
672     #endif /* ALLOW_SALT_PLUME */
673    
674     #ifdef ALLOW_RBCS
675     IF (useRBCS) THEN
676     CALL RBCS_ADD_TENDENCY(
677     U gT_arr,
678     I k, bi, bj, 1,
679     I myTime, 0, myThid )
680     ENDIF
681     #endif /* ALLOW_RBCS */
682    
683     #ifdef ALLOW_OBCS
684     IF (useOBCS) THEN
685     CALL OBCS_SPONGE_T(
686     U gT_arr,
687     I iMin,iMax,jMin,jMax, k, bi,bj,
688     I myTime, 0, myThid )
689     ENDIF
690     #endif /* ALLOW_OBCS */
691    
692     #ifdef ALLOW_BBL
693     IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
694     U gT_arr,
695     I iMin,iMax,jMin,jMax, k, bi,bj,
696     I myTime, 0, myThid )
697     #endif /* ALLOW_BBL */
698    
699     #ifdef ALLOW_MYPACKAGE
700     IF ( useMYPACKAGE ) THEN
701     CALL MYPACKAGE_TENDENCY_APPLY_T(
702     U gT_arr,
703     I iMin,iMax,jMin,jMax, k, bi,bj,
704     I myTime, 0, myThid )
705     ENDIF
706     #endif /* ALLOW_MYPACKAGE */
707    
708     #endif /* USE_OLD_EXTERNAL_FORCING */
709    
710     RETURN
711     END
712    
713     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
714     CBOP
715     C !ROUTINE: APPLY_FORCING_S
716     C !INTERFACE:
717     SUBROUTINE APPLY_FORCING_S(
718     U gS_arr,
719     I iMin,iMax,jMin,jMax, k, bi, bj,
720     I myTime, myIter, myThid )
721     C !DESCRIPTION: \bv
722     C *==========================================================*
723     C | S/R APPLY_FORCING_S
724     C | o Contains problem specific forcing for merid velocity.
725     C *==========================================================*
726     C | Adds terms to gS for forcing by external sources
727     C | e.g. fresh-water flux, climatalogical relaxation, etc ...
728     C *==========================================================*
729     C \ev
730    
731     C !USES:
732     IMPLICIT NONE
733     C == Global data ==
734     #include "SIZE.h"
735     #include "EEPARAMS.h"
736     #include "PARAMS.h"
737     #include "GRID.h"
738     #include "DYNVARS.h"
739     #include "FFIELDS.h"
740     #include "SURFACE.h"
741    
742     C !INPUT/OUTPUT PARAMETERS:
743     C gS_arr :: the tendency array
744     C iMin,iMax :: Working range of x-index for applying forcing.
745     C jMin,jMax :: Working range of y-index for applying forcing.
746     C k :: Current vertical level index
747     C bi,bj :: Current tile indices
748     C myTime :: Current time in simulation
749     C myIter :: Current iteration number
750     C myThid :: my Thread Id number
751     _RL gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
752     INTEGER iMin, iMax, jMin, jMax
753     INTEGER k, bi, bj
754     _RL myTime
755     INTEGER myIter
756     INTEGER myThid
757    
758     C !LOCAL VARIABLES:
759     C i,j :: Loop counters
760     C kSurface :: index of surface level
761     INTEGER i, j
762     #ifdef USE_OLD_EXTERNAL_FORCING
763     _RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
764     _RL tmpVar
765     #else
766     INTEGER kSurface
767     #endif
768     CEOP
769    
770     #ifdef USE_OLD_EXTERNAL_FORCING
771    
772     DO j=1-OLy,sNy+OLy
773     DO i=1-OLx,sNx+OLx
774     locVar(i,j) = gS(i,j,k,bi,bj)
775     ENDDO
776     ENDDO
777     CALL EXTERNAL_FORCING_S(
778     I iMin, iMax, jMin, jMax, bi, bj, k,
779     I myTime, myThid )
780     DO j=1-OLy,sNy+OLy
781     DO i=1-OLx,sNx+OLx
782     tmpVar = gS(i,j,k,bi,bj) - locVar(i,j)
783     gS(i,j,k,bi,bj) = locVar(i,j)
784     gS_arr(i,j) = gS_arr(i,j) + tmpVar
785     ENDDO
786     ENDDO
787    
788     #else /* USE_OLD_EXTERNAL_FORCING */
789    
790     IF ( fluidIsAir ) THEN
791     kSurface = 0
792     ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
793     kSurface = -1
794     ELSEIF ( usingPCoords ) THEN
795     kSurface = Nr
796     ELSE
797     kSurface = 1
798     ENDIF
799    
800     C-- Forcing term
801     #ifdef ALLOW_AIM
802     IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
803     U gS_arr,
804     I iMin,iMax,jMin,jMax, k, bi,bj,
805     I myTime, 0, myThid )
806     #endif /* ALLOW_AIM */
807    
808     #ifdef ALLOW_ATM_PHYS
809     IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
810     U gS_arr,
811     I iMin,iMax,jMin,jMax, k, bi,bj,
812     I myTime, 0, myThid )
813     #endif /* ALLOW_ATM_PHYS */
814    
815     #ifdef ALLOW_FIZHI
816     IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
817     U gS_arr,
818     I iMin,iMax,jMin,jMax, k, bi,bj,
819     I myTime, 0, myThid )
820     #endif /* ALLOW_FIZHI */
821    
822     #ifdef ALLOW_ADDFLUID
823     IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
824     IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
825     & .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
826     DO j=1,sNy
827     DO i=1,sNx
828     gS_arr(i,j) = gS_arr(i,j)
829     & + addMass(i,j,k,bi,bj)*mass2rUnit
830     & *( salt_addMass - salt(i,j,k,bi,bj) )
831     & *recip_rA(i,j,bi,bj)
832     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
833     C & *recip_deepFac2C(k)*recip_rhoFacC(k)
834     ENDDO
835     ENDDO
836     ELSE
837     DO j=1,sNy
838     DO i=1,sNx
839     gS_arr(i,j) = gS_arr(i,j)
840     & + addMass(i,j,k,bi,bj)*mass2rUnit
841     & *( salt_addMass - sRef(k) )
842     & *recip_rA(i,j,bi,bj)
843     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
844     C & *recip_deepFac2C(k)*recip_rhoFacC(k)
845     ENDDO
846     ENDDO
847     ENDIF
848     ENDIF
849     #endif /* ALLOW_ADDFLUID */
850    
851     C Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
852     IF ( k .EQ. kSurface ) THEN
853     DO j=1,sNy
854     DO i=1,sNx
855     gS_arr(i,j) = gS_arr(i,j)
856     & +surfaceForcingS(i,j,bi,bj)
857     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
858     ENDDO
859     ENDDO
860     ELSEIF ( kSurface.EQ.-1 ) THEN
861     DO j=1,sNy
862     DO i=1,sNx
863     IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
864     gS_arr(i,j) = gS_arr(i,j)
865     & +surfaceForcingS(i,j,bi,bj)
866     & *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
867     ENDIF
868     ENDDO
869     ENDDO
870     ENDIF
871    
872     IF (linFSConserveTr) THEN
873     DO j=1,sNy
874     DO i=1,sNx
875     IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
876     gS_arr(i,j) = gS_arr(i,j)
877     & +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
878     ENDIF
879     ENDDO
880     ENDDO
881     ENDIF
882    
883     #ifdef ALLOW_SHELFICE
884     IF ( useShelfIce )
885     & CALL SHELFICE_FORCING_S(
886     U gS_arr,
887     I iMin,iMax,jMin,jMax, k, bi,bj,
888     I myTime, 0, myThid )
889     #endif /* ALLOW_SHELFICE */
890    
891     #ifdef ALLOW_ICEFRONT
892     IF ( useICEFRONT )
893     & CALL ICEFRONT_TENDENCY_APPLY_S(
894     U gS_arr,
895     I k, bi, bj, myTime, 0, myThid )
896     #endif /* ALLOW_ICEFRONT */
897    
898     #ifdef ALLOW_SALT_PLUME
899     IF ( useSALT_PLUME )
900     & CALL SALT_PLUME_TENDENCY_APPLY_S(
901     U gS_arr,
902     I iMin,iMax,jMin,jMax, k, bi,bj,
903     I myTime, 0, myThid )
904     #endif /* ALLOW_SALT_PLUME */
905    
906     #ifdef ALLOW_RBCS
907     IF (useRBCS) THEN
908     CALL RBCS_ADD_TENDENCY(
909     U gS_arr,
910     I k, bi, bj, 2,
911     I myTime, 0, myThid )
912     ENDIF
913     #endif /* ALLOW_RBCS */
914    
915     #ifdef ALLOW_OBCS
916     IF (useOBCS) THEN
917     CALL OBCS_SPONGE_S(
918     U gS_arr,
919     I iMin,iMax,jMin,jMax, k, bi,bj,
920     I myTime, 0, myThid )
921     ENDIF
922     #endif /* ALLOW_OBCS */
923    
924     #ifdef ALLOW_BBL
925     IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
926     U gS_arr,
927     I iMin,iMax,jMin,jMax, k, bi,bj,
928     I myTime, 0, myThid )
929     #endif /* ALLOW_BBL */
930    
931     #ifdef ALLOW_MYPACKAGE
932     IF ( useMYPACKAGE ) THEN
933     CALL MYPACKAGE_TENDENCY_APPLY_S(
934     U gS_arr,
935     I iMin,iMax,jMin,jMax, k, bi,bj,
936     I myTime, 0, myThid )
937     ENDIF
938     #endif /* ALLOW_MYPACKAGE */
939    
940     #endif /* USE_OLD_EXTERNAL_FORCING */
941    
942     RETURN
943     END
  ViewVC Help
Powered by ViewVC 1.1.22